M atlab求解理论力学问题系列(二)典型机构的运动分析

高云峰

(清华大学航天航空学院,北京100084)

在传统的运动学教学中，为了避免复杂的求导而陷入纯数学演算，通常强调基点法和点的复合运动等具有物理含义的方法，求解系统在特定时刻或位置的运动信息[1]。这样处理的优点是强调物理概念，缺点是对系统的整体运动不容易了解，某些点的运动轨迹有时难以想象。这样训练出来的学生，对于机械运动缺乏了解，也难以胜任未来的设计工作。

利用运动学计算机辅助分析方法和相关软件，可以很容易演示系统的整个运动过程、求出任意点或刚体在任意时刻的速度、加速度、角速度、角加速度等运动量。Matlab是实现计算机辅助分析的工具之一，可以帮助学生加深理解机械运动的特点，为将来学生自己设计装置提供了方便。

本文着重介绍2个典型机构的运动轨迹、速度和加速度分析，以及利用Matlab中的符号推导功能，从位移函数求导得到速度和加速度的表达式，并进行画图和动画显示。

在传统点的复合运动中，要特别注意动点、动系的选择，如果选择不当，就分析不下去，因此对很多同学而言有一定的难度。而采用Matlab，可以采用统一的方法，只要列出一般位置的表达式就可以利用符号推导功能，得到速度、加速度等运动量的表达式。

1 M atlab中非线性方程的求解及动画演示

案例1：如图1，已知四连杆机构ABCD，AB杆长为a1，BC杆长为a2，CD杆长为a3，AD距离为a4。若AB杆以匀角速度ω0转动，初始θ0=0。求BC和CD杆的角度、角速度变化规律。

图1 四连杆机构

这个问题中尺寸不是特别给定的值，一般学生分析起来会有困难。而用Matlab可以进行规范化处理：建立坐标Axy，θ=ω0t已知，广义坐标为φ1和φ2。系统的约束方程为[2]

方程(1)是关于转角φ1和φ2的非线性方程组，通常没有解析解，下面给出一般的处理方法。

设有非线性方程组f i(x1,x2,···,x n),i=1,2,···,n，没有解析解，可用多种数值计算的方法求解，这里以牛顿−拉普森法为例，其主要步骤如下：

步骤(1)：给定计算中的允许误差ε，给出一组粗略的估计值

步骤(2)：计算f i(x∗)，判断|f i(x∗)|<ε是否成立，若成立，当前估计值即为一组数值解，求解结束。若不成立，到步骤(3)。

步骤(3)：计算非线性方程组的雅可比矩阵J(x)，代入当前估计值得到J(x∗)。多元函数f i(x1,x2,···,x n)的雅可比矩阵定义为

步骤(4)：类似一元函数的泰勒展开式，f(x)=f(x0)+f′(x0)(x−x0)+o(x−x0)，多元函数为

略去高阶小量，从而得到关于修正值dx=[dx1,dx2,···,dx n]T的一组线性方程组

步骤(5)：解线性方程组，得到修正值dx。

步骤(6)：得到新的估计值x∗=x∗+dx，返回到步骤(2)。

具体回到案例1，把方程(1)改写为

已知a1,a2,a3,a4和θ，估计值是计算雅可比矩阵，并代入当前估计值

得到一组关于修正值[dφ1,dφ2]T的线性方程组

类似代码X=inv(A)*B，可解出[dφ1,dφ2]T，总之按上述步骤1～6，可以解出任意时刻的角度。

在某些情况下，雅可比矩阵的行列式可能趋近于零或等于零，意味着方程(4)多解或无解。多解可能的原因是：机械系统由于设计导致在某个特定位置本身就是多解的，如图2中的曲柄−滑块机构，OA运动带动B运动通常不会有问题；但反过来B运动带动OA运动(如蒸汽机)，当OAB共线时，OA下一时刻可以向上也可以向下(考虑动力学时OA杆由于惯性继续运动，解是确定的，但仅考虑运动学时OA杆位置是多解的)。无解的可能是：机械系统由于尺寸关系导致运动发生冲突，例如图2中l

图2 无解或多解的情况

编程计算得到角度的变化关系后，可以算出任意时刻各铰的位置，以及BC杆上不同点的运动轨迹(图3)：很明显B点轨迹是圆，C点轨迹是圆的一部分(AB杆大范围运动时，CD杆只在小范围运动)，而在BC杆上不同的点轨迹就很复杂了。

图3 BC杆上不同点的运动轨迹

知道任意时刻系统的各铰位置，在屏幕上画出各铰点的位置并连接起来，就得到了四连杆机构在某一时刻的图象，延迟一定的时间后再画出下一时刻的图象，就形成了动画。本问题中动画的源代码见图4，其中plot函数表示画线段；h1是句柄，定义动画中每一帧图的特性，如颜色、线型、宽度等；set函数是读取数据；drawnow函数是画出当前帧的图案；pause函数表示暂停，让每一帧有一定的视觉残留以便更好地形成动画(类似放电影，1秒钟放24格图像)[3]。

图4 动画部分的源代码截图

2 函数求导获得角速度、角加速度

求出角度后，对约束方程(1)求导出现角速度，有

由于θ,φ1,φ2已在前面求出，因此得到关于 ˙φ1,˙φ2的一组线性方程组。类似X=inv(A)*B可解出角速度，从而可以获得角速度随时间或随θ变化的关系(图5)。

对式(5)进一步求导得到角加速度的方程

此时θ,φ1,φ2,˙φ1,˙φ2都是已知值，因此得到关于¨φ1,¨φ2的一组线性方程组。类似X=inv(A)*B可解出角加速度，从而可以获得角加速度随时间或随θ变化的关系(图6)。

这只是一类典型机构的例子，采用本文介绍的方法，任意机构的运动轨迹、刚体的角速度和角加速度、某点的速度和加速度都可以类似求出。而这些结果(图3、图5、图6)都是传统方法难以获得的。

图5 角速度与θ的关系

图6 角加速度与θ的关系

3 M atlab中的符号求导

如果运动学问题中没有具体数值，全是参数，用Matlab也可以处理。

案例2：半径为r的圆轮在水平桌面上作直线纯滚动，轮心速度v O的大小为常数。摇杆AB与桌面铰接，并靠在圆轮上(图7)。当摇杆与桌面夹角φ=60°时，试求摇杆的角速度和角加速度。

图7 一类接触问题

传统处理方法：取AB杆为动系，O为动点(图8)。通过速度和加速度分析(具体过程略)，可以得到

图8 速度分析

Mtlab处理本问题时，只需要把一般位置角度的表达式写出来，然后利用符号求导的方法，就可以得到角速度和角加速度。角度以逆时针方向为正，设初始时刻AB杆垂直，根据图9有

图9 运动学关系

根据式(8)有

对式(9)求一阶导数获得角速度，求二阶导数获得角加速度。程序的源代码见图10。

因此可以画出系统运动全程的角度、角速度、角加速度变化曲线，取归一化参数r=1，v O=1，得到的曲线关系见图12和图13。

图12 系统运动参量与时间的关系

图13 系统运动参量与角度的关系

源代码中具体涉及几个函数：(1)符号定义是“syms”。(2)函数求导的格式是diff(function,’variable’)，表示函数对某个变量求导。运行后得到的结果为

可以看到，传统方法的分析要很大的篇幅，而Matlab只需要几行就得到了结果(顺便说一下，前面式(5)和式(6)也可以利用符号推导从式(1)得到)。

如果要求特定位置的值，只需要把具体值代入，格式为subs(function,old,new)，表示用新的变量替换老的变量，这部分程序的源代码见图11。

图11 求特定位置的运动量

运行后屏幕显示结果为

借助Matlab还可以对题目进行深入研究。例如，用传统方法分析时，如果动点选为圆盘上的接触点，速度还可能会正确，但是加速度没有办法分析了，这是为什么呢？可以看看接触点相对AB杆的轨迹，这借助计算机很容易实现。

以φ=60°为系统初始位置，设圆盘上P点与AB杆接触，且注意负号表示顺时针方向，因此计算机推导出的结果与传统方法相同，见式(7)。

利用参数替换函数，还可以把式(1)和式(11)中的时间用角度替换(具体略)。当圆心运动Δx=v O t时，圆盘因为纯滚动的转角为θ=Δx/r(图14)，且P点在动系中有

图14 接触点的轨迹分析

在动系中画出的轨迹如图15所示，很像是旋轮线。如果P点的轨迹是旋轮线，则曲线的尖点应该接触AB杆；而实际上P点接触杆时，相对速度不为零，但轨迹的切线方向与杆平行，不过该曲线的曲率在详细分析前还不清楚。

这就回答了前面的问题：如果选P点为动点，在分析速度时认为相对速度沿杆方向(相信这样做的学生自己也不能说得很清楚，只是凭感觉)，碰巧是正确的；但是在加速度分析时，由于相对运动的轨迹不是太清楚，相对切向加速度未知，相对法向加速度也未知，再加上AB杆的角加速度未知，未知数太多求解不了。

注意在地面上看P点的轨迹是旋轮线，通过反转和平移，旋轮线与图15中的轨迹完全相同，这个结论是否出乎意料？

4 小结

在运动学中引入Matlab，可以方便了解系统整个运动过程，获得系统中任意点的运动轨迹，以及速度和加速度随时间或某角度的变化规律，直观明了，让学生更容易理解系统的运动。

本篇着重介绍了非线性方程组的求解方法，以及画图和动画的格式；另外介绍了参数方程的求导和替换。画图涉及到的Matlab函数是plot(画直线)，而句柄函数、set函数和drawon函数的组合就可以实现动画；另一个重要函数是diff(函数求导)和subs(变量替换)。掌握这几个命令后，就可以对一般的运动学问题进行全程的计算和动画演示，获得丰富的运动学信息。更进一步，借助Matlab可以更加深入研究运动学的问题，解决传统分析方法无法处理的很多问题。