又到了机械原理课程设计答辩季。很多同学拿到“平面六杆机构运动分析”这类题目后，第一反应是翻教材找图解法：按比例画速度多边形、加速度多边形，画完一张 A3 图纸；导师说“杆长改一下”，前面所有图全部重来。用 MATLAB 做六杆机构仿真，其实可以把这件事从“手工画图”变成“参数化计算”：改一个杆长，程序重新跑一遍，位置、速度、加速度曲线全部更新，整套数据和图表都能直接用于课程设计报告。
这篇文章不是简单贴一段“能跑的代码”，而是把六杆机构运动学分析背后的建模思路、MATLAB 数值解法、结果验证方法和常见坑一次讲清楚。无论你是机械设计制造及其自动化专业的学生，还是正在准备考研复试、机械创新竞赛，这篇文章都能帮你少走弯路。
本文以一个典型的“曲柄 AB + 连杆 BC + 摇杆 CD + 连杆 EF + 滑块 F”平面六杆滑块机构为例，从闭环矢量法建模开始，逐步给出完整 MATLAB 代码，并说明如何验证结果、如何排查迭代不收敛、角度跳变等常见问题。
1. 这篇文章真正要解决的问题
机械原理课程设计里，机构运动分析是必须完成的基础内容。传统做法是图解法和解析法并行：
- 图解法直观，但精度取决于作图比例尺和手工测量，改参数就要重新画；
- 解析法精度高，但四杆以上机构的手推公式会让人崩溃，尤其是速度、加速度方程，二阶求导出错概率很高；
- 用通用仿真软件，比如 ADAMS、Simulink Simscape Multibody，学习成本高，很多同学来不及在课设周期内掌握。
MATLAB 的价值在于：它既有数值计算能力，又有绘图能力，还能导出报告需要的图片和数据。六杆机构运动分析真正落地，靠的是一套“位置约束方程 + Jacobian 矩阵”框架：同一个 Jacobian 矩阵，既能做牛顿迭代求位置，又能联立求速度和加速度，代码量少、扩展性好。
读这篇文章，你最终能收获四件事：
- 理解六杆机构自由度计算和运动学建模思路；
- 掌握闭环矢量法建立位置约束方程的方法；
- 拿到一套完整可运行的 MATLAB 代码，直接改参数就能用；
- 知道怎么验证程序算对了，以及答辩时怎么解释关键问题。
下表是三种常见做法的对比：
| 做法 | 精度 | 改参数成本 | 可视化能力 | 学习成本 |
|---|---|---|---|---|
| 图解法 | 低 | 高 | 中 | 低 |
| 手推解析法 | 高 | 高 | 低 | 高 |
| MATLAB 数值法 | 高 | 低 | 高 | 中 |
2. 六杆机构运动学基础
2.1 什么是平面六杆机构
平面六杆机构指由六个构件（包含机架）通过平面低副连接而成的机构。相比四杆机构，六杆机构能实现更复杂的运动规律，比如间歇运动、急回特性、大行程滑块运动等。工程中常见的瓦特六杆机构、史蒂芬森六杆机构，以及本文这种“四杆机构 + 二级杆组”的组合式六杆机构，都属于这一范畴。
本文选取的机构可以描述为：
- 构件 1：曲柄 AB，绕固定铰链 A 匀速转动；
- 构件 2：连杆 BC；
- 构件 3：摇杆 CD，绕固定铰链 D 摆动；
- 构件 4：连杆 EF；
- 构件 5：滑块 F，沿水平导轨往复移动；
- 机架：固定支撑 A、D 和导轨。
其中 E 点不是独立构件，而是 BC 杆上的固定点位，通常用比例系数 m 表示，即 E = B + m × (C − B)。m=0.5 时 E 位于 BC 中点。EF 杆一端在 E 点与 BC 杆铰接，另一端与滑块 F 铰接。
2.2 自由度计算
平面机构自由度公式为：
F = 3n − 2PL − PH
其中 n 为活动构件数，PL 为低副数，PH 为高副数。本文机构活动构件数为 5（曲柄、连杆、摇杆、连杆 EF、滑块），低副数量为 7（A、B、C、D、E、F 转动副 + 滑块与导轨的移动副），没有高副。则：
F = 3×5 − 2×7 = 1
自由度等于 1，意味着给曲柄一个独立运动参数（转角 θ1），机构所有构件的运动就完全确定。这也是课程设计中“给一个原动件，整机运动确定”的标准表述。
2.3 运动学分析的三大任务
机构运动分析包括三类问题：
- 位置分析：已知曲柄转角 θ1，求各杆转角（θ2、θ3、θ4）和滑块位移 xF；
- 速度分析：求各杆角速度（ω2、ω3、ω4）和滑块速度 vF；
- 加速度分析：求各杆角加速度（α2、α3、α4）和滑块加速度 aF。
位置分析是基础，因为速度和加速度都需要在已知位置的前提下求解。加速度分析尤为重要，后续做惯性力分析、动态静力分析时，都要用到构件的角加速度和质心加速度。
3. 数学建模：闭环矢量法
3.1 为什么选择闭环矢量法
机械原理中有多种解析方法：复数法、矩阵法、杆组法。闭环矢量法的核心思想是：把每个构件看成矢量，首尾相接形成封闭环，然后投影到 x、y 方向得到约束方程。这种方法思路统一，适合用 MATLAB 数值求解，而且方便扩展到八杆、十杆机构。
3.2 建立位置约束方程
以固定铰链 A 为原点，x 轴水平向右，y 轴竖直向上。机构中 A、D 为固定点，其中 A = (0, 0)，D = (AD, 0)。设曲柄 AB 与 x 轴夹角为 θ1，连杆 BC 与 x 轴夹角为 θ2，摇杆 CD 与 x 轴夹角为 θ3，连杆 EF 与 x 轴夹角为 θ4，滑块 F 的水平坐标为 xF，导轨 y 坐标为 yF。
先从四杆机构 ABCD 的闭环出发，得到两个方程：
$$ l_1\cos\theta_1 + l_2\cos\theta_2 - l_3\cos\theta_3 - AD = 0 $$
$$ l_1\sin\theta_1 + l_2\sin\theta_2 - l_3\sin\theta_3 = 0 $$
再计算 E 点坐标。E 点在 BC 上，满足：
$$ E_x = l_1\cos\theta_1 + m, l_2\cos\theta_2 $$
$$ E_y = l_1\sin\theta_1 + m, l_2\sin\theta_2 $$
最后，EF 杆和滑块 F 形成第二组约束：
$$ E_x + l_4\cos\theta_4 - x_F = 0 $$
$$ E_y + l_4\sin\theta_4 - y_F = 0 $$
把四个方程整理成 G(z, θ1) = 0，其中未知向量 z = [θ2, θ3, θ4, xF]^T。这样就完成了位置分析的数学建模。
3.3 Jacobian 矩阵与速度、加速度方程
对 G(z, θ1) = 0 关于时间求导，可以得到速度方程：
J × ż = b_v
其中 J = ∂G/∂z 是 Jacobian 矩阵。对本机构，J 的解析形式为：
$$ J = \begin{bmatrix} -l_2 \sin\theta_2 & l_3 \sin\theta_3 & 0 & 0 \ l_2 \cos\theta_2 & -l_3 \cos\theta_3 & 0 & 0 \ -m l_2 \sin\theta_2 & 0 & -l_4 \sin\theta_4 & -1 \ m l_2 \cos\theta_2 & 0 & l_4 \cos\theta_4 & 0
