Matlab应用之机械零件设计.ppt

工程优化设计与Matlab实现,之工程应用篇,Part:,（一）机械零件设计：向量、受力和刚体平衡梁的挠度计算四连杆机构运动设计与分析凸轮轮廓（二）动力学与振动：轨迹单自由度系统两自由度系统,（三）优化设计（重点）：线性规划无约束优化单目标约束优化约束优化问题（四）工程统计,第5章机械零件设计,教学目标介绍运用Matlab进行各种类型机械零件设计的分析方法,学习要求能够运用Matlab的内置函数或是自编程序，进行简单的受力分析和机构运动分析。,目录,5.1向量、受力和刚体平衡5.2梁的挠度计算5.3连杆机构运动设计与分析5.4凸轮轮廓习题,5.1向量、受力和刚体平衡,5.1.1向量5.1.2力的分解与合成5.1.3刚体的平衡,5.1向量、受力和刚体平衡,5.1.1向量,对于向量a,向量a的点乘定义为,向量a的值则为,将向量a表示为,则点乘用函数表示为,其值则为,5.1向量、受力和刚体平衡,向量a方向余弦,上式为单位向量ua在a方向上的分量，可以写成,因而有,其中,单位向量ua,5.1向量、受力和刚体平衡,例5.1一个力系如下图所示，确定合力F的值及其方向余弦角的大小。,程序如下：F1=04060;F2=60-11070;F=norm(F1+F2)%求合力的值alpha=acos(F1+F2)/F)*180/pi,5.1.2力的合成及分解,将弧度转化为角度,例5.2一个力系如下图所示，分解作用在D点的合力F，并确定该力的方向余弦角。,程序如下：f=24-16-48;uf=f/norm(f);F=30*ufalpha=acos(uf)*180/pi,运行结果：F=12.8571-8.5714-25.7143alpha=64.6231106.6015148.9973,5.1向量、受力和刚体平衡,F=30,例5.3一个力系如下图所示，求FL和FR的合力F及其分量。,程序如下：rl=2.50-3;rr=2.50.5-3;F=35*rl/norm(rl)+25*rr/norm(rr)resultant=norm(F),运行结果FF=38.28153.1750-45.9378FH=59.8818,5.1向量、受力和刚体平衡,例5.4一个力系如下图所示，确定保持该力系平衡所需的力的分量、大小及其方向余弦角。,程序如下：F1=02750;F2=00-150;F4=-10000;r=-2-67;F3=400*r/norm(r);F=-(F1+F2+F3+F4)FH=norm(F)alpha=acos(F/FH)*180/pi,F=184.7998-20.6005-146.7994FH=236.9080alpha=38.735094.9885128.2904,5.1向量、受力和刚体平衡,5.1.3刚体的平衡,5.1向量、受力和刚体平衡,如下图所示杆系，设已知：G1=200;G2=100;L1=2;L2=sqrt(2);theta1=30*pi/180;theta2=45*pi/180;求其支撑反力Na,Nb,Nc。,两杆系统的受力图（左）分离体受力图（右）,例5.5双杆的平衡,X=0Nax+Ncx=0Y=0Nay+Ncy-G1=0;M=0Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-G1*L1/2*cos(theta1)=0;X=0Nbx-Ncx=0;Y=0Nby-Ncy-G2=0M=0Ncy*L2*cos(theta2)+Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0;这是一组包含六个未知数Nax,Nay,Nbx,Nby,Ncx,Ncy的六个线性代数方程,通常是要寻找简化的方法,但利用MATLAB工具，就可以列出矩阵方程AX=B,(其中X=Nax,Nay,Nbx,Nby,Ncx,NcyT，可用矩阵除法直接来解。,对杆件1和2：,5.1向量、受力和刚体平衡,程序如下：%给原始参数赋值G1=200;G2=100;L1=2;L2=sqrt(2);%将度化为弧度theta1=30*pi/180;theta2=45*pi/180;%则按此次序，系数矩阵A,B可写成下式A=1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-L1*sin(theta1),L1*cos(theta1);0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0,L2*sin(theta2),L2*cos(theta2)B=0;G1;G1*L1/2*cos(theta1);0;G2;-G2*L2/2*cos(theta2)X=AB;%用左除求解线性方程组,5.1向量、受力和刚体平衡,这样求解的方法不仅适用于全部静力学题目，而且可用于材料力学和结构力学中的超静定问题。因为那里只多了几个形变变量和变形协调方程，通常也是线性的，所以只不过是把矩阵方程扩大了几阶，解法没有什么差别。,5.1向量、受力和刚体平衡,例5.6长为L(m)的悬臂梁左端固定，在离固定端L1(m)处施加力P(N)，求它的转角和挠度。设梁E=200e9(N/m2)和I=2e-5(m4)为已知。,建模：材料力学中从弯矩求转角要经过一次不定积分，而从转角求挠度又要经过一次不定积分，在MATLAB中这却是非常简单的问题。可用cumsum函数作近似的不定积分，还可用更精确的函数cumtrapz来做不定积分。本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯矩为,5.2、梁的挠度计算,程序（compute_naodu.m）如下L=2;P=2000;L1=1.5;%给出常数E=200e9;I=2e-5;x=linspace(0,L,101);dx=L/100;%将L分100段n1=L1/dx+1;%确定x=L1处对应的下标M1=-P*(L1-x(1:n1);%第一段弯矩赋值M2=zeros(1,101-n1);%第二段弯矩赋值(全为零)M=M1,M2;%全梁的弯矩A=cumsum(M)*dx/(E*I);%对弯矩积分求转角Y=cumsum(A)*dx;%对转角积分求挠度subplot(3,1,1),plot(x,M),grid,title(弯矩图)%绘弯矩图subplot(3,1,2),plot(x,A),grid,title(转角图)%绘转角图subplot(3,1,3),plot(x,Y),grid,title(挠度图)%绘挠度图,5.2、梁的挠度计算,cumsum函数：元素的累计和例：A=1,2,3,4cumsum(A)ans=13610在这里用来近似求数值积分,所得的结果见右，注意几根曲线之间的积分关系。本题之所以简单，除了用cumsum来近似不定积分之外，还因为在x=0处，虽然弯矩最大而转角和挠度都为零，因此两次积分的积分常数恰好都为零。如果它不为零，程序中就得有确定积分常数的语句，在下例中将能看到。,5.2、梁的挠度计算,例5.7简支梁受左半均匀分布载荷q及右边L/4处集中力偶Mo作用，求其弯矩、转角和挠度。设L=2m,q=1000N/m,M0=900Nm,E=200e9,I=2e-6。,建模:此题解法基本上与上例相同，主要差别是要处理积分常数问题。支撑反力可由平衡方程求得，设Q=qL/2，则,5.2、梁的挠度计算,各段弯矩方程为：,5.2、梁的挠度计算,对M/EI积分，得转角A，再作一次积分，得挠度Y。每次积分都要出现待定一个常数,其中A0(x)=cumtrapz(M)*dx/EI，Y0(x)=cumtrapz(A0)*dx。,累计梯形积分函数,5.2梁的挠度计算,两个待定积分常数Ca和Cy可由边界条件Y(0)=0及Y(L)0确定：Y(0)=Y0(0)+Cy=0Y(L)=Y0(L)+Ca*L+Cy=0于是可得：,即,5.2梁的挠度计算,%输入已知参数L，q，Mo，E，IL=2;q=1000;Mo=900;E=200e9;I=2e-6;Na=(3*q*L2/8-Mo)/L;Nb=(q*L2/8+Mo)/L;x=linspace(0,L,101);dx=L/100;%用数组分三段列出M的表达式M1=Na*x(1:51)-q*x(1:51).2/2;M2=Nb*(L-x(52:76)-Mo;M3=Nb*(L-x(77:101);M=M1,M2,M3;%列写完整的M数组A0=cumtrapz(M)*dx/(E*I);%由M积分求转角Y0=cumtrapz(A0)*dx;%由转角积分求挠度C=0,1;L,1-Y0(1);-Y0(101);%左除求积分常数Ca,Cy,5.2梁的挠度计算,Ca=C(1),Cy=C(2),A=A0+Ca;Y=Y0+Ca*x+Cy;%转角与挠度全值subplot(3,1,1),plot(x,M),grid，subplot(3,1,2),plot(x,A),gridsubplot(3,1,3)plot(x,Y),grid,5.2梁的挠度计算,5.2梁的挠度计算,5.3连杆机构的运动设计及分析,连杆机构的运动分析就是在已知机构的运动尺寸和已知原动件的运动规律的前提下，确定机构中其他构件的位置、速度和加速度等运动参数。机构运动分析的方法有矢量方程图解法和解析法,5.3连杆机构的运动设计及分析,对于下图所示的铰链四杆机构，已知机构主动件AB与从动件CD的位置对应关系，根据各个构件长度在直角坐标系中的投影，建立机构的位置方程式,5.3.1给定连架杆对应位置的设计问题,整理上式可得因为连架杆的运动取决于各个构件的相对长度，设机构的相对杆件长度系数为,5.3连杆机构的运动设计及分析,5.3连杆机构的运动设计及分析,将以上系数代入机构的位置方程式，得到铰链四杆机构的位置参数方程当已知连架杆的对应位置关系，上式中有5个参数，说明在四杆机构常规设计中，能够满足两连架杆的对应位置关系最多为5组。,5.3连杆机构的运动设计及分析,例5.8：铰链四杆机构设计已知铰链四杆机构两连架杆AB、CD的初始位置角，它们三组的对应位置以及机架的长度L4=50mm。,分析：将已知参数代入,求解线性方程组，可得到R1、R2、R3，最后得出其它三个构件的长度。,5.3连杆机构的运动设计及分析,程序如下：%实现连架杆角位移(3组)的连杆机构运动设计%已知条件f0=0;p0=0;%连架杆初始位置角f=4590135*pi/180;%曲柄输入角p=5282112*pi/180;%摇杆输出角%杆件相对长度参数R1、R2和R3的系数矩阵a1=1-cos(f(1)+f0)cos(p(1)+p0);a2=1-cos(f(2)+f0)cos(p(2)+p0);a3=1-cos(f(3)+f0)cos(p(3)+p0);a=a1;a2;a3%线性方程组右边的常数矩阵b1=cos(f(1)-p(1)+(f0+p0);b2=cos(f(2)-p(2)+(f0+p0);b3=cos(f(3)-p(3)+(f0+p0);b=b1b2b3%线性方程组aR=b直接解法(采用求逆函数inv)R=inv(a)*b%或采用矩阵除法R=ab,5.3连杆机构的运动设计及分析,d=50;%机架长度x(1)=d/R(3);x(3)=d/R(2);x(2)=sqrt(x(1)2+x(3)2+d2-2*x(1)*x(3)*R(1);%检验解的精度en=norm(a*R-b);disp*计算结果*fprintf(1,曲柄长度a=%3.4fmmn,x(1);fprintf(1,连杆长度b=%3.4fmmn,x(2);fprintf(1,摇杆长度c=%3.4fmmn,x(3);fprintf(1,机架长度d=%3.4fmmn,d);dispfprintf(1,数值解的精度en=%3.4en,en);,5.3连杆机构的运动设计及分析,对于下图所示的铰链四杆机构,5.3.2四杆机构的位置、速度和加速度,对于一个特定的四杆机构，已知其各构件的长度和原动件2的运动规律，即为已知，而=0,故可求得未知方位角,。,5.3连杆机构的运动设计及分析,角位移方程的分量形式为：,上式对时间求解一阶导数，就能够得到角速度方程，如下所示：,其矩阵形式为：,5.3连杆机构的运动设计及分析,联立可求得：,对时间求解二阶导数，得到角加速度方程矩阵形式,求导中应用了下列公式：,5.3连杆机构的运动设计及分析,由上式可求得加速度：,求导中应用了下列公式：,5.3连杆机构的运动设计及分析,例5.9：一个四杆机构，已知其各构件的长度L1=304.8:mm,L2=101.6mm,L3=254.0mmL4=177.8mm,在30至360度之间。,编写程序如下：,(1)建立函数fourbarposition.m,functiont=fourbarposition(th,th2,L2,L3,L4,L1)t=L2*cos(th2)+L3*cos(th(1)-L4*cos(th(2)-L1;L2*sin(th2)+L3*sin(th(1)-L4*sin(th(2);,%给定已知量，各杆长L1,L2,L3,L4L1=304.8;L2=101.6;L3=254.0;L4=177.8;th2=1/6:1/6:2*pi;%曲柄输入角度从30至360度，步长为pi/6th34=zeros(length(th2),2);%建立一个N行2列的零矩阵，第一列存放3，第二列存放4options=optimset(display,off);form=1:length(th2)%建立for循环，求解3，4th34(m,:)=fsolve(fourbarposition,55,options,th2(m),L2,L3,L4,L1);%调用fsolve函数求解关于3，4end%解非线性方程组，结果保存在th34中,5.3连杆机构的运动设计及分析,y=L2*sin(th2)+L3*sin(th34(:,1);%连杆3的D端点Y坐标值x=L2*cos(th2)+L3*cos(th34(:,1);%连杆3的D端点X坐标值xx=L2*cos(th2);%连杆3的C端点X坐标值yy=L2*sin(th2);%连杆3的C端点Y坐标值figure(1)plot(x;xx,y;yy,k,0L1,00,k-,x,y,ko,xx,yy,ks)%绘制连杆3的几个位置点holdonth=linspace(0,2*pi,100);plot(L2*cos(th),L2*sin(th),k-)%绘制连杆2与连杆3连接点的轨迹title(连杆3的几个位置点)xlabel(水平方向)ylabel(垂直方向)axisequal%xy坐标均衡,5.3连杆机构的运动设计及分析,5.3连杆机构的运动设计及分析,%给定已知量，各杆长L1,L2,L3,L4L1=304.8;L2=101.6;L3=254.0;L4=177.8;th2=(0:2/72:2)*pi;%重新细分曲柄输入角度2，步长为5度th34=zeros(length(th2),2);options=optimset(display,off);form=1:length(th2)th34(m,:)=fsolve(fourbarposition,11,options,th2(m),L2,L3,L4,L1);end,5.3连杆机构的运动设计及分析,figure(2)plot(th2*180/pi,th34(:,1)*180/pi,th2*180/pi,th34(:,2)*180/pi)%绘制连杆3和摇杆4的角位移关于曲柄2的角位移图axis(03600170)%确定XY边界值grid%图形加网格xlabel(主动件转角theta_2(度)ylabel(从动件角位移(度)title(角位移线图)text(120,120,摇杆4角位移)text(150,40,连杆3角位移),5.3连杆机构的运动设计及分析,w2=250;%设定曲柄角速度fori=1:length(th2)A=-L3*sin(th34(i,1)L4*sin(th34(i,2);L3*cos(th34(i,1)-L4*cos(th34(i,2);B=w2*L2*sin(th2(i);-w2*L2*cos(th2(i);w=AB;w3(i)=w(1);w4(i)=w(2);end,5.3连杆机构的运动设计及分析,figure(3)plot(th2*180/pi,w3,th2*180/pi,w4);%绘制角速度图axis(0360-175200)text(50,160,摇杆4角速度(omega_4)text(220,130,连杆3角速度(omega_3)gridxlabel(主动件转角theta_2(度)ylabel(从动件角速度(radcdots-1)title(角速度线图),5.3连杆机构的运动设计及分析,fori=1:length(th2)C=-L3*sin(th34(i,1)L4*sin(th34(i,2);L3*cos(th34(i,1)-L4*cos(th34(i,2);D=w22*L2*cos(th2(i)+w3(i)2*L3*cos(th34(i,1)-w4(i)2*L4*cos(th34(i,2);w22*L2*sin(th2(i)+w3(i)2*L3*sin(th34(i,1)-w4(i)2*L4*sin(th34(i,2);a=inv(C)*D;a3(i)=a(1);a4(i)=a(2);end,5.3连杆机构的运动设计及分析,figure(4)plot(th2*180/pi,a3,th2*180/pi,a4);%绘制角加速度图axis(0360-7000065000)text(50,50000,摇杆4角加速度(alpha_4)text(220,12000,连杆3角加速度(alpha_3)gridxlabel(主动件转角theta_2(度)ylabel(从动件角加速度(radcdots-2)title(角加速度线图)disp曲柄转角连杆转角-摇杆转角-连杆角速度-摇杆角速度-连杆加速度-摇杆加速度ydcs=th2*180/pi,th34(:,1)*180/pi,th34(:,2)*180/pi,w3,w4,a3,a4;disp(ydcs),5.3连杆机构的运动设计及分析,运行结果：,5.3连杆机构的运动设计及分析,5.3连杆机构的运动设计及分析,5.3连杆机构的运动设计及分析,5.4凸轮轮廓,凸轮是把一种运动转化为另一种运动的装置。凸轮的廓线和从动件一起实现运动形式的转换。凸轮通常是为定轴转动，凸轮旋转运动可被转化成摆动、直线运动或是两者的结合。凸轮机构设计的内容之一是凸轮廓线的设计。,5.4凸轮轮廓,定义一个凸轮基圆rb作为最小的圆周半径。从动件的运动方程如下：,设凸轮的推程运动角和回程运动角均为，从动件的运动规律均为正弦加速度运动规律，则有：,上式是从动件的位移，h是从动件的最大位移，并且0。,5.4凸轮轮廓,如果假设凸轮的旋转速度是个常量，则速度、加速度a和瞬时加速度j（加速度对时间求导）分别如下所示。,定义无量纲位移S=s/h、无量纲速度V=/h、无量纲加速度A=a/h3和无量纲瞬时加速度J=j/h3。,位移S,5.4凸轮轮廓,速度V,加速度A,5.4凸轮轮廓,瞬时加速度J,5.4凸轮轮廓,例5.10：摆线凸轮的位移、速度、加速度、和瞬时加速度计算当60时，凸轮的无量纲位移、速度、加速度和瞬时加速度的值，并绘制图形。,程序如下(CamContour.m)：beta=60*pi/180;phi=linspace(0,beta,40);%0phibetaphi2=beta+phi;%betaphi22*betaph=phiphi2*180/pi;%将弧度转化为角度arg=2*pi*phi/beta;arg2=2*pi*(phi2-beta)/beta;,5.4凸轮轮廓,s=phi/beta-sin(arg)/2/pi1-(arg2-sin(arg2)/2/pi;v=(1-cos(arg)/beta-(1-cos(arg2)/beta;a=2*pi/beta2*sin(arg)-2*pi/beta2*sin(arg2);j=4*pi2/beta3*cos(arg)-4*pi2/beta3*cos(arg2);,5.4凸轮轮廓,subplot(2,2,1)plot(ph,s,k)xlabel(凸轮推程运动角(度)ylabel(位移(S)g=axis;g(2)=120;axis(g)%返回四维向量，x轴和y轴的坐标范围subplot(2,2,2)plot(ph,v,k,0120,00,k-)xlabel(凸轮推程运动角(度)ylabel(速度(V)g=axis;g(2)=120;axis(g),5.4凸轮轮廓,subplot(2,2,3)plot(ph,a,k,0120,00,k-)xlabel(凸轮推程运动角(度)ylabel(加速度(A)g=axis;%返回四维向量，x轴和y轴的坐标范围g(2)=120;axis(g)subplot(2,2,4)plot(ph,j,k,0120,00,k-)xlabel(凸轮推程运动角(度)ylabel(瞬时加速度(J)g=axis;g(2)=120;axis(g),5.4凸轮轮廓,运行结果,习题,已知F和W两个力，F1，F2，F3的方向已知。力的平衡方程为：F1+F2+F3=F+WF1=f1*ur1;F2=f2*ur2;F3=f3*ur3;f1f2f3ur1;ur2;ur3=Rf1f2f3*A=R,