实验二：曲柄滑块机构等效转动惯量及其倒数规律 题目：题目要求求解出下图不对心曲柄滑块机构等效力矩及其倒数的变化规律， 并画出曲线。 一、问题分析 要想找到转动惯量及其导数的关系，我们需要先根据等效公式列出等效转动 惯量的公式，然后对其求导，然后通过数值求解的方法对未知量分别进行求解， 然后将各自求得的未知变量的矩阵用循环的方式代入方程中求解得出一系列结 果，然后通过画图工具画出图形变化规律。 二、模型建立 首相我们通过平衡公式： 1*sin( 1) e l    l2*sin( )  可以得到：  =arcsin((l1*sin( 1) e) / l2)   又因为我们已知 x 点的横坐标公式： x  1*cos( 1) l l  ^ 2 (l1*sin( 1) e) ^ 2    所以我们可以得到 c 点的速度和加速度的表达式： dx 1 d  1 d  dt dx 1 d  dx dt 1* w vc *    ac   1* w dvc 1 d  d  d   * 1 dvc d dvc  1 dt d dt  d d    dt dt 2 d  2 dt d d  w1*d(  *  d  d  =w1* ) / d  此外，我们还已知 s2 点的横纵坐标的表达式： ) e  2 *sin( ls c  ys 2  xs 2  1*cos( 1) l   (l2*cos( )   ls 2*cos( ))  

然后我们可以求解 s2 点的 x,y 方向的速度、加速度的表达式： vs2x= 2 dxs dt 2 dxs d  * d  dt  1* w 2 dxs d      vs2y= 2 dys dt 2 dvs y dt  * 2 d dys  d dt  2 dvs y d  * d dt  1* w  1* w 2 dys d  2 dvs y d  2 as y  2 as x  2 dvs x dt 2 dvs x d  d dt  *  1* w 2 dvs x d  上式中：由于 w1 对方程的求解不产生影响，我们取 w1=1；先将各自的分 量通过符号求解求出公式，然后将 1 角划分 60 分，代入求解得到各分量的值， 然后代入以下转动惯量及其导数公式公式求解得到最终结果： 2 ^ 2  1^ 2 w 2 ^ 2 w 1^ 2 w ^ 2 vc 1^ 2 w Ja Js vs x 2* 2* 3* e J m m     vs y 2 ^ 2 d Je 1 d w  2 1^ 3 *[Js2*w 2*a 2 m3*vc*ac m2*(vs2x*as2x vs2 y*as2 y)]    三、求解代码 clear,clc l1=0.2;l2=0.5;ls2=0.2;e=0.05;w=1;;ja=3;js2=0.15;m2=5;m3=10; %给出已知变量 syms theta ; %定义符号变量 xc=l1*cos(theta)+sqrt(l2^2-(l1*sin(theta)-e)^2); %列出 C 点的横坐标公式 vc=w*diff(xc,theta); %通过符号求解的方法求解 C 点处的速度公式 ac=w*diff(vc,theta); %通过符号求解的方法求解 C 点处的加速度公式 bata=asin((l1*sin(theta)-e)/l2); %列出 Bata 的表达式 w2=diff(bata,theta); %求出 Bata 角的角速度表达式 a2=diff(w2,theta); ys2=ls2*sin(bata)+e; %列出 S2 点的 y 坐标方程 vys2=w*diff(ys2,theta); %求出 S2 点 Y 方向的速度表达式 xs2=l1*cos(theta)+(l2*cos(bata)-ls2*cos(bata)); %列出 S2 点的纵坐标表达式 vxs2=w*diff(xs2,theta); ays2=w*diff(vys2,theta); %求出 S2 点 Y 方向上的加速度表达式 axs2=w*diff(vxs2,theta); %求出 S2 点 X 方向上的加速度表达式 Je=[];DJe=[]; %求出 S2 点 X 方向的速度表达式 %求出 Bata 角的角加速度表达式 %建立空矩阵存储结果 theta=(0:6:360)*pi/180; %对 theta 角取值 vc=subs(vc,theta); ac=subs(ac,theta); w2=subs(w2,theta); vxs2=subs(vxs2,theta); vys2=subs(vys2,theta); a2=subs(a2,theta); %对 c 点的速度赋值 %对 c 点的加速度赋值 %对 bata 角的角速度赋值 %对 s2 点的 x 方向速度进行赋值 %对 s2 点的 y 方向速度进行赋值 %对 bata 角的角加速度进行赋值 

axs2=subs(axs2,theta); ays2=subs(ays2,theta); %对 X2 点的 x 方向的加速度进行赋值 %对 x2 点的 y 方向的加速度进行赋值 for i=1:60 %循环求解 je=vpa(ja+js2*w2(1,i)^2/w^2+m3*vc(1,i)^2/w^2+m2*(vxs2(1,i)^2+vys2(1,i)^2)/w^2); %将数值带入等效转动惯量的方程中求解 dje=vpa(2/w^3*(js2*w2(1,i)*a2(1,i)+m3*vc(1,i)*ac(1,i)+m2*(vxs2(1,i)*axs2(1,i)+vys2 (1,i)*ays2(1,i)))); %将数值代入 到等效转动惯量的求导公式中求解 Je=[Je je]; DJe=[DJe dje]; %记录等效转动惯量求导公式的结果 %记录等效转动惯量的结果 end [AX,H1,H2]=plotyy(Je,'-',DJe,'--') ; %画图 set(get(AX(1),'ylabel'),'string','Je/Kg.m^2'); %设置等效转动惯量的坐标属性 set(get(AX(2),'ylabel'),'string','(dJe/dtheta)/Kg.m^2');设置等效转动惯量的倒数纵坐 标属性 四、求解结果 五、实验小结 通过本次实验，我们进一步了解了曲柄滑块机构的变化规律，在列出求解微 分方程的同时，我们还采取数值计算的方法进行求解方程，加强了我们对实际方 程的求解能力，进一步拓展了我们的求解方法的知识面。