SymPy 物理力学教程：用 WrappingCylinder 与 WrappingPathway 建模阿特伍德机并验证绳索不可伸长-开发者社区

SymPy 物理力学教程：用 WrappingCylinder 与 WrappingPathway 建模阿特伍德机并验证绳索不可伸长

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy

本篇技术指南基于 SymPy 官方教程《Atwood Machine Example》（原始文档），完整演示如何在sympy.physics.mechanics框架下用WrappingCylinder表示滑轮、用WrappingPathway表示绕过滑轮的绳索，通过 Kane 方法推导阿特伍德机（Atwood machine）的运动方程，并自动验证绳索的不可伸长约束。读完本文，你将掌握 wrapping geometry 的完整建模流程：从质点位置设定、切点选取、绳长计算到力的生成与运动方程求解，并理解底层geodesic_length与to_loads的实现原理。

问题背景：为什么要用 wrapping pathway 建模阿特伍德机

阿特伍德机由两个质量分别为 $m_1$、$m_2$ 的质点组成，它们通过一根无质量、不可伸长的绳索连接，绳索跨过一个半径为 $r$ 的固定滑轮。当一个质量下降位移 $q$ 时，另一个质量上升同样的位移，因此绳索总长保持不变。

传统教材中这类问题通常手动写出约束方程，而 SymPy 的 wrapping_geometry.py 与 pathway.py 提供了专门用于沿曲面建模力的类：

WrappingCylinder：将滑轮抽象为理想圆柱体；
WrappingPathway：将绳索抽象为绕过该圆柱体的路径，自动计算测地线长度与端点受力。

本示例展示如何建立运动学、计算绳总长、验证不可伸长性，并最终用 Kane 方法推导运动方程。

定义变量与导入模块

首先导入所需的符号、参考系、点与类。系统只有一个广义坐标 $q(t)$，表示 $m_1$ 向下移动的位移；其时间导数 $u(t)=\dot q(t)$ 为广义速度。

>>> import sympy as sp >>> from sympy import Q, refine >>> from sympy.physics.mechanics import ( ... ReferenceFrame, ... Point, ... Particle, ... KanesMethod, ... dynamicsymbols, ... WrappingCylinder, ... WrappingPathway, ... Force, ... ) >>> >>> t = sp.symbols("t") >>> >>> # 常量参数：质量、重力加速度、滑轮半径、初始高度、绳张力 >>> >>> m1, m2, g, r, h, T = sp.symbols("m1 m2 g r h T", positive=True, real=True) >>> >>> # 广义坐标与广义速度：m1 沿 z 负方向下移 q(t) >>> >>> q = dynamicsymbols("q", real=True) # q(t) >>> u = dynamicsymbols("u", real=True) # u(t) = dq/dt

要点说明：

常量参数全部声明为positive=True, real=True，这为后续refine化简（如 $\sqrt{(h+q)^2}=h+q$）提供了必要的假设前提；
dynamicsymbols生成随时间变化的符号 $q(t)$，SymPy 会将其自动识别为关于 $t$ 的时变函数；
Force类定义于 loads.py，用于描述作用在点上的力（同时存储作用点与力向量）。

定义惯性参考系与滑轮中心

定义惯性参考系 $N$，将滑轮中心 $O$ 固定在原点，滑轮转轴沿 $\hat{\mathbf{N}}_x$ 方向：

>>> # 定义惯性参考系与滑轮中心 >>> >>> N = ReferenceFrame("N") >>> O = Point("O") >>> O.set_vel(N, 0) # 滑轮中心固定于原点

两个质点的位置与速度

设 $P_1$ 为 $m_1$ 的接触质点，$P_2$ 为 $m_2$ 的接触质点。初始时（$q=0$）两质点均在滑轮中心下方垂直距离 $h+r$ 处。当 $m_1$ 下移 $q$ 后：

$$ P_1: (x=0,; y=+r,; z=-(h+q)),, $$

$$ P_2: (x=0,; y=-r,; z=-(h-q)),. $$

在参考系 $N$ 中对位置向量求导即可得到速度：

>>> # 质量 m1 位于 P1：(x=0, y=+r, z=-(h+q)) >>> >>> P1 = Point("P1") >>> P1.set_pos(O, r * N.y + (-(h + q)) * N.z) >>> P1.vel(N) - Derivative(q(t), t)*N.z >>> M1 = Particle("M1", P1, m1) >>> >>> # 质量 m2 位于 P2：(x=0, y=-r, z=-(h-q)) >>> >>> P2 = Point("P2") >>> P2.set_pos(O, -r * N.y + (-(h - q)) * N.z) >>> P2.vel(N) Derivative(q(t), t)*N.z >>> M2 = Particle("M2", P2, m2)

从P1.vel(N)的输出可以看到，$m_1$ 的速度为 $-\dot q,\hat{\mathbf{N}}_z$（向下），$m_2$ 的速度为 $+\dot q,\hat{\mathbf{N}}_z$（向上），两者大小相等、方向相反，这正是不可伸长绳索的体现。

用 WrappingCylinder 建模滑轮

将滑轮抽象为半径 $r$、中心在 $O$、转轴沿 $\hat{\mathbf{N}}_x$ 的理想圆柱体：

>>> pulley = WrappingCylinder(r, O, N.x)

WrappingCylinder定义于 wrapping_geometry.py，其构造函数接收三个参数（见 wrapping_geometry.py）：

参数	类型	说明
`radius`	`Symbol`	圆柱半径，必须是正的常值符号（不能是动力学符号）
`point`	`Point`	圆柱轴线经过的点（此处为滑轮中心 $O$）
`axis`	`Vector`	圆柱轴线方向向量（构造时内部会调用`axis.normalize()`归一化）

WrappingCylinder继承自抽象基类WrappingGeometryBase（wrapping_geometry.py），该基类定义了统一的接口契约（point、geodesic_length、geodesic_end_vectors等抽象成员），用户也可以通过子类化创建自定义几何类型。

确定切点 T1 与 T2

由于每个质量都悬挂在滑轮最外侧/最内侧点正下方（即圆柱面上 $y=\pm r$、$z=0$），切点固定不变：

$$ T_1: (x=0,, y=+r,, z=0), \quad T_2: (x=0,, y=-r,, z=0). $$

在两点处放置点 $T_1$、$T_2$ 并令其速度为零：

>>> # P1 对应的切点（最外侧） >>> >>> T1 = Point("T1") >>> T1.set_pos(O, r * N.y + 0 * N.z) >>> T1.set_vel(N, 0) >>> >>> # P2 对应的切点（最内侧） >>> >>> T2 = Point("T2") >>> T2.set_pos(O, -r * N.y + 0 * N.z) >>> T2.set_vel(N, 0)

创建 WrappingPathway

有了两个切点 $T_1$、$T_2$ 和WrappingCylinder滑轮对象，即可构造WrappingPathway：

>>> wpath = WrappingPathway(T1, T2, pulley)

WrappingPathway定义于 pathway.py，其构造参数为：

参数	类型	说明
`attachment_1`	`Point`	路径的第一个端点（绳索一端锚点）
`attachment_2`	`Point`	路径的第二个端点
`geometry`	`WrappingGeometryBase`	路径所绕的几何体

从 pathway.py 的源码可以看出两点约束：

attachment_1、attachment_2必须恰好是两个Point（数量与类型错误都会抛出TypeError）；
geometry必须是WrappingGeometryBase的实例，且构造后不可变（再次赋值会抛出AttributeError）。

在内部，该对象会自动计算圆柱面上连接 $T_1$、$T_2$ 的测地线（最短路径）。在本例中，两点恰好在圆柱横截面的正对两侧，因此测地线就是半圆周长 $\pi r$，与 $q$ 无关。

计算各段绳长并验证不可伸长性

本节仅用于演示WrappingPathway的能力，并非得到正确加速度结果所必需。

绳索由三段组成：

段 1：从 $P_1$ 到 $T_1$ 的竖直段，长度 $L_1 = |P_1 - T_1| = h + q$（利用 $h+q>0$ 的假设化简）；
圆弧段：沿滑轮表面的 $T_1$ 到 $T_2$ 段，即测地线半圆周 $L_\text{curve} = \pi r$；
段 2：从 $T_2$ 到 $P_2$ 的竖直段，长度 $L_2 = |P_2 - T_2| = h - q$（利用 $h-q>0$ 的假设化简）。

总绳长：

$$ L_{\text{total}} = L_{1} + L_{\text{curve}} + L_{2} = (h + q) + \pi r + (h - q) = 2h + \pi r $$

与 $q$ 无关，因此 $\dfrac{dL_{\text{total}}}{dq}=0$。代码如下：

>>> # 段长：P1 到 T1 >>> >>> L1 = sp.sqrt((P1.pos_from(T1).dot(P1.pos_from(T1)))) >>> L1 = refine(L1, Q.positive(h + q)) # 强制假设 h+q > 0 >>> L1 h + q(t) >>> >>> # 段长：P2 到 T2 >>> >>> L2 = sp.sqrt((P2.pos_from(T2).dot(P2.pos_from(T2)))) >>> L2 = refine(L2, Q.positive(h - q)) # 强制假设 h-q > 0 >>> L2 h - q(t) >>> >>> # 滑轮上的圆弧段 >>> >>> L_curve = wpath.length >>> L_curve pi*r >>> >>> # 总长及其对 q 的导数 >>> >>> L_total = sp.simplify(L1 + L_curve + L2) >>> L_total 2*h + pi*r >>> dL_dq = sp.simplify(sp.diff(L_total, q)) >>> dL_dq 0

这里的关键是wpath.length。从源码看（pathway.py），WrappingPathway.length直接委托给几何体的geodesic_length(*self.attachments)：

@property def length(self): """Exact analytical expression for the pathway's length.""" return self.geometry.geodesic_length(*self.attachments)

而WrappingCylinder.geodesic_length（wrapping_geometry.py）用勾股定理计算测地线：一个直角边是两点沿圆柱轴线的平行距离，另一个直角边是两点在圆柱横截面上的圆弧长度（$r \times$ 中心角），斜边即测地线长度。由于 $T_1$、$T_2$ 在轴线方向无偏移，中心角为 $\pi$，故 $L_\text{curve} = r \cdot \pi$。

在 test_pathway.py 的test_static_pathway_on_cylinder_length测试中，可以找到对圆柱面上静态路径长度的参数化验证，例如 $(1,0,0)$ 与 $(0,1,0)$ 两个端点对应测地线长度 $\frac{1}{2}\pi r$，$(1,0,0)$ 与 $(-1,0,0)$ 对应 $\pi r$，与本例的半圆周结论相互印证。

定义重力载荷

每个质点受沿 $-\hat{\mathbf{N}}_z$ 方向的重力：

>>> grav1 = Force(P1, -m1 * g * N.z) >>> grav2 = Force(P2, -m2 * g * N.z)

收集所有载荷：to_loads 自动生成绳张力

系统中唯一的广义坐标是 $q$ 及其导数 $u$。绳索通过 wrapping pathway 将张力 $T$ 传递给两个质量。调用wpath.to_loads(T)会自动得到三个Force对象：

在 $P_1$ 处沿切向拉动质量 $m_1$ 的力；
在 $P_2$ 处拉动质量 $m_2$ 的力；
在滑轮中心 $O$ 处的等大反向反作用力。

>>> loads = wpath.to_loads(T) + [grav1, grav2]

从 pathway.py 的to_loads实现可以看到其内部逻辑：

pA, pB = self.attachments pO = self.geometry.point pA_force, pB_force = self.geometry.geodesic_end_vectors(pA, pB) pO_force = -(pA_force + pB_force) loads = [ Force(pA, force * pA_force), Force(pB, force * pB_force), Force(pO, force * pO_force), ]

即：两端点的力方向由几何体在端点处的测地线端向量geodesic_end_vectors决定（wrapping_geometry.py），而滑轮中心的反作用力恰好等于两端点力的负和，从而保证整个系统的合力与合力矩自洽。该行为在 test_pathway.py 的test_static_pathway_on_cylinder_to_loads测试中被逐一验证（例如端点 $(1,0,0)$ 与 $(0,1,0)$ 对应 $pA$ 受 $F\hat{\mathbf{N}}_y$、$pB$ 受 $F\hat{\mathbf{N}}_x$、$pO$ 受 $-F(\hat{\mathbf{N}}_x+\hat{\mathbf{N}}_y)$）。

建立运动学微分方程

声明通常的运动学关系 $u = \dot q$：

>>> kin_diff = [u - q.diff()]

用 Kane 方法建立并求解运动方程

以惯性系 $N$、一个坐标 $q$、一个速度 $u$ 及运动学关系 $u-\dot q=0$ 构造KanesMethod，两个质点M1、M2与loads列表描述了系统中的全部力：

>>> kane = KanesMethod(N, (q,), (u,), kd_eqs=kin_diff) >>> bodies = [M1, M2] >>> Fr, Frs = kane.kanes_equations(bodies, loads)

求解 $\ddot q$（即 $\dot u$）关于 $q$、$u$ 和 $T$ 的表达式。由于 $T$ 是未知反力，符号结果中会包含 $T$。化简后得到标准的二阶运动方程：

>>> [u, u_dot] = kane.rhs() >>> qdd = sp.simplify(u_dot) >>> sp.pprint(qdd, use_unicode=True) g⋅(m₁ - m₂) ─────────── m₁ + m₂

这正是阿特伍德机的经典加速度公式 $\ddot q = \dfrac{m_1 - m_2}{m_1 + m_2} g$。值得注意的是，张力 $T$ 在化简后的加速度表达式中被消去，这是质点在竖直方向只受重力与绳张力、且张力做功为零（绳索不可伸长）的必然结果。

数值验证

最后代入 $m_1=1$、$m_2=2$、$g=9.81$、$h=5.0$、$r=0.5$，数值确认 $\ddot q$ 与 $\frac{m_1 - m_2}{m_1 + m_2} g$ 一致：

>>> numeric_vals = {m1: 1.0, m2: 2.0, g: 9.81, h: 5.0, r: 0.5} >>> qdd_num = float(qdd.subs(numeric_vals)) >>> print(f"{qdd_num:.6f} m/s²") -3.270000 m/s²

代入公式验算：$\frac{1.0 - 2.0}{1.0 + 2.0} \times 9.81 = -\frac{9.81}{3} = -3.27$，与输出完全吻合（负号表示 $m_1$ 实际向上运动，因为 $m_2 > m_1$）。

源码级原理小结

本示例背后是sympy.physics.mechanics三个核心文件的协同工作：

几何抽象层wrapping_geometry.py：WrappingGeometryBase定义统一接口，WrappingCylinder实现圆柱面上的测地线长度（勾股定理式分解）与测地线端向量；同类还提供WrappingSphere、WrappingCone等几何体，可覆盖更复杂的曲面缠绕场景。
路径抽象层pathway.py：WrappingPathway组合两个锚点与一个几何体，对外暴露length（测地线长度）、extension_velocity（长度对时间的导数，即伸长速度）与to_loads(force)（生成两端点与几何中心三处载荷）。
动力学求解层kane.py（KanesMethod）与 loads.py：消费上述载荷列表，自动建立广义力并输出运动方程。

这种“几何 + 路径 + 求解器”的分层设计，使得同一个WrappingPathway可以无缝替换几何体（圆柱、球、锥），而无需改动载荷生成与运动方程求解的代码——这是本教程所选建模方式的扩展价值所在。

结论

本教程完整演示了在 SymPy 力学框架中建模阿特伍德机的过程：

用WrappingCylinder表示滑轮，用WrappingPathway表示绕过滑轮的绳索；
通过wpath.length自动计算圆弧段绳长，结合竖直段长度验证总绳长 $L_\text{total} = 2h + \pi r$ 与广义坐标 $q$ 无关，从而自动验证绳索不可伸长约束；
用wpath.to_loads(T)自动生成两端点与滑轮中心的三处张力载荷，再叠加重力；
用 Kane 方法推导出经典二阶运动方程，并恢复经典加速度公式 $\ddot q = \frac{m_1 - m_2}{m_1 + m_2} g$，数值验证一致。