Matlab(3):连续系统梁的基础理论及MATLAB编程3.1基础理论xy)(xP将原本分布在梁上的载荷转换成各节点上的集中载荷和弯矩,从而使问题简化因为梁的轴向无作用力,因此忽略轴向的力的影响iyFjyFiMjM3.1.1将元素载荷化为节点载荷首先将整个梁假设为一个长度为L的元素,再利用形函数将作用于整个元素上的载荷转为分别作用于节点上的集中力和弯矩,其值可由下式确定:{}∫=LTdxxP0])[(NF其中:(1){}F)(xP][N——为已经转换为作用于节点的集中力及弯矩;——为原本分布于梁元素上的载荷函数,为标量函数;——为该梁元素的形函数。14}{×∈⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=RMFMFjjyiiyF41][×∈RNiyFjyFiMjMxy3.1.2应变能由能量法可知梁的应变能为:∫=LdxEIMU0221⎟⎟⎠⎞⎜⎜⎝⎛∂∂=22xuEIM其中:(2)(3)——弹性模量;——横截面的惯性矩;——横截面的弯矩;——为距离原点x的点在t时刻的y方向的变形量。EIMuiyFjyFiMjMxydxuxPWL∫=0)(∫⎟⎟⎠⎞⎜⎜⎝⎛∂∂=LdxxuEIU022221将(3)式代入(2)式,可得应变能为:3.1.3外力所做的功外力所做的功可表示为:3.1.4总势能梁的总势能为应变能减去外力所做的功,可表示为:(4)(5)WUV−=∫⎥⎥⎦⎤⎢⎢⎣⎡−⎟⎟⎠⎞⎜⎜⎝⎛∂∂=LdxPuxuEIV022221将(4)、(5)式代入(6)式,有:(6)(7))}()]{([),(txNtxuφ=将表示为如下形式:u][Nux其中,为形函数,是坐标的函数,因此在静力学对微分时只对起作用。为时间的函数,因此在动力学中对t微分时,只对起作用。][Nxu}{φ}{φ(8){}14×∈⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=RMFMFFjjyiiy{}14×∈⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=Ruujjiiθθφ分别为梁的两个端点的y方向的变形和两个截面的转角jjiiuuθθ、、、(分离变量法)iyFjyFiMjMxy∫⎥⎥⎦⎤⎢⎢⎣⎡−⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂=LTTdxNPxNxNEIV02222}]{[}{}{21φφφ}{][φxNxu∂∂=∂∂tNtu∂∂=∂∂}{][φTTTNu][}{φ=}]{[][}{2φφNNuTT=0}{=∂∂φV将(8)式代入(7)式,有:(9)(9)式的推导中,用到如下公式:利用昀小势能法,有:(10)(9)式代入(10)式,可得:(11)dxNPdxxNxNEILTLT∫∫=⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂002222][}{φdxNPdxxNxNEILTLT∫∫=⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂002222][}{φ3.1.5刚度矩阵整理后,有:(12)在(12)式中,令:,4402222][×∈⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂=∫RdxxNxNEILTK}]{[}{XKF=可得有限元的力平衡方程式的标准形式:(13)140][}{×∈⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧==∫RMFMFdxNPjjyiiyLTF{}14][×∈⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡==RuuXjjiiθθφ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−=⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂=∫LEILEILEILEILEILEILEILEILEILEILEILEILEILEILEILEIdxxNxNEILT46266126122646612612][22232322232302222K刚度矩阵为:][K(14)]232231[23233222323322LxLxLxLxLxLxxLxLxN+−−+−+−=其中,形函数为:⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡jjiijjyiiyuuLEILEILEILEILEILEILEILEILEILEILEILEILEILEILEILEIMFMFθθ46266126122646612612222323222323}]{[}{XKF=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=jjyiiyMFMF}{F⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=jjiiuuXθθ][可得:将刚度矩阵代入(13)式:(15)⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡−jjyiiyjjiiMFMFLEILEILEILEILEILEILEILEILEILEILEILEILEILEILEILEIuu122232322232346266126122646612612θθ由(15)式可得节点的位移和转角:(16)3.2Matlab程序说明4402222][×∈⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂=∫RdxxNxNEILTKclearx=sym(‘x’);L=sym(‘L’);N=[1–3*(x^2)/(L^2)+2*(x^3)/(L^3),x–2*(x^2)/L+(x^3)/(L^2),….3*(x^2)/(L^2)–2*(x^3)/(L^3),-(x^2)/L+(x^3)/(L^2)];Ni=diff(N,2);Nt=transpose(Ni);kk=Nt*Ni;K=E*I*int(kk,0,‘L’);]232231[23233222323322LxLxLxLxLxLxxLxLxN+−−+−+−=几个Matlab命令说明sym命令x=sym(‘y’):建立符号变量x,变量的值为单引号内的字符或字符串y例如:建立多个符号变量:symsxyz注意:x、y、z之间不能加逗号,只能为空格键(1)x=sym(‘y’)得出x=y(2)x=sym(‘y+z’)得出x=y+z(3)y=sym(‘a*x^2+b*x+c’)得出y=a*x^2+b*x+csubs命令x=subs(y,’old’,’new’):将符号函数y中的变量“old”用新变量“new”替代(1)x=sym(‘y’);z=subs(x,‘y’,2)得出z=2(2)x=sym(‘y+z’);zz=subs(x,‘[y,z]’,[2,3])得出zz=5如果继续进行替换,如计算当x=4时的z的值,则可subs(z,‘x’,4)将得出结果:27(3)y=sym(‘a*x^2+b*x+c’);z=subs(y,‘[a,b,c]’,[1,2,3])得出z=x^2+2*x+3注意subs命令中引号的使用,例如上面(1)中,不能写成:z=subs(x,y,2)或z=subs(‘x’,y,2)或z=subs(‘x’,‘y’,2)•diff(N,‘x’,n)将函数N对x进行n次微分,命令也可简写为diff(N,n),Matlab默认值即是对x微分nnxN∂∂Matlab程序:x=sym(‘x’);fx=x^2+2*x+1;first_order=diff(fx,’x’,1)或者first_order=diff(fx,1)second_order=diff(fx,’x’,2)或者second_order=diff(fx,2)运行结果:first_order=2*x+2second_order=212)(2++=xxxf,xxf∂∂)(22)(xxf∂∂求:52),(223+++=xyyxyxyxf,yyxf∂∂),(22),(xyxf∂∂求:Matlab程序:x=sym(‘x’);y=sym(‘y’);或者symsxyfxy=(x^3)*(y^2)+2*(x^2)*y+x*y+5;first_order_y=diff(fxy,’y’,1)说明:此时必须指明是对y微分second_order_x=diff(fxy,’x’,2)或者second_order_x=diff(fxy,2)运行结果:first_order_y=2*(x^3)*y+2*(x^2)+xsecond_order_x=6*x*(y^2)+4*y•int(N,‘x’,L1,L2);int(N,‘x’,L1,L2)代表对函数N进行积分,积分区间为L1到L2。命令也可简写为diff(N,L1,L2),Matlab默认值既是对x积分∫21LLNdx3131103102==∫xdxxMatlab程序:x=sym(‘x’);fx=x^2;result=int(fx,‘x’,0,1)或者result=int(fx,0,1)运行结果:result=1/3注意:此时result得出的是符号值,若将变量result代入数值计算,昀好先采用numeric命令得出result的数值,格式为:numeric(变量名)例如:result=numeric(result)将得出result=0.3333yyyxyxdxxyyx2131)2131()(3102331032+=+=+∫Matlab程序:x=sym(‘x’);y=sym(‘y’);或者symsxyfxy=(x^2)*(y^3)+x*y;result=int(fxy,‘x’,0,1)或者result=int(fxy,0,1)运行结果:result=1/3*y^3+1/2*y要求变量result在y=6时的数值,可以:result=subs(result,‘y’,6)运行结果:result=753.3悬臂梁的静力分析——集中载荷如图所示,一个方形截面的梁,截面每边长为5cm,在左端约束固定后,在右端末Y轴方向施加一集中力FY。试用有限元法求解末端挠度。问题条件L=10m;b=h=5e-2m;E=3e10N/m2;FY=100Nbh截面放大图:10m100N3.3.1解题方法1.由于只有一个集中力,故此将整个梁分为两个元素,如下图所示:2.梁的惯性矩mbhI73223102.512)105)(105(12−−−×=××==232F3=100NM3=0F2=0M2=0121F1=0M1=0F2=0M2=0《单元1》《单元2》10m100N3.求解单元1的刚度矩阵《单元1》i=1,j=2,L=5m,E=3e11N/m2,I=5.2e-7m4将已知条件代入,可得到单元1的刚度矩阵位:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−−−−=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−=124800374406240037440374401497637440149766240037440124800374403744014976374401497646266126122646612612][2223232223231LEILEILEILEILEILEILEILEILEILEILEILEILEILEILEILEIK121F1=0M1=0F2=0M2=04.单元1节点的受力将单元1节点上所承受的集中力与弯矩写成矩阵形式,有:⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=0000}{2211MFMFF121F1=0M1=0F2=0M2=0单元1的受力图5.求解单元1节点上的挠度与转角将单元1的刚度矩阵和受力向量代入有限元的力平衡方程式,可以得到节点与节点的挠度与转角:}]{[}{XKF=12⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−−−−=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧2211221111248003744062400374403744014976374401497662400374401248003744037440149763744014976][0000θθθθuuuuK6.求解单元2的刚度矩阵《单元2》i=2,j=3,L=5m,E=3e11N/m2,I=5.2e-7m4将已知条件代入,可得到单元1的刚度矩阵位:][124800374406240037440374401497637440149766240037440124800374403744014976374401497646266126122646612612][122232