“旋转圆盘上的摆动模拟”听起来像个物理课建模题但真正落地做的时候你会发现它兼顾了刚体动力学、非惯性系分析和数值稳定性这几件麻烦事。我最早是因为一个同学在某个平台团建展示了一个演示圆盘转起来摆锤没有老老实实顺着惯性走而是在盘面上画出漂亮的进动花纹。当时几个人都觉得该直接套旋转坐标系方程结果自己动手一算方向反了、轨迹炸了、能量凭空涨起来的坑一个接一个。后来我干脆把整个模拟做成了一套可复现的代码和实验流程也顺手把几个容易混淆的点理清楚了。这篇就把完整思路、推导过程和实现细节都摊开聊一遍新手照着跑能出图老手也能在参数边界和参考系符号上互相验证一下。1. 从物理直觉到建模目标到底在模拟什么1.1 一个最容易说错的直觉看到“旋转圆盘上的摆”大部分人的第一反应是“摆当然跟着圆盘转”。这个直觉在慢速旋转时是对的在快速旋转时会暴露问题。真正困扰人的地方在于摆球并没有被固定在盘面上它通过绳子悬挂在旋转轴上方的支点。圆盘转动只会带动支点转动而摆球本身处于非惯性系中必须把离心力和科里奥利力都算进去。我一开始也犯了同样的错假设初始让摆球相对于圆盘静止就等同于让它随盘共转。实际上摆球一旦相对盘面有速度科里奥利力就会立刻改变它的运动方向摆动平面不可能老老实实跟着圆盘转。模拟的目的就是把这种看起来反直觉的进动关系定量展现出来。1.2 这个模拟解决什么问题说直接一点这个模拟解决三件事。第一验证在旋转参考系下动力学方程写得对不对尤其是科里奥利项的符号和交叉耦合关系。第二通过改变圆盘角速度、摆长、初始振幅直观对比不同参数下摆锤轨迹的差异。第三它也是理解傅科摆原理的绝佳小型模型——傅科摆本质上是地球这个旋转参考系对长摆施加科里奥利力的表现而这里不过是把地球换成了一个可以随意调速的圆盘。适合看这篇文章的人包括正在学理论力学、想让抽象公式变成真实曲线的学生也包括在写游戏物理、机器人仿真时需要处理旋转坐标系下物体运动的开发者。代码量不大但涉及的概念密度很高。2. 建立运动方程非惯性系里的一笔一划2.1 坐标系、变量和三个力把圆盘角速度设为 Ω方向竖直向上也就是沿着 z 轴。在圆盘参考系里摆球相对盘心的位置用 (x, y) 表示摆长 L重力加速度 g。因为是小角度摆动绳子在水平面的投影远小于 L可以用水平方向近似回复力。在非惯性系里写牛顿第二定律必须加上两类惯性力。一类是离心力m Ω² (x, y)它总是把摆球往外推。另一类是科里奥利力-2m Ω × v′。这里 v′ 是摆球相对圆盘的速度叉乘展开后科里奥利力在 x 方向的大小是 2mΩ vy′在 y 方向则是 -2mΩ vx′。注意这里的交叉耦合x 方向的方程里出现 y 方向速度y 方向的方程里出现 x 方向速度这是“旋转坐标系会让运动方向偏转”的数学根源。2.2 小角度下的一阶方程摆球的回复力水平分量近似为 -mg x / L 和 -mg y / L。把惯性力和回复力相加约掉质量 m得到x″ -(g/L)x Ω² x 2Ω y′ y″ -(g/L)y Ω² y - 2Ω x′我见过很多人在这里把 Ω² 项当成小量直接删掉其实要看场景。如果圆盘转速超过摆的自然振荡频率也就是 Ω 接近 sqrt(g/L) 甚至更大离心项会让摆的有效回复力变成负的摆球会直接“飞”到圆盘边缘。所以稳妥做法是保留这一项只在解释解析解时讨论 Ω² 远小于 g/L 的极限。如果确实处于慢转极限比如 Ω 只有 0.1 rad/s 而 sqrt(g/L) 是 3.13 rad/s那么可以省略 Ω² 项得到更简洁的形式x″ -ω² x 2Ω y′ y″ -ω² y - 2Ω x′其中 ω sqrt(g/L) 是摆的自然角频率。解这个线性方程组时设复变量 z x iy方程会化成一阶复数形式。最终可以得到两个特征频率ω Ω 和 ω - Ω。这一对频率之差正是摆动平面在圆盘参考系中缓慢进动的原因。也就是说即便你站在圆盘上看摆锤画的是近似直线往复那条直线本身也在以约 Ω 的角速度反方向旋转。2.3 什么时候必须放弃小角度近似小角度方程只适合振幅不超过大约 10° 到 15° 的情况。振幅一旦变大绳子张力在水平方向的分量不再与水平位移严格成正比摆动频率也会随振幅漂移轨迹会出现明显的谐波畸变。如果圆盘旋转带来的离心加速度和重力加速度同量级甚至会出现混沌特征小角度模型完全失效。我在第二次重构模拟时就踩了这个坑把初始振幅设为 0.5 米摆长只有 0.8 米结果轨迹图出现一整圈花瓣状张嘴曲线和理论预测完全对不上。后来换回小角度条件结果立刻吻合。所以做这个模拟先统一确认实验条件再决定用哪套方程。面向演示场景小角度足够面向工程仿真强烈建议直接上完整三维摆球模型。3. 数值求解从方程到可跑的代码3.1 参数选取和单位制模拟的输入参数一共四个摆长 L、重力加速度 g、圆盘角速度 Ω、初始状态。推荐一套容易出效果但数值也稳的参数L 1.0 米g 9.81 m/s²Ω 0.5 rad/s初始位置 x(0) 0.2 米y(0) 0初始速度全部为 0。此时自然频率约 3.13 rad/s进动周期 2π / Ω 约 12.57 秒摆动周期约 2 秒。你可以在一个进动周期内看到大约 6 个往复振荡画面既不过密也不过疏。角速度大于 1.2 rad/s 时离心项开始明显影响有效重力轨迹会出现外扩倾向需要调小时间步长。角速度小于 0.1 rad/s 时进动太慢得把模拟时长拉到 30 秒以上否则看不出效果。3.2 为什么我不用欧拉法模拟旋转坐标系下的振荡系统最忌讳用简单显式欧拉法。因为系统本身存在守恒性显式欧拉每个时间步都会引入误差能量会缓慢增长几分钟的模拟后摆锤“越摆越疯”轨迹半径明显放大。我在第一个版本里用过欧拉法步长 0.001 秒跑到 20 秒时能量已经变成初值的 1.3 倍这个结果一点参考价值都没有。换成经典四阶龙格库塔RK4后同样的步长能量漂移可以控制在千分之一以内轨迹闭合得很好。如果你不想手写直接用 scipy 的 solve_ivp指定 methodRK45 或者默认方法都可以。我为了讲清楚数值稳定性下面是手写 RK4 的版本。3.3 状态向量和驱动函数把状态定义为 s [x, y, vx, vy]右侧函数返回 [vx, vy, ax, ay]。这一步虽然简单但坐标顺序必须和方程完全对应。完整的 RK4 流程需要计算四个 k 向量然后把加权平均作为增量。代码如下import numpy as np import matplotlib.pyplot as plt L 1.0 g 9.81 Omega 0.5 omega_n np.sqrt(g / L) def rhs(t, state): x, y, vx, vy state ax -omega_n**2 * x Omega**2 * x 2 * Omega * vy ay -omega_n**2 * y Omega**2 * y - 2 * Omega * vx return np.array([vx, vy, ax, ay]) def rk4_step(f, t, state, dt): k1 f(t, state) k2 f(t 0.5 * dt, state 0.5 * dt * k1) k3 f(t 0.5 * dt, state 0.5 * dt * k2) k4 f(t dt, state dt * k3) return state (dt / 6.0) * (k1 2 * k2 2 * k3 k4) dt 0.005 t_max 16.0 n_steps int(t_max / dt) state np.array([0.2, 0.0, 0.0, 0.0]) traj [state.copy()] t_cur 0.0 for _ in range(n_steps): state rk4_step(rhs, t_cur, state, dt) t_cur dt traj.append(state.copy()) traj np.array(traj)跑完后 traj 的每一行是 [x, y, vx, vy]。注意这里坐标系是圆盘参考系画图直接取前两列就是盘面上摆锤的轨迹。4. 结果怎么看轨迹、进动和两种视角4.1 玫瑰线图形背后的频率关系运行上面的代码画 x-y 轨迹你会看到摆锤从 (0.2, 0) 出发画出一系列弯曲的环。由于进动角速度 Ω 远小于摆动频率 ω轨迹看起来像一条缓慢旋转的直线拉出的花瓣。当 Ω 和 ω 的比值接近简单分数时例如 ω/Ω 接近整数轨迹会闭合形成漂亮的玫瑰花状比值偏离整数时图案会逐渐填充一个圆环区域。这个现象用频域解释非常直观摆动本身是频率 ω 的振荡进动让振荡轴以角速度 Ω 旋转组合后的运动在两个方向上有不同的特征频率轨迹形状取决于 ω 和 Ω 的比例。我实测了一个数据点L 1.0g 9.81ω/Ω ≈ 6.26因为比值接近 6.26 不是精确整数轨迹在 16 秒内还没有闭合但已经能看到清晰的六瓣花型。4.2 圆盘视角和地面视角的对比摆锤在地面惯性系中的坐标和圆盘坐标之间是一个旋转矩阵关系。设圆盘坐标系中的位置为 (x, y)地面系坐标 (X, Y) 满足X x cos(Ωt) - y sin(Ωt) Y x sin(Ωt) y cos(Ωt)把两种视角的轨迹画在同一个图里对比是理解这个模拟最核心的一步。旋转盘上看摆动平面在逆时针进动地面系上看摆动平面几乎完全不动。这个“不动”才是真实物理惯性空间的表现。如果用傅科摆类比地面系就是地球外的惯性参考圆盘就是自转的地球摆锤振荡平面稳定地表观测者觉得它在转。实现这段对比代码很简单在 RK4 循环里同步计算惯性系坐标即可inertial_x x * np.cos(Omega * t) - y * np.sin(Omega * t) inertial_y x * np.sin(Omega * t) y * np.cos(Omega * t)4.3 画出随时间变化的进动角进动角可以直接从轨迹坐标反推每一时刻摆球相对盘心的方向角是 θ(t) atan2(y, x)。对一条近似线性的轨迹θ 并不完全等于摆动平面角度因为振幅会周期性过零点过零点时的 atan2 会跳变。更好的做法是记录轨迹长轴的方向或者对相邻极值点拟合直线。工程实现上我选择每摆动半周期找一次轨迹离原点最远的点这些远点连线方向代表该时刻的摆动平面方向。把这组方向角对时间作散点图再叠加一条斜率 -Ω 的直线你会看到模拟点和理论线高度重合。这也是验证方程与代码正确性的黄金标准。如果符号反了斜率变成正 Ω立即能发现科里奥利项写错了。5. 常见问题与现场排错我踩过的坑都在这5.1 快速定位问题速查遇到问题的操作顺序其实有规律可循。我把模拟过程中最常出现的几类和对应修复列成了一张表现象可能原因检查与修复轨迹半径越来越大欧拉法或步长太大换 RK4dt 降到 0.005 以下进动方向与理论相反科里奥利项符号写反检查叉乘顺序确认 Ω 正方向轨迹高频抖动时间步长过大频率分辨率不够dt 至少小于振荡周期的 1/500摆球快速飞向盘边缘离心项大于回复项降低 Ω确保 Ω 远小于 sqrt(g/L)图形在某个方向被拉伸初始速度非零导致额外漂移初始速度设为 0或调整释放条件地面系和盘面系轨迹完全一致忘记做坐标旋转检查旋转矩阵变换是否接入绘图5.2 能量漂移的排查运行 20 秒以上的模拟能量是一个很好的探测指标。小角度无阻尼摆的能量应该严格守恒数值方法如果稳定能量曲线基本是一条水平直线。我早期用 dt 0.02 跑 RK420 秒内能量漂移约 0.5%视觉上轨迹虽看不出问题但放大频率谱会看到伪峰。把 dt 降到 0.005 后漂移降到 0.02% 以下完全可忽略。检查能量时注意旋转参考系下的“机械能”并不等于地面系机械能因为离心力引入的势能是 -mΩ²r²/2。所以如果要讨论守恒性最好切换到地面惯性系计算总能量否则会因为坐标变换引入虚假的增减项一度让我误判是数值发散。5.3 让人挠头的科里奥利符号科里奥利力的叉乘方向是初学者最容易搞错的地方。我记得某次重新推导时把 Ω × v′ 和 v′ × Ω 搞混了结果所有轨迹图对称翻转又花了半天排查。提供一个可靠记忆法在北半球运动的物体总被“推向运动方向的右侧”科里奥利力等效为一种偏转力。在这套公式里如果 Ω 向上且 v′ 沿 x 正向科里奥利加速度的 y 分量是负的对应物体向 y 负向偏转也就是右侧。如果你设定的坐标系 x 向右、y 向上右手定则核对一遍再从数学上验证能省很多冤枉时间。6. 从这个小模型跳到更广的物理空间6.1 加上阻尼会怎么样摆和圆盘之间如果引入空气阻尼运动方程会多一项 -γ v′。这会让轨迹逐渐收敛摆动平面进动也不再维持恒定转速会出现衰减包络。工程上这项很有用因为实际装置不可能完全无阻尼预测残余振动的持续时间需要它。把阻尼系数设为 0.05 s⁻¹轨迹在 20 秒后几乎缩到原点附近进动周期也变模糊了。这说明阻尼不仅影响振幅还影响进动的可观测性。6.2 从圆盘到真实世界的三类延伸第一类是傅科摆教学演示地球上 Ω 很小摆长又特别长进动周期可以达到十几个小时。这里的圆盘模型瞬间把时间尺度压缩到十几秒非常适合课堂演示。第二类是旋转机械中的摆式减振器例如风力机叶片或直升机旋翼上的悬挂质量会利用旋转系中的科里奥利效应抵消振动。在这种场景下摆的安装点不在旋转轴上而是在半径 r₀ 处方程里还要加一项 r₀ 相关的几何耦合模拟方式可以在这个 μJavaScript基础上直接扩展。第三类是太空飞行器姿态控制中的角动量交换装置旋转参考系下的动力学方程结构类似只是重力回复项被弹簧或铰链力矩替代。你在圆盘摆动里积累的 RK4 和参考系理解迁移过去基本无缝衔接。6.3 还能改什么值得一试的变体想玩得更深入试三个变体。变体一是摆长与半径同量级的三维球摆模型去掉小角度近似使用球坐标或者四元数表示姿态轨迹会从平面进动变成完整的三维曲面进动视觉效果猛增。变体二是在圆盘上加一个径向槽把摆球约束在槽内此时只保留一维运动科里奥利力变成一种“等效阻尼”这属于旋转坐标系下另一类经典问题。变体三是给支点施加周期性垂直振动相当于参数激励系统可以在特定 Ω 区间出现摆动共振和小孩荡秋千的原理同源。每个变体都建议保留和我一样的调试流程先杜绝数值误差再对照解析解最后才调到花哨参数上去看现象。这样无论模型多复杂你都有一条可靠的排查路径。我个人最后再啰嗦一句做旋转圆盘摆动这类模拟最难的不是会调 solve_ivp而是时刻清楚你正站在哪个参考系里描述数据。圆盘坐标、地面坐标、惯性空间坐标这三个视角中的同一段轨迹可能画出来完全不一样但物理本质只有一条。把这条主线抓住了后续不管拓展到多复杂的旋转系统都不会慌。
