第 43 卷 第 3 期 力 学 与 实 践 2021 年 6 月 Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 高云峰1) (清华大学航天航空学院, 北京 100084) 摘要 运动学计算机辅助分析可以帮助人们方便了解机械系统的整个运动过程，获得系统中任意点的运 动轨迹、速度和加速度的变化规律，直观明了。Matlab 软件是实现计算机辅助分析的工具之一。本文着重介绍 了 2 个典型机构的运动分析，一个案例侧重从数值计算的角度介绍 Matlab 如何求解机构运动涉及的非线性方 程组和如何进行动画演示；另一个案例着重介绍如何利用 Matlab 中的符号推导功能，求解全参数化的运动学 问题。 关键词 非线性代数方程组，牛顿 − 拉普森方法，运动轨迹，动画显示，符号推导 KINEMATICS ANALYSIS OF TYPICAL MECHANISMS GAO Yunfeng1) (School of Aerospace Engineering, Tsinghua University, Beijing 100084, China) Abstract Computer-aided analysis of kinematics can help understand the whole movement process of a mechanical system, along with the motion trajectory, the velocity and the acceleration of the system at any time or position. Matlab can be used to realize the computer-aided analysis. This paper focuses on the kinematics analysis of two typical mechanisms. In one case, the focus is on how to solve the nonlinear algebraic equations by Matlab and how to perform the animation demonstration. In another case, the focus is on how to solve the fully parametric kinematics problem by using the symbolic derivation function in Matlab. Key words nonlinear algebraic equations, Newton–Raphson method, motion trajectory, animation demonstration, symbolic derivation 在传统的运动学教学中，为了避免复杂的求导 而陷入纯数学演算，通常强调基点法和点的复合运 动等具有物理含义的方法，求解系统在特定时刻或 位置的运动信息[1]。这样处理的优点是强调物理概 念，缺点是对系统的整体运动不容易了解，某些点的 运动轨迹有时难以想象。这样训练出来的学生，对于 机械运动缺乏了解，也难以胜任未来的设计工作。 利用运动学计算机辅助分析方法和相关软件，可 以很容易演示系统的整个运动过程、求出任意点或刚 体在任意时刻的速度、加速度、角速度、角加速度等 运动量。Matlab 是实现计算机辅助分析的工具之一， 可以帮助学生加深理解机械运动的特点，为将来学 生自己设计装置提供了方便。 本文着重介绍 2 个典型机构的运动轨迹、速度 和加速度分析，以及利用 Matlab 中的符号推导功 能，从位移函数求导得到速度和加速度的表达式，并 进行画图和动画显示。 在传统点的复合运动中，要特别注意动点、动系 的选择，如果选择不当，就分析不下去，因此对很多 同学而言有一定的难度。而采用 Matlab，可以采用 2020–06–01 收到第 1 稿，2020–07–23 收到修改稿。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (二) 典型机构的运动分析. 力学与实践, 2021, 43(3): 414-419 Gao Yunfeng. Kinematics analysis of typical mechanisms. Mechanics in Engineering, 2021, 43(3): 414-419
第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 415 统一的方法，只要列出一般位置的表达式就可以利 用符号推导功能，得到速度、加速度等运动量的表 达式。 1 Matlab 中非线性方程的求解及动画演示 案例 1：如图 1，已知四连杆机构 ABCD，AB 杆长为 a1，BC 杆长为 a2，CD 杆长为 a3，AD 距离 为 a4。若 AB 杆以匀角速度 ω0 转动，初始 θ0 = 0。 求 BC 和 CD 杆的角度、角速度变化规律。 a a a a ϕ ϕ y A B D C x θ 图 1 四连杆机构 这个问题中尺寸不是特别给定的值，一般学生 分析起来会有困难。而用 Matlab 可以进行规范化处 理：建立坐标 Axy，θ = ω0t 已知，广义坐标为 φ1 和 φ2。系统的约束方程为[2] a1 cos θ + a2 cos φ1 + a3 cos φ2 − a4 = 0 a1 sin θ + a2 sin φ1 + a3 sin φ2 = 0 } (1) 方程 (1) 是关于转角 φ1 和 φ2 的非线性方程组，通 常没有解析解，下面给出一般的处理方法。 设有非线性方程组 fi(x1, x2, · · · , xn), i = 1, 2, · · · , n，没有解析解，可用多种数值计算的方法求解， 这里以牛顿−拉普森法为例，其主要步骤如下： 步骤 (1)：给定计算中的允许误差 ε，给出一组 粗略的估计值 x ∗ = [x ∗ 1 , x∗ 2 , · · · , x∗ n] T。 步骤 (2)：计算 fi(x ∗ )，判断 |fi(x ∗ )| < ε 是否 成立，若成立，当前估计值即为一组数值解，求解结 束。若不成立，到步骤 (3)。 步骤 (3)：计算非线性方程组的雅可比矩阵 J(x)，代入当前估计值得到 J(x ∗ )。多元函数 fi(x1, x2, · · · , xn) 的雅可比矩阵定义为 J(x) = ∂f1/∂x1 ∂f1/∂x2 · · · ∂f1/∂xn ∂f2/∂x1 ∂f2/∂x2 · · · ∂f2/∂xn . . . . . . . . . . . . ∂fn/∂x1 ∂fn/∂x2 · · · ∂fn/∂xn 步骤 (4)：类似一元函数的泰勒展开式，f(x) = f(x0) + f ′ (x0)(x − x0) + o(x − x0)，多元函数为 fi(x) = fi(x ∗ ) + J(x ∗ )dx + o(dx) 略去高阶小量，从而得到关于修正值 dx = [dx1, dx2, · · · , dxn] T 的一组线性方程组 J(x ∗ )dx + fi(x ∗ ) = 0 步骤 (5)：解线性方程组，得到修正值 dx。 步骤 (6)：得到新的估计值 x ∗ = x ∗ + dx，返回 到步骤 (2)。 具体回到案例 1，把方程 (1) 改写为 f1(φ1, φ2) = a2 cos φ1 + a3 cos φ2 + a1 cos θ − a4 f2(φ1, φ2) = a2 sin φ1 + a3 sin φ2 + a1 sin θ } (2) 已知 a1, a2, a3, a4 和 θ，估计值是 (φ ∗ 1 , φ∗ 2 )。计算雅可 比矩阵，并代入当前估计值 J(φ ∗ 1 , φ∗ 2 ) = [ −a2 sin φ ∗ 1 −a3 sin φ ∗ 2 a2 cos φ ∗ 1 a3 cos φ ∗ 2 ] (3) 得到一组关于修正值 [dφ1, dφ2] T 的线性方程组 [ −a2 sin φ ∗ 1 −a3 sin φ ∗ 2 a2 cos φ ∗ 1 a3 cos φ ∗ 2 ] [ dφ1 dφ2 ] + [ f1(φ ∗ 1 , φ∗ 2 ) f2(φ ∗ 1 , φ∗ 2 ) ] = [ 0 0 ] (4) 类似代码 X=inv(A)*B，可解出 [dφ1, dφ2] T，总之按 上述步骤 1 ∼ 6，可以解出任意时刻的角度。 在某些情况下，雅可比矩阵的行列式可能趋近 于零或等于零，意味着方程 (4) 多解或无解。多解可 能的原因是：机械系统由于设计导致在某个特定位 置本身就是多解的，如图 2 中的曲柄−滑块机构，OA 运动带动 B 运动通常不会有问题；但反过来 B 运 动带动 OA 运动 (如蒸汽机)，当 OAB 共线时，OA 下一时刻可以向上也可以向下 (考虑动力学时 OA 杆 vB B B l l A A r r O O θ ϕ ω 图 2 无解或多解的情况
416 力 学 与 实 践 2021 年 第 43 卷 由于惯性继续运动，解是确定的，但仅考虑运动学 时 OA 杆位置是多解的)。无解的可能是：机械系统 由于尺寸关系导致运动发生冲突，例如图 2 中 l < r 时 OA 杆就不能连续转动，本问题中若四杆长度满 足 a2 + a3 < a1 + a4，当 θ = π 时方程 (1) 的确无解 (物理解释是 AB 杆运动范围有限制，不能达到这一 位置)。利用方程无解这一特点，未来自己在设计装 置时，可以发现可能的运动冲突。 编程计算得到角度的变化关系后，可以算出任 意时刻各铰的位置，以及 BC 杆上不同点的运动轨 迹 (图 3)：很明显 B 点轨迹是圆，C 点轨迹是圆的 一部分 (AB 杆大范围运动时，CD 杆只在小范围运 动)，而在 BC 杆上不同的点轨迹就很复杂了。 知道任意时刻系统的各铰位置，在屏幕上画出 各铰点的位置并连接起来，就得到了四连杆机构在 某一时刻的图象，延迟一定的时间后再画出下一时 刻的图象，就形成了动画。本问题中动画的源代码 见图 4，其中 plot 函数表示画线段；h1 是句柄，定义 A B D C 图 3 BC 杆上不同点的运动轨迹 图 4 动画部分的源代码截图 动画中每一帧图的特性，如颜色、线型、宽度等；set 函数是读取数据；drawnow 函数是画出当前帧的图 案；pause 函数表示暂停，让每一帧有一定的视觉残 留以便更好地形成动画 (类似放电影，1 秒钟放 24 格图像) [3]。 2 函数求导获得角速度、角加速度 求出角度后，对约束方程 (1) 求导出现角速度， 有 −a2φ˙ 1 sin φ1 − a3φ˙ 2 sin φ2 − a1ω0 sin θ = 0 a2φ˙ 1 cos φ1 + a3φ˙ 2 cos φ2 + a1ω0 cos θ = 0 } (5) 由于 θ, φ1, φ2 已在前面求出，因此得到关于 φ˙ 1, φ˙ 2 的一组线性方程组。类似 X=inv(A)*B 可解出角速 度，从而可以获得角速度随时间或随 θ 变化的关系 (图 5)。 对式 (5) 进一步求导得到角加速度的方程 −a2φ¨1 sin φ1 − a3φ¨2 sin φ2 − a2φ˙ 2 1 cos φ1− a3φ˙ 2 2 cos φ2 − a1ω 2 0 s cos θ = 0 a2φ¨1 cos φ1 + a3φ¨2 cos φ2 − a2φ˙ 2 1 sin φ1− a3φ˙ 2 2 sin φ2 − a1ω 2 0 sin θ = 0 (6) 此时 θ, φ1, φ2, φ˙ 1, φ˙ 2 都是已知值， 因此得到关于 φ¨1, φ¨2 的一组线性方程组。类似 X=inv(A)*B 可解 出角加速度，从而可以获得角加速度随时间或随 θ 变化的关系 (图 6)。 这只是一类典型机构的例子，采用本文介绍的 方法，任意机构的运动轨迹、刚体的角速度和角加速 度、某点的速度和加速度都可以类似求出。而这些结 果 (图 3、图 5、图 6) 都是传统方法难以获得的
第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 417 0 1 2 3 4 5 6 7 θ/rad -2 1 0 1 2 dθ/dt dϕ/dt dϕ/dt 䀂䙏ᓖ/(rad.s-1 ) 图 5 角速度与 θ 的关系 0 1 2 3 4 5 6 7 θ/rad d θ/dt d ϕ/dt d ϕ/dt 䀂࣐䙏ᓖ/(rad.s-2 ) -6 -4 -2 0 2 4 6 图 6 角加速度与 θ 的关系 3 Matlab 中的符号求导 如果运动学问题中没有具体数值，全是参数，用 Matlab 也可以处理。 案例 2：半径为 r 的圆轮在水平桌面上作直线 纯滚动，轮心速度 vO 的大小为常数。摇杆 AB 与 桌面铰接，并靠在圆轮上 (图 7)。当摇杆与桌面夹角 φ = 60◦ 时，试求摇杆的角速度和角加速度。 O vO A B ϕ 图 7 一类接触问题 传统处理方法：取 AB 杆为动系，O 为动点 (图 8)。通过速度和加速度分析 (具体过程略)，可 以得到 ωAB = vO 2r (顺时针), εAB = √ 3v 2 O/4r 2 (逆时针) (7) ve vO ωAB vr O A B 图 8 速度分析 Mtlab 处理本问题时，只需要把一般位置角度的 表达式写出来，然后利用符号求导的方法，就可以得 到角速度和角加速度。角度以逆时针方向为正，设初 始时刻 AB 杆垂直，根据图 9 有 x = r + vOt tan ( 1 2 φ ) = r/x (8) ωAB O r x A B ϕ 图 9 运动学关系 根据式 (8) 有 φ = 2 arctan[r/(r + vOt)] (9) 对式 (9) 求一阶导数获得角速度，求二阶导数获 得角加速度。程序的源代码见图 10
418 力 学 与 实 践 2021 年 第 43 卷 图 10 函数求导数的源代码 源代码中具体涉及几个函数：(1) 符号定义是 “syms”。(2) 函数求导的格式是 diff (function, ’variable’)，表示函数对某个变量求导。运行后得到的结 果为 ωAB = − 2rv0 2r 2 + v 2 0 t 2 + 2rv0t (10) εAB = 4rv2 0 (r + v0t) (2r 2 + v 2 0 t 2 + 2rv0t) 2 (11) 可以看到，传统方法的分析要很大的篇幅，而 Matlab 只需要几行就得到了结果 (顺便说一下，前 面式 (5) 和式 (6) 也可以利用符号推导从式 (1) 得 到)。 如果要求特定位置的值，只需要把具体值代入， 格式为 subs (function, old, new)，表示用新的变量替 换老的变量，这部分程序的源代码见图 11。 图 11 求特定位置的运动量 运行后屏幕显示结果为 when Φ = 60◦ dΦ/dt = −v0/(2 ∗ r) when Φ = 60◦ d2Φ/dt2 = (3∧(1/2) ∗ v0∧2)/(4 ∗ r ∧2) 注意负号表示顺时针方向，因此计算机推导出 的结果与传统方法相同，见式 (7)。 利用参数替换函数，还可以把式 (1) 和式 (11) 中的时间用角度替换 (具体略)。 因此可以画出系统运动全程的角度、角速度、角 加速度变化曲线，取归一化参数 r = 1，vO = 1，得 到的曲线关系见图 12 和图 13。 0 0.5 1.0 1.5 2.0 t/s φ/rad, ω/(rad.s-1), ε/(rad.s-2 ) -0.5 0 0.5 1.5 1.0 2.0 䀂ᓖφ 䀂䙏ᓖω 䀂࣐䙏ᓖε 图 12 系统运动参量与时间的关系 0.8 1.0 1.2 1.4 1.6 -0.6 -0.4 -0.2 0 0.2 0.4 䀂䙏ᓖω 䀂࣐䙏ᓖε φ/rad ω/(rad.s-1), ε/(rad.s-2 ) 图 13 系统运动参量与角度的关系 借助 Matlab 还可以对题目进行深入研究。例如， 用传统方法分析时，如果动点选为圆盘上的接触点， 速度还可能会正确，但是加速度没有办法分析了，这 是为什么呢？可以看看接触点相对 AB 杆的轨迹，这 借助计算机很容易实现。 以 φ = 60◦ 为系统初始位置，设圆盘上 P 点与 AB 杆接触，且 AP = √ 3r。当圆心运动 ∆x = vOt 时，圆盘因为纯滚动的转角为 θ = ∆x/r (图 14)，且 P 点在动系中有 xP = r(1 − cos θ) yP = √ 3r + ∆x + r sin θ = √ 3r + r(θ + sin θ) (12)
第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 419 在动系中画出的轨迹如图 15 所示，很像是旋轮 线。如果 P 点的轨迹是旋轮线，则曲线的尖点应该 接触 AB 杆；而实际上 P 点接触杆时，相对速度不 为零，但轨迹的切线方向与杆平行，不过该曲线的曲 率在详细分析前还不清楚。 yr Dx vO xr B A O P θ θ ϕ 图 14 接触点的轨迹分析 xr/m yr/m A ᧕䀖⛩䖘䘩 ABᵶ 图 15 接触点的轨迹 这就回答了前面的问题：如果选 P 点为动点， 在分析速度时认为相对速度沿杆方向 (相信这样做 的学生自己也不能说得很清楚，只是凭感觉)，碰巧 是正确的；但是在加速度分析时，由于相对运动的轨 迹不是太清楚，相对切向加速度未知，相对法向加速 度也未知，再加上 AB 杆的角加速度未知，未知数 太多求解不了。 注意在地面上看 P 点的轨迹是旋轮线，通过反 转和平移，旋轮线与图 15 中的轨迹完全相同，这个 结论是否出乎意料？ 4 小结 在运动学中引入 Matlab，可以方便了解系统整 个运动过程，获得系统中任意点的运动轨迹，以及速 度和加速度随时间或某角度的变化规律，直观明了， 让学生更容易理解系统的运动。 本篇着重介绍了非线性方程组的求解方法，以 及画图和动画的格式；另外介绍了参数方程的求导 和替换。画图涉及到的 Matlab 函数是 plot (画直 线)，而句柄函数、set 函数和 drawon 函数的组合就 可以实现动画；另一个重要函数是 diff (函数求导) 和 subs (变量替换)。掌握这几个命令后，就可以对 一般的运动学问题进行全程的计算和动画演示，获 得丰富的运动学信息。更进一步，借助 Matlab 可以 更加深入研究运动学的问题，解决传统分析方法无 法处理的很多问题。 参 考 文 献 1 李俊峰等. 理论力学. 北京: 清华大学出版社, 2003 2 高云峰. 理论力学辅导与习题集. 北京: 清华大学出版社, 2003 3 甘勤涛, 胡仁喜, 程政田等. Matlab 数学计算与工程分析 —— 从入门到精通. 北京: 机械工业出版社, 2019 (责任编辑: 胡 漫)