Matlab连杆机构运动学仿真：四杆与曲柄滑块gif动画

在机械原理课程里，四杆机构是最经典的运动学分析对象，但我见过很多同学卡在同一步：位置方程能列出来，手算也能算几个特殊位置，可是一旦要把整周运动画成位移曲线、看清连杆的真实轨迹，或者让机构真正“动起来”，就不知道下一步怎么办了。

这篇文章想给出一个明确判断：用 Matlab 做连杆机构运动学仿真的真正门槛，不在于编程，而在于把机构简图转化成可求解的数学方程。只要你掌握了闭环矢量方程、逐帧绘图、gif 动画输出这条主线，曲柄滑块、四杆、五杆、六杆在你眼里都会是同一道题的三种变形。读完本文，你可以独立复现曲柄滑块机构和铰链四杆机构的完整运动学仿真，并且生成动态 gif，为后续做六杆机构、机构优化和 Simscape 动力学分析打好底子。

1. 为什么用 Matlab 做连杆机构运动学仿真

做机构运动学仿真，可选方案很多：Adams、SolidWorks Motion、Simulink Simscape，都自带可视化。但 Matlab 纯脚本方案到目前为止仍然有不可替代的价值。

第一，它让你面对数学本身。在 Adams 或 SolidWorks 里，你拖几个约束、点一下仿真，就能看到动画，但内部的位置方程是怎么解出来的，对使用者来说是黑盒。Matlab 里你必须自己把机构几何关系写成方程，再自己解方程。这个过程看起来多绕了一步，实际上恰恰是理解机构运动学最快的一条路。

第二，它的自动化能力强。Matlab 的矩阵运算让“一个循环跑 200 个位置”变得极其自然。你不再需要手动取点，也不需要像 CAD 软件那样逐个位置做几何约束，而是用一段脚本批量计算整周运动。

第三，它的绘图和 gif 输出链路非常成熟。逐帧生成图片、用getframe捕获、再通过imwrite合成 gif，这套流程可以复用到任何仿真结果上，不局限于连杆机构，也可以用于齿轮啮合、凸轮轮廓、机械臂轨迹等项目。

所以，本文的定位不是“用 Matlab 替代专业机械仿真软件”，而是先用 Matlab 把机构运动分析的数学原理讲透。Computed Aided Engineering 工具再强大，也无法替代你对机构本身的理解。

2. 连杆机构运动学基础：从四杆到滑块的本质

2.1 运动学与动力学的边界

机构学里，运动学研究的是位置、速度、加速度与时间或输入角度的关系，不涉及力和质量。动力学则研究力、力矩、质量与运动之间的因果关系。

可以这样类比：运动学像录像分析，你只看运动员的肢体轨迹、速度变化；动力学像肌肉发力分析，你要研究哪些力产生了这些运动。连杆机构仿真通常先做运动学，得到运动规律后再做动力学。本文只覆盖运动学，但代码框架对动力学扩展同样有效。

2.2 平面机构的自由度

平面机构自由度用 Grübler-Kutzbach 公式计算：

F = 3(n - 1) - 2P_L - P_H

其中n是构件数，P_L是低副数（转动副、移动副），P_H是高副数。

• 四杆机构：n=4, P_L=4，得到F=1，只需要一个输入，就能确定整个机构运动。
• 曲柄滑块机构：n=4, P_L=4（三个转动副、一个移动副），自由度也是 1。
• 纯铰链五杆机构：n=5, P_L=5，F=2，需要两个输入才能确定运动，这就是五杆比四杆复杂的原因。
• 六杆机构：典型瓦特六杆n=6, P_L=7，F=1，但它的闭环更多，求解时往往要拆成多个闭环依次求解。

理解自由度是写仿真代码的第一步。自由度为 1 时，你遍历输入角度就能生成机构的完整运动；自由度为 2 时，你需要在两个输入之间建立某种约束，否则动画会“乱动”。

2.3 平面连杆机构的闭环矢量方程

平面连杆机构可以看成一系列矢量首尾相接形成的闭环。以铰链四杆机构为例，A 为固定铰，D 为固定铰，AB 是曲柄，BC 是连杆，CD 是摇杆，闭环矢量方程为：

AB + BC = AD + DC

写成复数形式：

a*e^(i*theta1) + b*e^(i*theta2) = d + c*e^(i*theta3)

实部和虚部分别相等，得到两个位置方程。已知输入角theta1，未知数是theta2和theta3，两个方程两个未知数，可解。

曲柄滑块机构其实是四杆机构的变体：当摇杆长度趋于无穷大时，摇杆末端的圆弧轨迹退化为直线，就得到滑块。因此曲柄滑块的解析解更简单，适合零基础入门。

2.4 Grashof 条件

不是任意四杆都能让曲柄整周回转。设四杆长度为a, b, c, d，其中s为最短杆，l为最长杆，如果满足：

s + l <= 其余两杆之和

并且最短杆为机架或与机架相邻，则存在整周回转构件。这就是 Grashof 条件。

• 最短杆为连架杆时，得到曲柄摇杆机构。
• 最短杆为机架时，得到双曲柄机构。
• 最短杆为连杆时，得到双摇杆机构。

写仿真代码之前先检查 Grashof 条件，可以避免在某个角度出现“根号内为负”或“arccos 越界”的尴尬错误。

3. 环境准备与仿真主流程

3.1 环境要求

本文代码全部使用 Matlab 基础函数，不需要 Simulink，也不需要额外的工具箱。理论上 R2016b 之后的版本都能运行。如果用 R2020a 之后的版本，可以使用exportgraphics输出 gif；本文默认采用兼容性最好的getframe + imwrite方案。

操作系统不限，Windows、macOS、Linux 均可。需要注意两点：

• 动画捕获getframe需要图形界面环境，不建议在完全 headless 的服务器上直接运行。
• 中文注释在部分旧版 Matlab 编辑器里可能乱码，建议源码文件统一保存为 UTF-8 编码，或者直接使用英文注释。

3.2 仿真主流程

后面的示例都遵守同一个主流程：

1. 确定机构拓扑和杆长参数。
2. 根据闭环矢量方程建立位置求解公式。
3. 遍历输入角度，逐帧求解未知角度（或滑块位移）。
4. 在需要时对时间求导，得到速度和加速度。
5. 每一帧重绘机构简图，捕获画面，写入 gif。

把这个流程记熟，后面的四杆、五杆、六杆只是位置方程更多、迭代更复杂，主线不变。

4. 示例一：曲柄滑块机构仿真（解析法 + gif）

4.1 数学模型

曲柄滑块机构中，设曲柄长度r，连杆长度l，曲柄转角theta，滑块位移x。几何关系为：

x = r*cos(theta) + sqrt(l^2 - r^2*sin(theta)^2)

对时间求导，得到滑块速度：

v = -r*omega*sin(theta) - (r^2*omega*sin(theta)*cos(theta)) / sqrt(l^2 - r^2*sin(theta)^2)

加速度可以直接继续求导，也可以像本文代码一样用gradient数值微分。解析式容易写错，数值微分适合验证。

4.2 完整代码

保存为slider_crank.m：

% 文件: slider_crank.m % 曲柄滑块机构运动学仿真，输出 gif 动画 clear; clc; close all; % 参数定义 r = 0.10; % 曲柄长度 m l = 0.30; % 连杆长度 m omega = 2*pi; % 曲柄角速度 rad/s，约每秒一转 N = 200; % 一个周期采样点数 theta = linspace(0, 2*pi, N); % 曲柄转角 t = theta / omega; % 对应时间 dt = t(2) - t(1); % 时间步长 % 位置解析解 x = r*cos(theta) + sqrt(l^2 - r^2*sin(theta).^2); % 速度解析解 v = -r*omega*sin(theta) ... - (r^2*omega*sin(theta).*cos(theta)) ./ sqrt(l^2 - r^2*sin(theta).^2); % 加速度：数值微分验证 a = gradient(v, dt); % 创建画布 figure('Position', [100 100 900 420]); for i = 1:N % 左图：机构动画 subplot(1, 2, 1); cla; hold on; axis equal; xlim([-0.45 0.45]); ylim([-0.35 0.35]); xA = 0; yA = 0; % 曲柄固定铰 xB = r*cos(theta(i)); yB = r*sin(theta(i)); % 曲柄与连杆连接点 xC = x(i); yC = 0; % 滑块位置 % 曲柄 plot([xA xB], [yA yB], 'b-o', 'LineWidth', 2); % 连杆 plot([xB xC], [yB yC], 'r-o', 'LineWidth', 2); % 滑块（矩形） rect_pos = [xC-0.02 -0.02; xC+0.02 -0.02; xC+0.02 0.02; xC-0.02 0.02]; patch('Vertices', rect_pos, 'Faces', [1 2 3 4], ... 'FaceColor', [0.7 0.7 0.7], 'EdgeColor', 'k'); % 滑道 plot([-0.45 xC], [0 0], 'k--'); % 固定铰标记 plot(xA, yA, 'ko', 'MarkerFaceColor', 'k'); title(sprintf('t = %.3f s', t(i))); xlabel('x (m)'); ylabel('y (m)'); % 右图：滑块位移曲线 subplot(1, 2, 2); plot(t(1:i), x(1:i), 'b-', 'LineWidth', 1.5); xlim([0 t(end)]); ylim([min(x) max(x)]); xlabel('时间 (s)'); ylabel('滑块位移 x (m)'); title('位移曲线'); grid on; % 捕获当前帧，写入 gif drawnow; frame = getframe(gcf); im = frame2im(frame); [imind, cm] = rgb2ind(im, 256); if i == 1 imwrite(imind, cm, 'slider_crank.gif', 'gif', ... 'Loopcount', inf, 'DelayTime', 0.03); else imwrite(imind, cm, 'slider_crank.gif', 'gif', ... 'WriteMode', 'append', 'DelayTime', 0.03); end end disp('动画已保存为 slider_crank.gif');

4.3 代码逻辑说明

位置解只用了一行x = r*cos(theta) + sqrt(l^2 - r^2*sin(theta).^2)，这就是解析法的好处。速度求导容易出错，所以代码里保留了完整的解析表达式。加速度用gradient(v, dt)做数值微分，是一种快速校验手段，如果解析速度公式写错，加速度曲线会出现明显跳变。

动画部分，我用了cla清空坐标区再重新绘制。这种方式简单直接，缺点是性能一般。如果采样点数较大，可以改用先创建图形对象、再更新XData/YData的方式，后面第 9 节会说明。

4.4 运行结果与验证

直接运行脚本，会在当前目录生成slider_crank.gif。用浏览器或图片查看器打开，可以看到曲柄带动连杆推动滑块往复运动。

从运动学角度验证几个特征：

• 滑块位移范围应该在[l-r, l+r] = [0.2, 0.4]之间。
• 滑块速度在位移中点附近最大，在两端接近 0。
• 一个周期内滑块往复一次，位移曲线与单缸发动机活塞运动规律一致。

如果位移范围不对，优先检查r和l的定义；如果 gif 文件没有正常生成，检查当前目录是否可写、循环里是否匹配了if i == 1和else分支。

5. 示例二：铰链四杆机构仿真（数值解 + 轨迹绘制）

5.1 数学模型

铰链四杆机构的位置求解比曲柄滑块复杂一些，因为theta2和theta3不能直接显式解出，但可以通过代数消元得到解析表达式。

设固定铰 A 为(0,0)，D 为(d,0)。曲柄 AB 长度为a，连杆 BC 长度为b，摇杆 CD 长度为c。由矢量闭环：

B = (a*cos(theta1), a*sin(theta1))

C 点既要满足C = B + b*(cos(theta2), sin(theta2))，又要满足|C - D| = c。展开后得到：

A*cos(theta2) + B*sin(theta2) = K

其中：

A = xB - d B = yB K = (c^2 - (xB-d)^2 - yB^2 - b^2) / (2*b)

两边同除R = sqrt(A^2 + B^2)，化为：

cos(theta2 - phi) = K / R phi = atan2(B, A)

于是：

theta2 = phi ± acos(K / R)

两个解对应机构的两种装配模式，通常称为“开式”和“交叉式”。在动画中必须保证每一帧选择同一个装配模式，否则机构会出现瞬间翻转。最简单的策略是让当前帧的解与上一帧的角度差值最小。

求出theta2后，C 点坐标已知，摇杆角度：

theta3 = atan2(yC, xC - d)

5.2 完整代码

保存为four_bar.m：

% 文件: four_bar.m % 铰链四杆机构运动学仿真，输出 gif 动画和摇杆摆角曲线 clear; clc; close all; % 杆长参数 a = 0.10; % 曲柄 AB b = 0.30; % 连杆 BC c = 0.25; % 摇杆 CD d = 0.28; % 机架 AD % Grashof 条件检查 s = min([a b c d]); l = max([a b c d]); rest = sum([a b c d]) - s - l; if s + l <= rest + 1e-10 disp('满足 Grashof 条件，曲柄可以整周回转。'); else disp('不满足 Grashof 条件，机构可能存在装配死角。'); end omega = 2*pi; % 曲柄角速度 N = 200; theta1 = linspace(0, 2*pi, N); t = theta1 / omega; dt = t(2) - t(1); theta2 = zeros(1, N); theta3 = zeros(1, N); xC_all = zeros(1, N); yC_all = zeros(1, N); % 初始猜测，选择开式装配模式 theta2_prev = 0.5; for i = 1:N xB = a * cos(theta1(i)); yB = a * sin(theta1(i)); A = xB - d; B = yB; K = (c^2 - (xB-d)^2 - yB^2 - b^2) / (2*b); R = sqrt(A^2 + B^2); phi = atan2(B, A); % 判断装配是否可达 if abs(K / R) > 1 error('theta1 = %.2f 时机构无法装配，请检查杆长参数', theta1(i)); end alpha = acos(K / R); cand1 = phi + alpha; cand2 = phi - alpha; % 选择与上一帧角度差值最小的解，保持装配模式连续 [~, idx] = min(abs([cand1 cand2] - theta2_prev)); cand = [cand1 cand2]; theta2(i) = cand(idx); theta2_prev = theta2(i); % 由 theta2 计算 C 点坐标 xC_all(i) = xB + b*cos(theta2(i)); yC_all(i) = yB + b*sin(theta2(i)); % 摇杆角度（DC 从 D 指向 C） theta3(i) = atan2(yC_all(i), xC_all(i) - d); end % 角速度（数值微分） omega3 = gradient(theta3, dt); % 绘制动画 figure('Position', [100 100 950 430]); for i = 1:N % 左图：机构动画 subplot(1, 2, 1); cla; hold on; axis equal; xlim([-0.15 0.55]); ylim([-0.35 0.35]); xB = a*cos(theta1(i)); yB = a*sin(theta1(i)); xC = xC_all(i); yC = yC_all(i); % 固定铰 plot(0, 0, 'ko', 'MarkerFaceColor', 'k'); plot(d, 0, 'ko', 'MarkerFaceColor', 'k'); % 曲柄 AB plot([0 xB], [0 yB], 'b-o', 'LineWidth', 2); % 连杆 BC plot([xB xC], [yB yC], 'r-o', 'LineWidth', 2); % 摇杆 CD plot([d xC], [0 yC], 'g-o', 'LineWidth', 2); % C 点运动轨迹 plot(xC_all(1:i), yC_all(1:i), 'm.', 'MarkerSize', 4); title(sprintf('t = %.3f s', t(i))); xlabel('x (m)'); ylabel('y (m)'); % 右图：摇杆摆角和角速度 subplot(1, 2, 2); yyaxis left; plot(t(1:i), theta3(1:i)*180/pi, 'g-', 'LineWidth', 1.5); ylabel('摇杆摆角 (deg)'); yyaxis right; plot(t(1:i), omega3(1:i), 'r--', 'LineWidth', 1.2); ylabel('摇杆角速度 (rad/s)'); xlabel('时间 (s)'); xlim([0 t(end)]); grid on; title('摇杆运动曲线'); drawnow; frame = getframe(gcf); im = frame2im(frame); [imind, cm] = rgb2ind(im, 256); if i == 1 imwrite(imind, cm, 'four_bar.gif', 'gif', ... 'Loopcount', inf, 'DelayTime', 0.03); else imwrite(imind, cm, 'four_bar.gif', 'gif', ... 'WriteMode', 'append', 'DelayTime', 0.03); end end disp('动画已保存为 four_bar.gif');

5.3 运行结果与验证

运行后生成four_bar.gif。左图是机构动画，C 点会画出一条闭合的连杆曲线，这条曲线在机械原理中叫“连杆曲线”，是四杆机构最重要的输出轨迹之一。右图同时显示摇杆摆角与摇杆角速度。

验证要点：

• 曲柄转角从 0 到 360 度变化时，摇杆摆动角度应该是连续周期变化，不会出现突变。
• 摇杆角速度曲线在换向点附近接近 0，符合实际物理规律。
• 如果动画在某帧突然翻转，说明装配模式选择逻辑失效，应检查theta2_prev的初始值以及角度连续性判断。

这段代码展示了一个重要思路：当位置方程不能直接显式求解时，先消元化成A*cos(theta2) + B*sin(theta2) = K的形式，再求解。这个思路比直接调用fsolve更稳定，也不需要优化工具箱，适合零基础复现。

6. 通用 gif 输出封装与动画绘制技巧

前面两个示例都重复了一段 gif 写入代码。在实际项目中，建议封装成独立函数，避免每个仿真脚本都复制一遍。

保存为save_gif.m：

function save_gif(fig, filename, frame_idx, delay) % 将当前 figure 保存为 gif % fig: figure 句柄 % filename: 输出文件名 % frame_idx: 当前帧序号，第 1 帧创建文件，后续帧追加 % delay: 帧间延迟，单位秒 frame = getframe(fig); im = frame2im(frame); [imind, cm] = rgb2ind(im, 256); if frame_idx == 1 imwrite(imind, cm, filename, 'gif', ... 'Loopcount', inf, 'DelayTime', delay); else imwrite(imind, cm, filename, 'gif', ... 'WriteMode', 'append', 'DelayTime', delay); end end

调用方式：

% 在动画循环内 save_gif(gcf, 'my_sim.gif', i, 0.03);

如果你的 Matlab 版本是 R2020a 及以上，还可以使用更简洁的exportgraphics：

exportgraphics(gcf, 'my_sim.gif', 'Append', i ~= 1, 'Resolution', 100);

两种方式的区别：

动画绘制的几个通用技巧：

• 每一帧都要固定xlim和ylim，否则画面会随机构位置缩放而抖动。
• axis equal要放在xlim之前，避免坐标轴比例失真。
• DelayTime适合设为 0.02 到 0.1 秒。太小时动画闪动，太大时看起来卡顿。
• 如果一帧里既要画机构，又要画曲线，优先用subplot分开展示。

7. 从四杆到五杆、六杆：通用化建模思路

7.1 五杆机构自由度是关键

纯铰链五杆机构的自由度是 2，因此严格意义上不能像四杆机构那样只给一个曲柄输入就得到确定运动。实际工程中的单自由度五杆机构，通常通过以下方式实现：

• 引入一个移动副，例如五杆加滑块，自由度降为 1。
• 两个输入之间增加齿轮或带传动约束，形成“齿轮五杆机构”。
• 让两个输入保持固定比例关系，例如曲柄同步旋转。

仿真这类机构时，闭合矢量方程仍然是核心，只是未知角度从 2 个变成 3 个或更多，需要把多个环路方程联立成一个方程组，再用牛顿-拉夫森迭代求解。

7.2 六杆机构的闭环拆分

典型六杆机构如瓦特六杆、史蒂芬森六杆，通常由两个闭环组成。求解策略有两种：

• 顺序求解：先解第一个四杆环，再把结果作为第二个环的已知条件。编程简单，但不是所有六杆都能拆成标准四杆。
• 联立求解：把所有闭环的位置方程写成F(X) = 0的形式，用牛顿迭代一次性求解所有未知角度。通用性强，但需要给出合适的初值。

联立求解的伪代码框架如下：

% 伪代码：多闭环牛顿-拉夫森迭代 % X 为未知角度向量，例如 [theta2; theta3; theta5; theta6] % F(X) 为位置残差向量，每个闭环提供实部、虚部两个方程 function X = solve_mechanism(X0, params, tol) X = X0; for k = 1:100 F = closed_loop_residual(X, params); J = closed_loop_jacobian(X, params); dX = J \ (-F); X = X + dX; if norm(dX) < tol break; end end end

这里的核心工作量在求 Jacobian。如果手动推导太繁琐，可以用 MATLAB 的符号工具箱生成 Jacobian 的解析表达式，再转成数值函数，这是一条很实用的工程路径。

7.3 从四杆到六杆本质没有变

四杆机构用“已知一个角度，求两个角度”的消元法；五杆、六杆用“已知若干输入，求多个未知角度”的牛顿迭代。数学形式从二维扩展到多维，但思路一致：列出闭环矢量方程，实部虚部分别为等式，解非线性方程组，逐帧更新，输出动画。这就是我在文章开头说的“同一道题的三种变形”。

8. 常见问题与排查方法

仿真代码跑不通时，优先看模型问题而不是语法问题。下面是连杆机构仿真中最常见的几类问题。