1. 为什么想到做“旋转圆盘上的摆动模拟”事情源于某次科普演示的准备。当时要给一群同学讲“非惯性参考系里的运动”PPT 上公式推了一遍又一遍大家点头都点得很勉强。说实话光看科里奥利力的表达式谁都觉得抽象——一个虚拟力既没有施力物体方向又垂直于速度光靠文字描述很难建立直觉。于是我想做一个小工具把一个摆放在旋转圆盘上用数值模拟把摆锤的轨迹画出来。转盘转着摆锤一边周期摆动一边被“看不见的力”牵着走轨迹会形成各种漂亮的图案。这种“把虚拟力变成可见轨迹”的效果比任何公式都直观。项目不大但涉及的东西很典型旋转参考系建模、运动方程推导、数值积分、双坐标系可视化几乎把物理模拟入门的关键环节都踩了一遍。这个模拟主要解决两类问题。第一类是教学内容可视化比如用轨迹展示科里奥利力的作用方向解释傅科摆的摆平面为什么会进动第二类是建模预习在做更复杂的旋转机械、飞行器姿态仿真之前先用这个简单模型理解旋转坐标系里的动力学方程长什么样。适合的人群也很明确物理或力学方向的课程设计、对数值模拟感兴趣的 Python 初学者、以及所有想把“悬在空中的公式”落到屏幕上的人。我最终实现的版本是用 Python 写的数学部分手动推导积分用 RK4绘图用 Matplotlib效果分两个视角同时显示一个视角固定在转盘上跟摆相对静止地观察一个视角站在地面上俯视全局。这个双视角设计在后来的调试中帮了大忙后面会详细说。2. 运动方程是怎么推出来的惯性系、科里奥利力与小角度近似2.1 先看惯性系里的摆要模拟一个物理过程第一步是写运动方程。摆的模型我选了最常见的“理想单摆”摆长固定为 L摆锤质量集中一点悬挂点固定在转盘正上方的竖直轴上绳始终可以自由转动摆锤在水平面附近活动。站在地面惯性系看摆锤受到的力只有两个重力 m g 和绳张力 T。重力竖直向下张力沿绳方向指向悬挂点。如果摆角很小比如不超过 10°水平方向的合力近似与水平位移成正比、方向相反也就是经典的简谐振动a -(g / L) * r这里 r 是摆锤相对悬挂点竖直投影点的水平位移矢量。定义固有频率 ω_0 sqrt(g/L)那么惯性系中水平运动就是一个以 ω_0 为角频率的简谐振动。在没有任何阻尼和驱动力的情况下摆锤应该在同一个竖直平面里来回摆动轨迹是一条直线段。这个结论本身很无聊但它是后面一切的基础。2.2 换到旋转参考系虚拟力登场现在难点来了。我们要模拟的是“站在转盘上看到的摆动”也就是以转盘为参考系去描述摆锤的运动。这个参考系和地面之间有一个恒定的旋转角速度 Ω它不是惯性系所以牛顿第二定律不能直接硬套。这里需要用到加速度合成公式。假设有一个相对转盘静止的点它在地面系里实际在做圆周运动有向心加速度如果物体相对转盘还有速度还会额外产生科里奥利加速度。把惯性系下的真实加速度和旋转参考系下的加速度之间的关系写出来就是a_惯性 a_相对 2 Ω × v_相对 Ω × (Ω × r)公式里第一项是“相对转盘的加速度”第二项来自科里奥利效应第三项是向心加速度项。在匀速旋转、悬挂点就在转轴上的情况下位置矢量 r 在两个参考系里是同一个量所以可以把惯性系中的真实加速度等于摆的回复加速度代入a_相对 2 Ω × v_相对 Ω × (Ω × r) -ω_0² r移项之后就得到旋转参考系下的运动方程a_相对 -ω_0² r - 2 Ω × v_相对 - Ω × (Ω × r)这样写可能还不够直观我把两个矢量叉乘展开。设 Ω 的方向沿竖直轴也就是垂直纸面转盘坐标系里摆锤的水平坐标为 (x, y)速度为 (vx, vy)则科里奥利加速度-2 Ω × v (2 Ω vy, -2 Ω vx)离心加速度-Ω × (Ω × r) (Ω² x, Ω² y)于是最终得到一组二维的常微分方程组x -ω_0² x 2 Ω vy Ω² xy -ω_0² y - 2 Ω vx Ω² y2.3 生活化的解释如果你觉得公式抽象可以用一个生活场景帮自己理解。想象你在一个正在转弯的公交车里站着车转弯时你感觉自己被“甩”向车厢外侧这就是离心项如果你在车厢里朝车头方向走会感觉自己被一股力量推向侧边那就是科里奥利项。旋转圆盘上的摆也是这样站在转盘上的人看到的“异常偏转”背后就是这两个虚拟力在起作用。这一组方程中x 和 y 是相互耦合的x 方向的运动会通过科里奥利项影响 y 方向反之亦然。这种耦合正是轨迹会出现旋转、椭圆、花瓣图案的原因。如果 Ω 0方程退化为两个独立的简谐振动轨迹就是一条直线一旦 Ω 不等于 0摆锤的运动立刻变得“立体”起来。需要注意这组方程成立的前提是小角度近似。如果摆角超过约 15°sin θ 与 θ 的差别就不能忽略了水平方向回复力也不再严格等于 -ω_0² x。实际模拟时我会把初始摆角控制在 10° 以内这样线性模型足够精确。如果你想要更大角度的真实模拟就必须把完整的非线性项写进去后面第 5 节会讲这个坑。3. 数值积分与可视化完整 Python 实现3.1 无量纲化很关键直接代入实际物理量也能跑但效果不好调。摆长 0.25 m 时 ω_0 约 6.26 rad/s转盘角速度如果设成 0.5 rad/s两个频率差了十几倍画出来的图案会很扁参数扫描时也容易算错数量级。我习惯先把方程无量纲化令 ω_0 1时间单位取 1/ω_0长度单位取初始摆幅 A这样只剩下一个自由参数 Ω/ω_0。无量纲化后的方程是x -x 2 ω vy ω² xy -y - 2 ω vx ω² y其中 ω Ω / ω_0。程序里我固定 ω_0 1只改 ω就能观察到所有典型的运动形态。这个习惯非常重要强烈建议你也这么做——它会让参数扫描变得异常轻松。3.2 RK4 求解器实现ODE 积分的方案我选了经典的龙格-库塔法也就是 RK4。为什么不用更简单的欧拉法因为摆动运动是周期性的欧拉法每步都会引入相位误差和能量漂移几百步之后轨迹会明显“发胖”看起来像加了阻尼一样而 RK4 是四阶精度误差随步长 h 的 4 次方衰减对它来说 h 0.01 就能得到非常干净的轨迹。核心代码很简洁。把状态向量设为 [x, y, vx, vy]函数 f(state) 返回导数import numpy as np def derivatives(state, omega): x, y, vx, vy state ax -x 2 * omega * vy omega**2 * x ay -y - 2 * omega * vx omega**2 * y return np.array([vx, vy, ax, ay]) def rk4_step(state, dt, omega): def f(s): return derivatives(s, omega) k1 f(state) k2 f(state 0.5 * dt * k1) k3 f(state 0.5 * dt * k2) k4 f(state dt * k3) return state (dt / 6.0) * (k1 2 * k2 2 * k3 k4)主循环就更好写了def simulate(omega, dt0.01, total_time80.0, x01.0): state np.array([x0, 0.0, 0.0, 0.0]) steps int(total_time / dt) traj np.zeros((steps, 2)) for i in range(steps): traj[i] state[:2] state rk4_step(state, dt, omega) return traj我用它跑了 ω 0.3、初始位移 x 1、初速度 0 的情况总时长 80 个时间单位相当于 80/2π ≈ 12.7 个摆周期。轨迹是一圈一圈逆时针进动的椭圆非常漂亮。3.3 双视角动画转盘坐标系 vs 地面坐标系数值算出来的 (x, y) 是相对转盘的位置。但人站在地面上看摆锤的真实轨迹还需要把坐标旋转回去。转盘转了 Ω t 角度所以惯性系坐标是X x cos(ω t) - y sin(ω t)Y x sin(ω t) y cos(ω t)我用 Matplotlib 的 FuncAnimation 画了两个并排的子图左图显示转盘参考系中的轨迹和摆锤位置右图显示地面参考系的轨迹。左图能看到轨迹在转盘坐标系中缓慢进动右图则能看到接近直线的往复摆动小 ω 时。两个图对比起来非常震撼。import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation # 预积分整个时间序列 states np.zeros((steps, 4)) s np.array([1.0, 0.0, 0.0, 0.0]) for i in range(steps): states[i] s s rk4_step(s, dt, omega) # 计算惯性系坐标 theta omega * np.arange(steps) * dt X_inert states[:, 0] * np.cos(theta) - states[:, 1] * np.sin(theta) Y_inert states[:, 0] * np.sin(theta) states[:, 1] * np.cos(theta) fig, (ax1, ax2) plt.subplots(1, 2, figsize(11, 5)) ax1.set_xlim(-1.5, 1.5) ax1.set_ylim(-1.5, 1.5) ax1.set_aspect(equal) ax1.set_title(转盘参考系) line1, ax1.plot([], [], lw1) dot1, ax1.plot([], [], o, markersize8) ax2.set_xlim(-1.5, 1.5) ax2.set_ylim(-1.5, 1.5) ax2.set_aspect(equal) ax2.set_title(地面惯性系) line2, ax2.plot([], [], lw1) dot2, ax2.plot([], [], o, markersize8) def update(frame): line1.set_data(states[:frame, 0], states[:frame, 1]) dot1.set_data(states[frame, 0], states[frame, 1]) line2.set_data(X_inert[:frame], Y_inert[:frame]) dot2.set_data(X_inert[frame], Y_inert[frame]) return line1, dot1, line2, dot2 ani FuncAnimation(fig, update, framesrange(0, steps, 5), interval20, blitTrue) plt.tight_layout() plt.show()这个代码跑起来之后左图可以看到摆锤轨迹像一个旋转的椭圆右图则几乎是一根直线在来回扫动。有趣的是如果你只看左图会误以为摆锤在绕着一个无形的中心打转而切换到右图你才明白地面上的真实运动其实非常“规矩”。这就是参考系选择对观测者判断的影响也是这个模拟最有说服力的地方。跑了大概两三分钟后我得到了一组稳定漂亮的进动轨迹。把这个双视角动画保存下来直接扔进科普演示里效果比 PPT 上画十条力线都好。如果你只是想要静态图也可以把 update 逻辑去掉直接 plot 全轨迹。4. 摆轨迹的形态频率比 ω 才是真正的“导演”4.1 慢速旋转摆平面的进动参数 ω即 Ω/ω_0是决定轨迹形态的唯一无量纲数。我做的第一个实验是 ω 0.1也就是转盘转得非常慢摆的固有频率是转盘频率的 10 倍。初始位移沿 x 方向无初速度。在转盘参考系里轨迹看上去仍然近似是一条直线段但这条直线的方向会慢慢旋转——大约转完一整圈需要的时间是 2π/ω也就是一个转盘周期。这个现象就是傅科摆的摆平面进动在旋转参考系中摆的运动平面不再固定而是以某个角速度缓慢旋转。公式上解释很简单当 ω 很小时离心项 ω² x 小到可以忽略方程近似为两个相同频率的振动通过科里奥利项耦合耦合的结果就是本征模式变成沿两个正交方向的圆偏振轨迹按 ω 的频率进动。在地面惯性系里看摆动平面反而基本上是固定的转盘每转一圈摆平面几乎回到原方向。左图右图一对比进动的真正来源就一目了然——它不是某种神秘的“额外力”而是观测者本身在旋转。4.2 中等转速椭圆进动与花瓣轨迹把 ω 调大一些到 0.5情况开始变得复杂。转盘参考系里轨迹不再是直线段而是一个不断旋转并变形的椭圆如果给的初始条件不是纯径向位移而是带一点切向初速度轨迹就会变成漂亮的多瓣玫瑰线。比如初始状态设为 x0 1, y0 0, vx 0, vy 0.3给一个沿 y 方向的初始速度ω 0.5积分 60 个时间单位后轨迹是一个逆时针进动的椭圆长轴也在缓慢旋转。而如果初始速度更大一些vy 0.8轨迹会出现明显的交叉花瓣。交叉点的位置、每片花瓣的夹角都由 ω 决定网上那些“利萨茹花样摆”的动图本质就是这一类运动。这个阶段我做了几组参数扫描结果整理成了一张表方便以后回看ω Ω/ω_0初始条件转盘参考系中的典型轨迹地面系中的轨迹0.1x01v0近似直线段缓慢进动固定方向的直线摆动0.3x01v0椭圆定向旋转轻微弯曲的直线0.5x01v0明显进动的椭圆轨迹出现弧形偏转0.5x01vy0.8多瓣交叉轨迹曲线缠绕0.9x01v0轨迹剧烈缠绕近似圆形包络摆锤做复杂的曲线运动看到 ω 0.9 这行轨迹的包络几乎变成了一个圆摆周期和转盘周期接近科里奥利力与回复力“共振”运动变得很不稳定。而一旦 ω 超过 1离心项 ω² x 的等效“反弹簧”效应压过了回复力线性模型下轨迹会指数发散。这是模拟的一个边界需要在下一节专门讲一下。4.3 离心项的本质等效刚度在变软这组方程里最微妙的是离心项。把 x 方向的项合并一下x -(1 - ω²) x 2 ω vy也就是说旋转让摆的恢复力从 1 变成了 1 - ω²。ω 越大等效刚度越小摆的“表现频率”越低。当 ω 接近 1 时表现频率趋近于 0意味着在转盘上观察摆几乎不再摆动而是慢悠悠地漂移超过 1 之后平衡点翻转摆会被“离心力”推向无穷远。真实物理中这对应的是圆锥摆摆角变大线性模型已经失真但数值模拟里你依然能看到这个临界点——它像一面警示牌提醒你模型的使用边界。我在代码里加了一个简单检查如果轨迹半径超过初始半径的 20 倍就主动终止模拟并提示“omega 超出线性模型适用范围”。这个小功能帮我在演示时避免了不少尴尬场面。5. 模拟中的坑与排查经验5.1 数值发散第一排查项时间步长我第一次跑这个模拟设 dt0.1总时长 50结果到 20 个时间单位之后轨迹就“炸”了摆锤飞到屏幕外轨迹乱成一团。第一反应以为是物理模型写错了查了一遍方程没问题最后才反应过来是 RK4 的步长太大。摆动系统的周期是 2π ≈ 6.28我设 dt0.1相当于一个周期只采样 63 步这对 RK4 来说其实已经够用但问题是科里奥利耦合项的频率可能更高尤其在 ω 接近 1 时系统表现得像高频振荡必须有更细的步长。经验法则对于 RK4一个周期至少采样 100 步才谈得上轨迹平滑我最终固定 dt0.01也就是约 628 步一个周期轨迹非常稳定。调试时我建议先画能量的时间曲线。理想单摆机械能守恒如果能量曲线抬升或下降就是数值误差过大的信号。5.2 “非线性坑”小角度假设失效时怎么办前文一直在强调小角度近似。我之前在演示环节手一抖把初始位移 x0 设成了 1.2相对于摆长来说已经超过 15°结果转盘参考系的轨迹和线性模型预期的花瓣图案差别很大还出现了不规则的漂移。原因在于大摆角下水平回复力不再严格与位移成正比还要考虑竖直方向的位置变化相当于刚性摆的自由振荡频率会随振幅变化。这就引入了非线性频率调制轨迹形态自然偏离“标准答案”。如果你需要模拟大摆角建议把运动方程改成完整的二维球面摆形式或者直接引入坐标变换和拉格朗日力学。我后来在扩展版里加了 sin θ 形式的非线性项但只在 ω 较小时使用否则数值刚性会很强。另一个相关的坑是阻尼。真实的摆是有空气阻尼的阻尼项会让轨迹从旋转的椭圆慢慢收缩成一条线最终停在平衡位置。如果你在动画里看到轨迹越来越“实心”先检查方程里是不是不小心带了阻尼系数而不是急着怀疑代码。5.3 动画帧与真实时间的混淆Matplotlib 的 FuncAnimation 里参数 interval 单位是毫秒是“每帧的显示间隔”不是物理时间 dt。初始我把 interval 设成 dt*1000以为这样能保证动画和真实时间同步结果转盘参考系里的进动慢得让人怀疑程序卡了。实际上动画的播放速度由 frame 之间的间隔决定而不是物理时间步长。我的做法是每 5 个物理积分步绘制一帧interval 设 20 ms这样视觉上大约 100 步/秒折算成物理时间是 1 个时间单位/秒既能看清轨迹的细节又不会拖沓。想让动画快一点就增大 frames 的步长想慢一点就减小步长——和 dt 没有直接关系。5.4 最容易忽略的参考系坐标 vs 视图坐标最后一个坑很隐蔽。我在绘制转盘参考系轨迹时直接 plot(states[:,0], states[:,1])这没错。但如果你在旁边加一个表示圆盘边缘的圆形线框要注意圆盘上的固定点比如某个标记在动画中也是旋转的如果不画出来观众很容易产生“转盘没转、只有摆在动”的错觉。我在右图里花了点功夫除了轨迹还画了一个跟随转盘转动的半径线。动画播放时这根线的旋转能立刻告诉观众转盘本身的运动状态。一个小小的装饰对理解整个物理现象起到了决定性的作用。后来我给同学们展示时他们几乎第一眼就能看出“左边是在转盘上看右边是在地上看”。6. 从玩具模型到真实系统扩展与改进方向6.1 加入阻尼项真实摆不可能永远摆动。加入线性阻尼后方程变成x -x - 2 ζ ω_0 x 2 ω vy ω² xy -y - 2 ζ ω_0 y - 2 ω vx ω² y其中 ζ 是阻尼比通常取 0.01~0.1。修改后的模拟轨迹会从椭圆逐渐收缩成一个点这在演示“能量耗散与旋转参考系的关系”时很有用。代码只需在 derivatives 函数里加一项即可五秒钟的事。6.2 悬挂点不在旋转轴上如果悬挂点不在圆盘中心而是固定在圆盘边缘的一根立柱上位置矢量就不再是单纯的 r方程里会出现额外的圆周运动项。物理图像会完全改变摆锤相对地面的运动是“转盘公转 自身摆动”的叠加轨迹更像是旋转的摆轮。这个方向适合用来模拟旋转机械中的摆式阻尼器等工程对象但方程推导要重新来过不能直接套上一节的公式。6.3 从二维到三维平面模拟把运动限制在水平面实际上圆锥摆有竖直方向的运动。想更真实地模拟建议直接上三维刚体摆模型状态向量扩展到 6 维用四元数或欧拉角描述姿态。这个时候手写 RK4 开始吃力我会改用成熟的积分库比如 scipy 的 solve_ivp或者干脆用物理引擎它们内置了更稳健的自适应步长算法能处理刚性方程。不过我要提醒一句换成重工具之前先把二维线性模型吃透。我在这个项目里最大的收获不是写出来的代码而是彻底弄懂了科里奥利项和离心项的物理意义。6.4 一个实用小技巧最后分享一个小技巧。如果你想把参数 ω 连续变化观察轨迹形态的演化可以做一个滑杆交互控件。滑块从 0.05 拖到 0.9转盘参考系里的轨迹慢慢从直线变成椭圆再变成花瓣地面系里从直线变成弧线再变成乱麻整个过程非常直观。我在课程演示中就用这个滑杆让学生自己探索效果比任何静态图片都好。做物理模拟最忌讳的就是程序跑出漂亮的图却说不清每一笔轨迹背后的力。这个旋转圆盘摆动模拟的精髓恰恰在于小角度、线性化、RK4 这些看似朴素的选择每一步都让你更靠近物理的核心。试着把 ω 拉大盯着转盘参考系里的轨迹变化你会真正“看见”科里奥利力在跳舞。