1. 从玩具到数学万花尺的工程化仿真之旅小时候很多人都玩过一种叫“万花尺”的绘图玩具一个大尺环套着一个小齿轮齿轮上有孔用笔插进孔里推动小齿轮沿着大尺环的内壁滚动就能画出各种繁复而美丽的曲线。这种看似简单的玩具背后隐藏的数学原理却相当迷人——它就是圆内旋轮线。作为一名长期与MATLAB打交道的工程师我最近把儿时的玩具搬进了数字世界用MATLAB完整仿真了万花尺的绘制过程并深入到了“异形齿轮”的拓展领域。这不仅仅是一次怀旧更是一次将几何、运动学和数值计算融会贯通的绝佳实践。对于机械设计、动画生成乃至艺术创作领域的朋友来说理解如何用代码“驱动”一支虚拟的笔画出确定性的复杂图案具有很高的参考价值。本文将手把手带你复现这个过程并探讨如何将其原理应用于非圆齿轮的仿真让你不仅得到漂亮的图更能透彻理解每一行代码背后的“为什么”。2. 核心原理拆解圆内旋轮线的参数方程推导万花尺画出的曲线在数学上称为“圆内旋轮线”或“内摆线”。它的生成可以抽象为一个经典的几何运动模型一个半径为r的小圆在一个半径为R的大圆内部无滑动地滚动。我们关注的是固定在小圆上、距离小圆圆心为d的一个点P的运动轨迹。这个点P就对应万花尺上笔尖的位置。2.1 建立数学模型与坐标系为了用代码描述这个运动我们首先需要建立参数方程。设大圆圆心为原点O(0,0)。假设小圆初始时刻位于大圆的最右侧接触点小圆圆心C位于(R-r, 0)笔尖点P位于(R-rd, 0)。当小圆在大圆内逆时针滚动时它滚过的弧长等于大圆内壁被滚过的弧长无滑动条件。设小圆相对于自身圆心转过的角度为θ这是我们的核心运动参数那么小圆圆心C绕大圆圆心O转过的角度φ满足关系r * θ R * φ。由此可得φ (r / R) * θ。2.2 推导笔尖点P的坐标现在我们需要在固定的世界坐标系中求出点P的坐标。这个过程可以分解为两个变换的叠加小圆圆心C的位置C点绕原点O作圆周运动其半径为(R - r)当前角度为φ。因此C的坐标为Cx (R - r) * cos(φ)Cy (R - r) * sin(φ)笔尖P相对于圆心C的位置在小圆自身的坐标系中笔尖P位于距离圆心C为d的位置。由于小圆自身旋转了θ角度并且是在大圆内“滚动”而非“自转”其相对于大圆切向的方向需要仔细考虑。从C点指向P点的向量其方向角并不是简单的θ。因为小圆在滚动时其与地面的接触点速度为零这意味着小圆上每一点的运动是滚动公转和滑动自转的合成。正确的处理方式是P点相对于圆心C的位置是由小圆自身的旋转θ决定的但这个旋转是发生在随着圆心C移动的动坐标系中。并且由于滚动方向小圆在大圆内逆时针滚动与自转方向相反从接触点看P点在小圆坐标系中的角度应为-θ此处符号取决于坐标系定义是推导的关键。综合以上点P的坐标最终参数方程为设θ为自变量φ (r / R) * θ; Px (R - r) * cos(φ) d * cos(φ - θ); Py (R - r) * sin(φ) d * sin(φ - θ);其中(φ - θ)项描述了笔尖相对于动坐标系随小圆圆心平动并旋转的坐标系的方向。这个方程是仿真的基石。注意方程中的符号φ - θ还是θ - φ直接影响了曲线的旋向和形态。我强烈建议在首次实现时通过绘制小圆圆心C和笔尖P的瞬时位置向量来验证几何关系的正确性。一个快速的验证方法是当d 0时P点应与C点重合轨迹是一个圆当d r时P点应位于小圆周上应能画出标准的内摆线。3. MATLAB仿真实现从方程到动态绘图理解了原理用MATLAB实现就变得清晰了。我们的目标是1画出最终的静态曲线2动态展示绘制过程重现万花尺的运作。3.1 基础静态曲线绘制首先我们定义参数并计算轨迹点。这里的关键是参数θ的采样。为了让曲线闭合且光滑小圆需要滚动足够多的圈数直到笔尖回到起点。这要求R/r是一个有理数这样小圆滚动若干圈后公转角度φ将是2π的整数倍。我们通常取R和r为整数比。% 万花尺圆内旋轮线基础仿真 clear; clc; close all; % 1. 定义参数 R 100; % 大圆半径 r 31; % 小圆半径 d 65; % 笔尖到小圆圆心的距离 n_cycles 10; % 小圆滚动圈数根据R/r比值决定实际圈数 % 2. 计算采样点 % 小圆自转的总角度 圈数 * 2π theta_total n_cycles * 2 * pi; % 生成theta数组点数足够多以保证曲线光滑 num_points 5000; theta linspace(0, theta_total, num_points); % 3. 应用参数方程计算轨迹 phi (r / R) * theta; % 小圆圆心公转角度 Px (R - r) * cos(phi) d * cos(phi - theta); Py (R - r) * sin(phi) d * sin(phi - theta); % 4. 绘制静态结果 figure(Position, [100, 100, 800, 800]); plot(Px, Py, b-, LineWidth, 1.5); axis equal; grid on; box on; title(sprintf(圆内旋轮线 (R%d, r%d, d%d), R, r, d)); xlabel(X); ylabel(Y);运行这段代码你就能立即得到一幅精美的旋轮线图案。调整R,r,d三个参数你会发现图案千变万化当d r时笔尖在小圆内部曲线相对平滑。当d r时笔尖恰在小圆圆周上画出经典的内摆线可能有尖点。当d r时笔尖在小圆之外曲线会形成 loops环图案更为复杂。R与r的比值决定了图案的“瓣数”或对称性。若比值为整数n则曲线有n个尖点或n重对称性若为分数图案可能不闭合或周期极长。3.2 加入动态仿真与绘图过程静态图很美但动态过程更能揭示其生成机理。我们可以利用MATLAB的循环和drawnow命令来制作动画。% 动态绘制万花尺过程 figure(Position, [100, 100, 800, 800]); hold on; axis equal; grid on; box on; xlim([-R*1.2, R*1.2]); ylim([-R*1.2, R*1.2]); title(万花尺绘制动态仿真); xlabel(X); ylabel(Y); % 预先计算所有轨迹点同上 % ... % 初始化绘图对象 h_trace plot(Px(1), Py(1), b-, LineWidth, 1.5); % 轨迹线 h_point plot(Px(1), Py(1), ro, MarkerSize, 8, MarkerFaceColor, r); % 当前笔尖位置 h_circle_big rectangle(Position, [-R, -R, 2*R, 2*R], ... Curvature, [1, 1], EdgeColor, k, LineWidth, 1, LineStyle, --); h_circle_small plot(0, 0, g--, LineWidth, 1); % 小圆轮廓稍后更新 % 动态绘制循环 for i 1:50:num_points % 跳步绘制以提高动画速度 % 更新轨迹线绘制到当前点 set(h_trace, XData, Px(1:i), YData, Py(1:i)); % 更新笔尖位置 set(h_point, XData, Px(i), YData, Py(i)); % 计算并更新小圆轮廓 current_phi phi(i); current_theta theta(i); Cx (R - r) * cos(current_phi); Cy (R - r) * sin(current_phi); % 生成小圆上的点 alpha linspace(0, 2*pi, 100); small_circle_x Cx r * cos(alpha); small_circle_y Cy r * sin(alpha); set(h_circle_small, XData, small_circle_x, YData, small_circle_y); % 刷新图形 drawnow; % 可加入 pause(0.01) 控制速度 end hold off;这个动态仿真就像在电脑里运行了一个虚拟的万花尺。你可以清晰地看到小圆如何在大圆内滚动笔尖如何一步步勾勒出复杂图案。这对于教学和理解运动合成非常有帮助。3.3 参数探索与图案库生成万花尺的魅力在于参数微调就能产生截然不同的图案。我们可以写一个简单的脚本批量生成不同参数下的曲线建立一个“数字万花尺图案库”。% 参数探索固定R和d变化r R 100; d 60; r_values [19, 29, 33, 41, 47, 51]; % 尝试不同的r值 figure(Position, [50, 50, 1200, 600]); for idx 1:length(r_values) r r_values(idx); % 计算轨迹确保闭合 n lcm(R, r) / r; % 计算最小公倍数除以r得到所需圈数近似 theta linspace(0, n*2*pi, 3000); phi (r / R) * theta; Px (R - r) * cos(phi) d * cos(phi - theta); Py (R - r) * sin(phi) d * sin(phi - theta); subplot(2, 3, idx); plot(Px, Py, LineWidth, 1); axis equal; axis off; title(sprintf(r %d, r)); end sgtitle(不同小圆半径r下的万花尺图案 (R100, d60));通过这样的批量实验你可以直观地感受到数学的魔力并可能发现一些特别漂亮的参数组合。这本身就是一种创造性的探索。4. 进阶挑战异形齿轮非圆齿轮的万花尺仿真传统的万花尺基于两个圆形齿轮内啮合。那么如果大尺环或小齿轮不是标准的圆形呢这就是“异形齿轮”的范畴。仿真非圆齿轮的万花尺核心挑战在于无滑动滚动条件变得复杂因为齿轮的节曲线相当于等效滚动圆半径不再是常数。4.1 非圆齿轮节曲线与运动关系假设我们有一个椭圆形的外齿轮大尺环和一个圆形的小齿轮。椭圆没有恒定的曲率半径因此小圆圆心到接触点的距离即瞬时滚动半径是变化的。设椭圆参数方程为x_e a * cos(t), y_e b * sin(t)其中a b。 椭圆上某一点t的曲率半径ρ(t)公式为ρ(t) ( (a^2 * sin^2(t) b^2 * cos^2(t))^(3/2) ) / (a * b)在无滑动滚动条件下小圆半径r的圆心C不再沿一个固定半径的圆运动而是沿着椭圆的一条等距线运动该等距线与椭圆的距离为(R_eff - r)但这里的R_eff是时变的等于椭圆的曲率半径ρ(t)。这使得运动关系从一个简单的线性比例φ (r/R)*θ变成了一个需要数值积分求解的微分关系。具体来说设椭圆上的接触点参数为t小圆自转角度为θ。无滑动条件意味着小圆滚过的弧长等于椭圆上被滚过的弧长r * dθ ds其中ds是椭圆的弧长微分ds sqrt( (dx_e/dt)^2 (dy_e/dt)^2 ) * dt。 因此有关系式dθ (1/r) * ds。 同时小圆圆心C的位置是椭圆上该接触点沿内法线方向偏移(ρ(t) - r)的距离假设小圆在椭圆内部滚动。4.2 MATLAB数值仿真实现由于解析解很难获得我们采用数值方法逐步积分求解。% 异形齿轮椭圆外齿轮 圆形小齿轮万花尺仿真 clear; clc; close all; % 1. 定义椭圆参数和小圆参数 a 80; % 椭圆长半轴 b 50; % 椭圆短半轴 r 15; % 小圆半径 d 25; % 笔尖距离小圆圆心距离 % 2. 定义椭圆上的参数t数组对应接触点 num_steps 5000; t linspace(0, 2*pi, num_steps); % 椭圆参数角从0到2π % 3. 预计算椭圆几何量 x_e a * cos(t); y_e b * sin(t); dx_dt -a * sin(t); % x对t的导数 dy_dt b * cos(t); % y对t的导数 ds_dt sqrt(dx_dt.^2 dy_dt.^2); % 弧长微分 ds/dt % 计算曲率半径 ρ(t) numerator (a^2 * sin(t).^2 b^2 * cos(t).^2).^(3/2); denominator a * b; rho_t numerator ./ denominator; % 曲率半径 % 4. 数值积分求解小圆自转角度 theta theta zeros(size(t)); for i 2:num_steps % 弧长微元 ds ds ds_dt(i) * (t(i) - t(i-1)); % 根据无滑动条件r * dθ ds d_theta ds / r; theta(i) theta(i-1) d_theta; end % 5. 计算小圆圆心C的轨迹 % 椭圆上点(x_e, y_e)处的单位内法向量 N (-dy_dt, dx_dt) / norm % 注意方向要保证是“内”法线对于参数化方程 (a*cos(t), b*sin(t))内法线方向需根据实际情况判断。 % 一个简单方法是对于标准椭圆从原点指向椭圆上点的向量是 (x_e, y_e)其内法线大致指向椭圆中心。 % 更稳健的方法是计算单位切向量 T (dx_dt, dy_dt)/ds_dt然后通过旋转90度得到法向量。 % 这里我们计算单位切向量然后顺时针旋转90度得到内法向量假设逆时针参数化。 T [dx_dt./ds_dt; dy_dt./ds_dt]; % 2 x N 矩阵 % 顺时针旋转90度: (x, y) - (y, -x) N [T(2, :); -T(1, :)]; % 这给出了近似内法向量 % 圆心C 椭圆上的点 (ρ(t) - r) * 内法向量 Cx x_e (rho_t - r) .* N(1, :); Cy y_e (rho_t - r) .* N(2, :); % 6. 计算笔尖P的轨迹 % 笔尖相对于圆心C的向量长度为d方向角需要考虑小圆自身的旋转theta % 在动坐标系随C点平动且随小圆旋转中笔尖的方向角为 -theta滚动方向导致 % 需要将这个向量转换到世界坐标系。动坐标系的x轴方向是椭圆切向方向Ty轴方向是法向N。 % 因此笔尖在世界坐标系的偏移量为 d * [cos(-theta); sin(-theta)] 在 (T, N) 基下的坐标 % 即 offset d * (cos(theta).*T (-sin(theta)).*N) 这里需要仔细推导。 % 更直接的方法笔尖相对于C的位置在小圆自身的固定坐标系中角度为 -theta。 % 但小圆自身的固定坐标系相对于世界坐标系的旋转角度是多少这等于接触点切线与世界坐标系x轴的夹角。 % 设切向角为 alpha atan2(dy_dt, dx_dt) alpha atan2(dy_dt, dx_dt); % 那么笔尖P的坐标为 Px Cx d * cos(alpha - theta); Py Cy d * sin(alpha - theta); % 7. 绘图 figure(Position, [100, 100, 900, 900]); subplot(2,2,1); plot(x_e, y_e, k--, LineWidth, 1); hold on; plot(Cx, Cy, g-, LineWidth, 1); axis equal; grid on; title(椭圆齿轮节曲线与小圆圆心轨迹); legend(椭圆节曲线, 小圆圆心轨迹); subplot(2,2,2); plot(Px, Py, b-, LineWidth, 1.5); axis equal; grid on; title(异形齿轮万花尺轨迹); subplot(2,2,[3,4]); plot(x_e, y_e, k--, LineWidth, 0.5); hold on; plot(Cx, Cy, g:, LineWidth, 0.5); plot(Px, Py, b-, LineWidth, 2); axis equal; grid on; title(综合视图); legend(椭圆齿轮, 小圆圆心路径, 笔尖轨迹);这段代码实现了椭圆外齿轮下的万花尺仿真。其中最关键也是最容易出错的部分是坐标变换如何将笔尖相对于小圆圆心的位置正确地转换到世界坐标系。这涉及到对瞬时滚动接触点几何关系的深刻理解。上面的代码提供了一种基于瞬时切向角alpha的转换方法在大多数情况下是有效的但对于曲率变化剧烈的非圆齿轮可能需要更严谨的微分几何推导。4.3 仿真中的常见问题与调试技巧轨迹不闭合或错位这通常是因为数值积分误差累积或者对滚动条件dθ ds / r的理解有误。确保t参数从0到2π完整遍历椭圆一周并且最终计算出的theta(end)应该是2π的整数倍因为小圆也滚动回了初始方位。如果不是检查ds的计算和积分过程。笔尖轨迹出现预期外的扭结或跳跃这几乎总是坐标变换错误导致的。重点检查计算笔尖坐标的那一行代码。一个有效的调试方法是将d设为0此时笔尖轨迹应与小圆圆心轨迹(Cx, Cy)完全重合。如果不重合说明你的坐标变换逻辑有根本错误。计算速度慢对于复杂的非圆齿轮曲线如高阶傅里叶级数描述的曲线计算曲率半径和法向量可能很耗时。可以预先计算好所有几何量并向量化避免在循环中进行复杂函数调用。实操心得在仿真异形齿轮时我强烈建议从一个非常简单的非圆齿轮开始比如一个偏心圆圆心不在原点的圆。偏心圆本质上还是一个圆其曲率半径恒定因此其万花尺轨迹应该与标准圆内旋轮线类似只是中心偏移了。用这个简单案例验证你的仿真代码可以快速定位出是运动学关系错误还是坐标变换错误。5. 工程应用延伸与创意拓展完成了基础圆形和简单异形齿轮的仿真我们可以思考一下这个项目的更多可能性。5.1 在机械设计中的应用万花尺仿真的核心是共轭齿廓的接触运动仿真。这在非圆齿轮设计中非常重要。非圆齿轮用于实现变传动比常见于流量计、印刷机械和某些自动化机构中。通过本项目的仿真框架你可以验证齿廓设计给定一个非圆齿轮的节曲线仿真其与另一个齿轮圆或非圆的啮合运动检查是否存在干涉、传动是否平滑。生成复杂凸轮轮廓将笔尖视为从动件的端点小圆的运动规律作为输入生成的轨迹就可以作为凸轮的轮廓线。这为设计特定运动规律的凸轮提供了一种图形化方法。研究运动合成理解一个复杂平面轨迹如何分解为两个简单运动的合成公转自转这对于机构学分析很有帮助。5.2 在计算机图形与艺术中的应用万花尺图案本身具有强烈的数学美感可以很容易地拓展为创意生成工具。多色与多层绘制在仿真中可以设置多个不同d值的“笔尖”同时绘图或者让d值随时间周期性变化产生色彩斑斓、层次丰富的图案。三维空间拓展将原理扩展到三维让小圆在一个球面或其它曲面内滚动笔尖在三维空间中运动可以生成极其复杂的空间曲线用于三维装饰或雕塑设计。动画生成将绘制过程录制成视频本身就是一段优美的数学可视化材料。可以结合音乐生成视听作品。5.3 性能优化与交互式工具开发基础的MATLAB脚本已经足够好用但我们可以做得更专业。GUI开发利用MATLAB的App Designer或GUIDE创建一个图形用户界面让用户可以通过滑块实时调整R,r,d等参数并立即看到图案变化。这对于教学和艺术创作非常友好。参数化扫描与优化编写脚本自动搜索产生“美观”图案的参数组合。可以定义一些美学指标如图案的对称性、复杂度、自相交次数等然后使用优化算法如遗传算法在参数空间中寻找最优解。生成高分辨率矢量图MATLAB默认输出位图对于需要印刷或激光雕刻的图案可以优化代码直接生成SVG或EPS格式的矢量图形确保无限放大不失真。这个项目从一个小小的玩具出发贯穿了数学建模、数值计算、图形可视化和工程思维。它完美地展示了如何用计算思维去解构和重塑一个物理世界中的现象。无论你是想重温童年的乐趣还是寻找一个有趣的编程练习或是需要为某个工程问题寻找灵感万花尺的MATLAB仿真都是一个绝佳的起点。我最深的体会是在调试异形齿轮仿真的那一周里对微分几何和坐标变换的理解比过去几年都要深刻——有时候最好的学习就是动手去实现一个具体而微的项目。
1. 从玩具到数学万花尺的工程化仿真之旅小时候很多人都玩过一种叫“万花尺”的绘图玩具一个大尺环套着一个小齿轮齿轮上有孔用笔插进孔里推动小齿轮沿着大尺环的内壁滚动就能画出各种繁复而美丽的曲线。这种看似简单的玩具背后隐藏的数学原理却相当迷人——它就是圆内旋轮线。作为一名长期与MATLAB打交道的工程师我最近把儿时的玩具搬进了数字世界用MATLAB完整仿真了万花尺的绘制过程并深入到了“异形齿轮”的拓展领域。这不仅仅是一次怀旧更是一次将几何、运动学和数值计算融会贯通的绝佳实践。对于机械设计、动画生成乃至艺术创作领域的朋友来说理解如何用代码“驱动”一支虚拟的笔画出确定性的复杂图案具有很高的参考价值。本文将手把手带你复现这个过程并探讨如何将其原理应用于非圆齿轮的仿真让你不仅得到漂亮的图更能透彻理解每一行代码背后的“为什么”。2. 核心原理拆解圆内旋轮线的参数方程推导万花尺画出的曲线在数学上称为“圆内旋轮线”或“内摆线”。它的生成可以抽象为一个经典的几何运动模型一个半径为r的小圆在一个半径为R的大圆内部无滑动地滚动。我们关注的是固定在小圆上、距离小圆圆心为d的一个点P的运动轨迹。这个点P就对应万花尺上笔尖的位置。2.1 建立数学模型与坐标系为了用代码描述这个运动我们首先需要建立参数方程。设大圆圆心为原点O(0,0)。假设小圆初始时刻位于大圆的最右侧接触点小圆圆心C位于(R-r, 0)笔尖点P位于(R-rd, 0)。当小圆在大圆内逆时针滚动时它滚过的弧长等于大圆内壁被滚过的弧长无滑动条件。设小圆相对于自身圆心转过的角度为θ这是我们的核心运动参数那么小圆圆心C绕大圆圆心O转过的角度φ满足关系r * θ R * φ。由此可得φ (r / R) * θ。2.2 推导笔尖点P的坐标现在我们需要在固定的世界坐标系中求出点P的坐标。这个过程可以分解为两个变换的叠加小圆圆心C的位置C点绕原点O作圆周运动其半径为(R - r)当前角度为φ。因此C的坐标为Cx (R - r) * cos(φ)Cy (R - r) * sin(φ)笔尖P相对于圆心C的位置在小圆自身的坐标系中笔尖P位于距离圆心C为d的位置。由于小圆自身旋转了θ角度并且是在大圆内“滚动”而非“自转”其相对于大圆切向的方向需要仔细考虑。从C点指向P点的向量其方向角并不是简单的θ。因为小圆在滚动时其与地面的接触点速度为零这意味着小圆上每一点的运动是滚动公转和滑动自转的合成。正确的处理方式是P点相对于圆心C的位置是由小圆自身的旋转θ决定的但这个旋转是发生在随着圆心C移动的动坐标系中。并且由于滚动方向小圆在大圆内逆时针滚动与自转方向相反从接触点看P点在小圆坐标系中的角度应为-θ此处符号取决于坐标系定义是推导的关键。综合以上点P的坐标最终参数方程为设θ为自变量φ (r / R) * θ; Px (R - r) * cos(φ) d * cos(φ - θ); Py (R - r) * sin(φ) d * sin(φ - θ);其中(φ - θ)项描述了笔尖相对于动坐标系随小圆圆心平动并旋转的坐标系的方向。这个方程是仿真的基石。注意方程中的符号φ - θ还是θ - φ直接影响了曲线的旋向和形态。我强烈建议在首次实现时通过绘制小圆圆心C和笔尖P的瞬时位置向量来验证几何关系的正确性。一个快速的验证方法是当d 0时P点应与C点重合轨迹是一个圆当d r时P点应位于小圆周上应能画出标准的内摆线。3. MATLAB仿真实现从方程到动态绘图理解了原理用MATLAB实现就变得清晰了。我们的目标是1画出最终的静态曲线2动态展示绘制过程重现万花尺的运作。3.1 基础静态曲线绘制首先我们定义参数并计算轨迹点。这里的关键是参数θ的采样。为了让曲线闭合且光滑小圆需要滚动足够多的圈数直到笔尖回到起点。这要求R/r是一个有理数这样小圆滚动若干圈后公转角度φ将是2π的整数倍。我们通常取R和r为整数比。% 万花尺圆内旋轮线基础仿真 clear; clc; close all; % 1. 定义参数 R 100; % 大圆半径 r 31; % 小圆半径 d 65; % 笔尖到小圆圆心的距离 n_cycles 10; % 小圆滚动圈数根据R/r比值决定实际圈数 % 2. 计算采样点 % 小圆自转的总角度 圈数 * 2π theta_total n_cycles * 2 * pi; % 生成theta数组点数足够多以保证曲线光滑 num_points 5000; theta linspace(0, theta_total, num_points); % 3. 应用参数方程计算轨迹 phi (r / R) * theta; % 小圆圆心公转角度 Px (R - r) * cos(phi) d * cos(phi - theta); Py (R - r) * sin(phi) d * sin(phi - theta); % 4. 绘制静态结果 figure(Position, [100, 100, 800, 800]); plot(Px, Py, b-, LineWidth, 1.5); axis equal; grid on; box on; title(sprintf(圆内旋轮线 (R%d, r%d, d%d), R, r, d)); xlabel(X); ylabel(Y);运行这段代码你就能立即得到一幅精美的旋轮线图案。调整R,r,d三个参数你会发现图案千变万化当d r时笔尖在小圆内部曲线相对平滑。当d r时笔尖恰在小圆圆周上画出经典的内摆线可能有尖点。当d r时笔尖在小圆之外曲线会形成 loops环图案更为复杂。R与r的比值决定了图案的“瓣数”或对称性。若比值为整数n则曲线有n个尖点或n重对称性若为分数图案可能不闭合或周期极长。3.2 加入动态仿真与绘图过程静态图很美但动态过程更能揭示其生成机理。我们可以利用MATLAB的循环和drawnow命令来制作动画。% 动态绘制万花尺过程 figure(Position, [100, 100, 800, 800]); hold on; axis equal; grid on; box on; xlim([-R*1.2, R*1.2]); ylim([-R*1.2, R*1.2]); title(万花尺绘制动态仿真); xlabel(X); ylabel(Y); % 预先计算所有轨迹点同上 % ... % 初始化绘图对象 h_trace plot(Px(1), Py(1), b-, LineWidth, 1.5); % 轨迹线 h_point plot(Px(1), Py(1), ro, MarkerSize, 8, MarkerFaceColor, r); % 当前笔尖位置 h_circle_big rectangle(Position, [-R, -R, 2*R, 2*R], ... Curvature, [1, 1], EdgeColor, k, LineWidth, 1, LineStyle, --); h_circle_small plot(0, 0, g--, LineWidth, 1); % 小圆轮廓稍后更新 % 动态绘制循环 for i 1:50:num_points % 跳步绘制以提高动画速度 % 更新轨迹线绘制到当前点 set(h_trace, XData, Px(1:i), YData, Py(1:i)); % 更新笔尖位置 set(h_point, XData, Px(i), YData, Py(i)); % 计算并更新小圆轮廓 current_phi phi(i); current_theta theta(i); Cx (R - r) * cos(current_phi); Cy (R - r) * sin(current_phi); % 生成小圆上的点 alpha linspace(0, 2*pi, 100); small_circle_x Cx r * cos(alpha); small_circle_y Cy r * sin(alpha); set(h_circle_small, XData, small_circle_x, YData, small_circle_y); % 刷新图形 drawnow; % 可加入 pause(0.01) 控制速度 end hold off;这个动态仿真就像在电脑里运行了一个虚拟的万花尺。你可以清晰地看到小圆如何在大圆内滚动笔尖如何一步步勾勒出复杂图案。这对于教学和理解运动合成非常有帮助。3.3 参数探索与图案库生成万花尺的魅力在于参数微调就能产生截然不同的图案。我们可以写一个简单的脚本批量生成不同参数下的曲线建立一个“数字万花尺图案库”。% 参数探索固定R和d变化r R 100; d 60; r_values [19, 29, 33, 41, 47, 51]; % 尝试不同的r值 figure(Position, [50, 50, 1200, 600]); for idx 1:length(r_values) r r_values(idx); % 计算轨迹确保闭合 n lcm(R, r) / r; % 计算最小公倍数除以r得到所需圈数近似 theta linspace(0, n*2*pi, 3000); phi (r / R) * theta; Px (R - r) * cos(phi) d * cos(phi - theta); Py (R - r) * sin(phi) d * sin(phi - theta); subplot(2, 3, idx); plot(Px, Py, LineWidth, 1); axis equal; axis off; title(sprintf(r %d, r)); end sgtitle(不同小圆半径r下的万花尺图案 (R100, d60));通过这样的批量实验你可以直观地感受到数学的魔力并可能发现一些特别漂亮的参数组合。这本身就是一种创造性的探索。4. 进阶挑战异形齿轮非圆齿轮的万花尺仿真传统的万花尺基于两个圆形齿轮内啮合。那么如果大尺环或小齿轮不是标准的圆形呢这就是“异形齿轮”的范畴。仿真非圆齿轮的万花尺核心挑战在于无滑动滚动条件变得复杂因为齿轮的节曲线相当于等效滚动圆半径不再是常数。4.1 非圆齿轮节曲线与运动关系假设我们有一个椭圆形的外齿轮大尺环和一个圆形的小齿轮。椭圆没有恒定的曲率半径因此小圆圆心到接触点的距离即瞬时滚动半径是变化的。设椭圆参数方程为x_e a * cos(t), y_e b * sin(t)其中a b。 椭圆上某一点t的曲率半径ρ(t)公式为ρ(t) ( (a^2 * sin^2(t) b^2 * cos^2(t))^(3/2) ) / (a * b)在无滑动滚动条件下小圆半径r的圆心C不再沿一个固定半径的圆运动而是沿着椭圆的一条等距线运动该等距线与椭圆的距离为(R_eff - r)但这里的R_eff是时变的等于椭圆的曲率半径ρ(t)。这使得运动关系从一个简单的线性比例φ (r/R)*θ变成了一个需要数值积分求解的微分关系。具体来说设椭圆上的接触点参数为t小圆自转角度为θ。无滑动条件意味着小圆滚过的弧长等于椭圆上被滚过的弧长r * dθ ds其中ds是椭圆的弧长微分ds sqrt( (dx_e/dt)^2 (dy_e/dt)^2 ) * dt。 因此有关系式dθ (1/r) * ds。 同时小圆圆心C的位置是椭圆上该接触点沿内法线方向偏移(ρ(t) - r)的距离假设小圆在椭圆内部滚动。4.2 MATLAB数值仿真实现由于解析解很难获得我们采用数值方法逐步积分求解。% 异形齿轮椭圆外齿轮 圆形小齿轮万花尺仿真 clear; clc; close all; % 1. 定义椭圆参数和小圆参数 a 80; % 椭圆长半轴 b 50; % 椭圆短半轴 r 15; % 小圆半径 d 25; % 笔尖距离小圆圆心距离 % 2. 定义椭圆上的参数t数组对应接触点 num_steps 5000; t linspace(0, 2*pi, num_steps); % 椭圆参数角从0到2π % 3. 预计算椭圆几何量 x_e a * cos(t); y_e b * sin(t); dx_dt -a * sin(t); % x对t的导数 dy_dt b * cos(t); % y对t的导数 ds_dt sqrt(dx_dt.^2 dy_dt.^2); % 弧长微分 ds/dt % 计算曲率半径 ρ(t) numerator (a^2 * sin(t).^2 b^2 * cos(t).^2).^(3/2); denominator a * b; rho_t numerator ./ denominator; % 曲率半径 % 4. 数值积分求解小圆自转角度 theta theta zeros(size(t)); for i 2:num_steps % 弧长微元 ds ds ds_dt(i) * (t(i) - t(i-1)); % 根据无滑动条件r * dθ ds d_theta ds / r; theta(i) theta(i-1) d_theta; end % 5. 计算小圆圆心C的轨迹 % 椭圆上点(x_e, y_e)处的单位内法向量 N (-dy_dt, dx_dt) / norm % 注意方向要保证是“内”法线对于参数化方程 (a*cos(t), b*sin(t))内法线方向需根据实际情况判断。 % 一个简单方法是对于标准椭圆从原点指向椭圆上点的向量是 (x_e, y_e)其内法线大致指向椭圆中心。 % 更稳健的方法是计算单位切向量 T (dx_dt, dy_dt)/ds_dt然后通过旋转90度得到法向量。 % 这里我们计算单位切向量然后顺时针旋转90度得到内法向量假设逆时针参数化。 T [dx_dt./ds_dt; dy_dt./ds_dt]; % 2 x N 矩阵 % 顺时针旋转90度: (x, y) - (y, -x) N [T(2, :); -T(1, :)]; % 这给出了近似内法向量 % 圆心C 椭圆上的点 (ρ(t) - r) * 内法向量 Cx x_e (rho_t - r) .* N(1, :); Cy y_e (rho_t - r) .* N(2, :); % 6. 计算笔尖P的轨迹 % 笔尖相对于圆心C的向量长度为d方向角需要考虑小圆自身的旋转theta % 在动坐标系随C点平动且随小圆旋转中笔尖的方向角为 -theta滚动方向导致 % 需要将这个向量转换到世界坐标系。动坐标系的x轴方向是椭圆切向方向Ty轴方向是法向N。 % 因此笔尖在世界坐标系的偏移量为 d * [cos(-theta); sin(-theta)] 在 (T, N) 基下的坐标 % 即 offset d * (cos(theta).*T (-sin(theta)).*N) 这里需要仔细推导。 % 更直接的方法笔尖相对于C的位置在小圆自身的固定坐标系中角度为 -theta。 % 但小圆自身的固定坐标系相对于世界坐标系的旋转角度是多少这等于接触点切线与世界坐标系x轴的夹角。 % 设切向角为 alpha atan2(dy_dt, dx_dt) alpha atan2(dy_dt, dx_dt); % 那么笔尖P的坐标为 Px Cx d * cos(alpha - theta); Py Cy d * sin(alpha - theta); % 7. 绘图 figure(Position, [100, 100, 900, 900]); subplot(2,2,1); plot(x_e, y_e, k--, LineWidth, 1); hold on; plot(Cx, Cy, g-, LineWidth, 1); axis equal; grid on; title(椭圆齿轮节曲线与小圆圆心轨迹); legend(椭圆节曲线, 小圆圆心轨迹); subplot(2,2,2); plot(Px, Py, b-, LineWidth, 1.5); axis equal; grid on; title(异形齿轮万花尺轨迹); subplot(2,2,[3,4]); plot(x_e, y_e, k--, LineWidth, 0.5); hold on; plot(Cx, Cy, g:, LineWidth, 0.5); plot(Px, Py, b-, LineWidth, 2); axis equal; grid on; title(综合视图); legend(椭圆齿轮, 小圆圆心路径, 笔尖轨迹);这段代码实现了椭圆外齿轮下的万花尺仿真。其中最关键也是最容易出错的部分是坐标变换如何将笔尖相对于小圆圆心的位置正确地转换到世界坐标系。这涉及到对瞬时滚动接触点几何关系的深刻理解。上面的代码提供了一种基于瞬时切向角alpha的转换方法在大多数情况下是有效的但对于曲率变化剧烈的非圆齿轮可能需要更严谨的微分几何推导。4.3 仿真中的常见问题与调试技巧轨迹不闭合或错位这通常是因为数值积分误差累积或者对滚动条件dθ ds / r的理解有误。确保t参数从0到2π完整遍历椭圆一周并且最终计算出的theta(end)应该是2π的整数倍因为小圆也滚动回了初始方位。如果不是检查ds的计算和积分过程。笔尖轨迹出现预期外的扭结或跳跃这几乎总是坐标变换错误导致的。重点检查计算笔尖坐标的那一行代码。一个有效的调试方法是将d设为0此时笔尖轨迹应与小圆圆心轨迹(Cx, Cy)完全重合。如果不重合说明你的坐标变换逻辑有根本错误。计算速度慢对于复杂的非圆齿轮曲线如高阶傅里叶级数描述的曲线计算曲率半径和法向量可能很耗时。可以预先计算好所有几何量并向量化避免在循环中进行复杂函数调用。实操心得在仿真异形齿轮时我强烈建议从一个非常简单的非圆齿轮开始比如一个偏心圆圆心不在原点的圆。偏心圆本质上还是一个圆其曲率半径恒定因此其万花尺轨迹应该与标准圆内旋轮线类似只是中心偏移了。用这个简单案例验证你的仿真代码可以快速定位出是运动学关系错误还是坐标变换错误。5. 工程应用延伸与创意拓展完成了基础圆形和简单异形齿轮的仿真我们可以思考一下这个项目的更多可能性。5.1 在机械设计中的应用万花尺仿真的核心是共轭齿廓的接触运动仿真。这在非圆齿轮设计中非常重要。非圆齿轮用于实现变传动比常见于流量计、印刷机械和某些自动化机构中。通过本项目的仿真框架你可以验证齿廓设计给定一个非圆齿轮的节曲线仿真其与另一个齿轮圆或非圆的啮合运动检查是否存在干涉、传动是否平滑。生成复杂凸轮轮廓将笔尖视为从动件的端点小圆的运动规律作为输入生成的轨迹就可以作为凸轮的轮廓线。这为设计特定运动规律的凸轮提供了一种图形化方法。研究运动合成理解一个复杂平面轨迹如何分解为两个简单运动的合成公转自转这对于机构学分析很有帮助。5.2 在计算机图形与艺术中的应用万花尺图案本身具有强烈的数学美感可以很容易地拓展为创意生成工具。多色与多层绘制在仿真中可以设置多个不同d值的“笔尖”同时绘图或者让d值随时间周期性变化产生色彩斑斓、层次丰富的图案。三维空间拓展将原理扩展到三维让小圆在一个球面或其它曲面内滚动笔尖在三维空间中运动可以生成极其复杂的空间曲线用于三维装饰或雕塑设计。动画生成将绘制过程录制成视频本身就是一段优美的数学可视化材料。可以结合音乐生成视听作品。5.3 性能优化与交互式工具开发基础的MATLAB脚本已经足够好用但我们可以做得更专业。GUI开发利用MATLAB的App Designer或GUIDE创建一个图形用户界面让用户可以通过滑块实时调整R,r,d等参数并立即看到图案变化。这对于教学和艺术创作非常友好。参数化扫描与优化编写脚本自动搜索产生“美观”图案的参数组合。可以定义一些美学指标如图案的对称性、复杂度、自相交次数等然后使用优化算法如遗传算法在参数空间中寻找最优解。生成高分辨率矢量图MATLAB默认输出位图对于需要印刷或激光雕刻的图案可以优化代码直接生成SVG或EPS格式的矢量图形确保无限放大不失真。这个项目从一个小小的玩具出发贯穿了数学建模、数值计算、图形可视化和工程思维。它完美地展示了如何用计算思维去解构和重塑一个物理世界中的现象。无论你是想重温童年的乐趣还是寻找一个有趣的编程练习或是需要为某个工程问题寻找灵感万花尺的MATLAB仿真都是一个绝佳的起点。我最深的体会是在调试异形齿轮仿真的那一周里对微分几何和坐标变换的理解比过去几年都要深刻——有时候最好的学习就是动手去实现一个具体而微的项目。
