
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, positiveTrue, realTrue) # 广义坐标与广义速度m1 沿 z 负方向下移 q(t) q dynamicsymbols(q, realTrue) # q(t) u dynamicsymbols(u, realTrue) # u(t) dq/dt要点说明常量参数全部声明为positiveTrue, realTrue这为后续refine化简如 $\sqrt{(hq)^2}hq$提供了必要的假设前提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$ 的接触质点。初始时$q0$两质点均在滑轮中心下方垂直距离 $hr$ 处。当 $m_1$ 下移 $q$ 后$$ P_1: (x0,; yr,; z-(hq)),, $$$$ P_2: (x0,; y-r,; z-(h-q)),. $$在参考系 $N$ 中对位置向量求导即可得到速度 # 质量 m1 位于 P1(x0, yr, z-(hq)) 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(x0, 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参数类型说明radiusSymbol圆柱半径必须是正的常值符号不能是动力学符号pointPoint圆柱轴线经过的点此处为滑轮中心 $O$axisVector圆柱轴线方向向量构造时内部会调用axis.normalize()归一化WrappingCylinder继承自抽象基类WrappingGeometryBasewrapping_geometry.py该基类定义了统一的接口契约point、geodesic_length、geodesic_end_vectors等抽象成员用户也可以通过子类化创建自定义几何类型。确定切点 T1 与 T2由于每个质量都悬挂在滑轮最外侧/最内侧点正下方即圆柱面上 $y\pm r$、$z0$切点固定不变$$ T_1: (x0,, yr,, z0), \quad T_2: (x0,, y-r,, z0). $$在两点处放置点 $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_1Point路径的第一个端点绳索一端锚点attachment_2Point路径的第二个端点geometryWrappingGeometryBase路径所绕的几何体从 pathway.py 的源码可以看出两点约束attachment_1、attachment_2必须恰好是两个Point数量与类型错误都会抛出TypeErrorgeometry必须是WrappingGeometryBase的实例且构造后不可变再次赋值会抛出AttributeError。在内部该对象会自动计算圆柱面上连接 $T_1$、$T_2$ 的测地线最短路径。在本例中两点恰好在圆柱横截面的正对两侧因此测地线就是半圆周长 $\pi r$与 $q$ 无关。计算各段绳长并验证不可伸长性本节仅用于演示WrappingPathway的能力并非得到正确加速度结果所必需。绳索由三段组成段 1从 $P_1$ 到 $T_1$ 的竖直段长度 $L_1 |P_1 - T_1| h q$利用 $hq0$ 的假设化简圆弧段沿滑轮表面的 $T_1$ 到 $T_2$ 段即测地线半圆周 $L_\text{curve} \pi r$段 2从 $T_2$ 到 $P_2$ 的竖直段长度 $L_2 |P_2 - T_2| h - q$利用 $h-q0$ 的假设化简。总绳长$$ 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)) # 强制假设 hq 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.pyWrappingPathway.length直接委托给几何体的geodesic_length(*self.attachments)property def length(self): Exact analytical expression for the pathways length. return self.geometry.geodesic_length(*self.attachments)而WrappingCylinder.geodesic_lengthwrapping_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 q0$ 构造KanesMethod两个质点M1、M2与loads列表描述了系统中的全部力 kane KanesMethod(N, (q,), (u,), kd_eqskin_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_unicodeTrue) g⋅(m₁ - m₂) ─────────── m₁ m₂这正是阿特伍德机的经典加速度公式 $\ddot q \dfrac{m_1 - m_2}{m_1 m_2} g$。值得注意的是张力 $T$ 在化简后的加速度表达式中被消去这是质点在竖直方向只受重力与绳张力、且张力做功为零绳索不可伸长的必然结果。数值验证最后代入 $m_11$、$m_22$、$g9.81$、$h5.0$、$r0.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.pyWrappingGeometryBase定义统一接口WrappingCylinder实现圆柱面上的测地线长度勾股定理式分解与测地线端向量同类还提供WrappingSphere、WrappingCone等几何体可覆盖更复杂的曲面缠绕场景。路径抽象层pathway.pyWrappingPathway组合两个锚点与一个几何体对外暴露length测地线长度、extension_velocity长度对时间的导数即伸长速度与to_loads(force)生成两端点与几何中心三处载荷。动力学求解层kane.pyKanesMethod与 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$数值验证一致。相关参考原始教程文档atwoods_machine_example.rst本教程配套示意图见 atwood_machine.svg力学教程索引mechanics/index.rst内含滚动圆盘、非最小坐标摆、多自由度完整约束系统等更多 Kane 方法示例几何实现wrapping_geometry.py路径实现pathway.py载荷实现loads.py测试用例test_pathway.py、test_wrapping_geometry.py【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考