SymPy 多体动力学实战：用 Kane 方法从接触点建模滚动圆盘（Rolling Disc）

SymPy 多体动力学实战用 Kane 方法从接触点建模滚动圆盘Rolling Disc【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy导读本文围绕 SymPysympy.physics.mechanics模块中的经典教学案例——滚动圆盘rolling disc展开完整讲解如何从接触点contact point向上构造运动学避免引入额外的广义速度仅用 3 个广义坐标、3 个广义速度、圆盘质量/半径与局部重力即可用 Kane 方法得到运动方程。读完本文你将掌握ReferenceFrame.orientnew的 3-1-2Z、X、Y简单旋转序列、Point.v2pt_theory二点速度理论、inertia惯性并矢构造、KanesMethod求解Fr Fr* 0的完整流程并能基于源码理解质量矩阵与 forcing 的底层形成机制。本文依据仓库文档 rollingdisc_example_kane.rst 编写并结合 kane.py 源码与 test_kane.py 测试用例进行纵深佐证。1. 物理模型与建模策略滚动圆盘案例的物理假设为圆盘无限薄infinitely thin与地面仅在1 个点接触圆盘在地面上无滑动滚动rolling without slip。SymPy 官方教程将同一系统用三种方法建模用于展示sympy.physics.mechanics模块的多种功能见索引文档 rollingdisc_example.rstKane 方法本文主题rollingdisc_example_kane.rstKane 方法 约束力auxiliary speedsrollingdisc_example_kane_constraints.rstLagrange 方法rollingdisc_example_lagrange.rst。本案例最关键的设计决策是运动学从接触点接触点速度恒为零向上构造这样就“不需要引入广义速度”removing the need to introduce generalized speeds。因此描述整个系统只需要3 个构型变量configuration variablesq1, q2, q33 个速度变量speed variablesu1, u2, u3圆盘质量m、半径r以及局部重力g。值得注意的一点是质量m最终会从运动方程中消去note that mass will drop out因为重力与惯性力均与质量成正比。这在最终的rhs表达式中可以直观验证——三个方程均不含m。2. 环境准备符号、动力学符号与输出设置建模的第一步是导入所需符号与力学模块。原文档给出的代码如下 from sympy import symbols, sin, cos, tan from sympy.physics.mechanics import * q1, q2, q3, u1, u2, u3 dynamicsymbols(q1 q2 q3 u1 u2 u3) q1d, q2d, q3d, u1d, u2d, u3d dynamicsymbols(q1 q2 q3 u1 u2 u3, 1) r, m, g symbols(r m g) mechanics_printing(pretty_printFalse)关键点逐一说明from sympy.physics.mechanics import *导入了ReferenceFrame、Point、dynamicsymbols、RigidBody、KanesMethod、inertia、dot等符号dynamicsymbols(q1 q2 q3 u1 u2 u3)创建的是关于时间t的动力学符号即q1(t)这种函数形式第二个参数1表示创建它们的一阶时间导数q1d(t)等mechanics_printing(pretty_printFalse)关闭 pretty printing使后续输出为常规 ASCII 矩阵形式。从源码 kane.py 可以看到KanesMethod内部通过dynamicsymbols._t表示时间如self._qdot self.q.diff(dynamicsymbols._t)所有动力学符号的时间微分都在此基础上进行因此在使用时必须保持符号的“动力学”属性否则q1.diff(t)恒为 0。3. 旋转运动学3-1-2Z、X、Y简单旋转序列旋转运动学通过一系列简单旋转simple rotations构造。每一次简单旋转创建一个新参考系下一次旋转基于新参考系的基向量定义。本案例采用3-1-2 旋转序列即依次绕Z、X、Y轴旋转 N ReferenceFrame(N) Y N.orientnew(Y, Axis, [q1, N.z]) L Y.orientnew(L, Axis, [q2, Y.x]) R L.orientnew(R, Axis, [q3, L.y]) w_R_N_qd R.ang_vel_in(N) R.set_ang_vel(N, u1 * L.x u2 * L.y u3 * L.z)各参考系的物理含义N惯性参考系地面Y绕N.z旋转角度q1得到的“偏航yaw”系L绕Y.x旋转角度q2得到的“侧倾/lean”系圆盘倾倒角R绕L.y旋转角度q3得到的“滚动/roll”系圆盘绕自身轴自转。orientnew(Y, Axis, [q1, N.z])的含义是以Axis方式、绕N.z轴旋转q1角度创建新参考系Y其余依此类推。特别需要强调的是角速度的表示方式 R.set_ang_vel(N, u1 * L.x u2 * L.y u3 * L.z)圆盘相对惯性系的角速度w_R_N直接用第二个参考系Llean frame的基向量分解为三个广义速度分量u1, u2, u3。这正是原文档强调“需要定义中间参考系、而非直接使用 body-three 姿态”的原因——若使用 body-three 姿态角速度表达式的符号形式会复杂得多。而先通过w_R_N_qd R.ang_vel_in(N)计算出“由坐标时间导数表达的角速度”是为了后面构造运动微分方程kd时做对照。4. 平动运动学从接触点向上构造平动运动学是本文档最核心的建模技巧。步骤为创建接触点C并显式设置其在惯性系N中的速度为零用locatenew构造从接触点到圆盘质心Dmc的位置向量沿L.z方向长度为半径r用v2pt_theory由C的速度直接推导Dmc的速度。 C Point(C) C.set_vel(N, 0) Dmc C.locatenew(Dmc, r * L.z) Dmc.v2pt_theory(C, N, R) r*u2*L.x - r*u1*L.yv2pt_theory(otherpoint, outframe, interframe)是 SymPy 力学模块的“二点理论”已知一点速度、两点位置向量与中间参考系的角速度即可推出另一点速度其数学依据为v_Dmc v_C w_R_N × r_C→Dmc。输出r*u2*L.x - r*u1*L.y即质心速度在L系中的分解——只与u1, u2有关因为L.z方向法向的分量由滚动约束自动满足。这一节展示了“从接触点向上”建模的优势由于接触点速度已知为零无需像其他建模方式那样引入额外的约束变量Dmc的速度自然满足无滑动滚动条件。5. 惯性张量用inertia()构造惯性并矢圆盘以Llean frame为参考系定义惯性张量 I inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2) mprint(I) m*r**2/4*(L.x|L.x) m*r**2/2*(L.y|L.y) m*r**2/4*(L.z|L.z)inertia(frame, ixx, iyy, izz, ixy0, iyz0, izx0)定义于 functions.py其签名默认将非对角项置零。对于无限薄圆盘绕盘面内直径L.x、L.z的转动惯量为m*r**2/4绕盘面法向对称轴L.y的转动惯量为m*r**2/2。原文档指出圆盘滚动过程中其惯性张量在 lean frame 中不发生变化这使得最终方程更简洁——因为惯性并矢在L系中是对角的且为常量不会引入额外的时间导数项。mprint输出中的(L.x|L.x)等即并矢dyadic的标准打印形式。6. 运动微分方程Kinematic Differential Equations接下来建立广义坐标时间导数与广义速度之间的关系 kd [dot(R.ang_vel_in(N) - w_R_N_qd, uv) for uv in L]这段代码的含义圆盘角速度有两种表达——由广义速度定义的R.ang_vel_in(N)即u1*L.x u2*L.y u3*L.z和由坐标时间导数定义的w_R_N_qd。二者做差后分别点乘L.x, L.y, L.z得到三个标量方程构成运动微分方程组kd即速度定义与坐标导数之间的一致性约束形式为f(q_dot, u, q) 0且对q_dot与u均线性。从源码 kane.py 的_initialize_kindiffeq_matrices可以看到KanesMethod会将这些方程线性化为标准形式k_ku(q,t)*u(t) k_kqdot(q,t)*q(t) f_k(q,t) 0并通过linear_eq_to_matrix提取系数矩阵k_ku与k_kqdot再求解得到qdot_u_map即q到u的映射字典。源码还会校验kd在u与q上的线性性若非线性会抛出ValueError提示“The provided kinematic differential equations are nonlinear in ...”。同时源码要求运动微分方程条数与广义坐标数相等否则报错。7. 力列表与刚体定义 ForceList [(Dmc, - m * g * Y.z)] BodyD RigidBody(BodyD, Dmc, R, m, (I, Dmc)) BodyList [BodyD]主动力列表ForceList仅包含作用在质心Dmc上的重力-m*g*Y.z沿惯性竖直方向向下RigidBody(BodyD, Dmc, R, m, (I, Dmc))构造刚体需要指定质心点Dmc、刚体固连参考系R、质量m以及“惯性张量 张量作用点”的二元组(I, Dmc)刚体列表BodyList只含这一个刚体。8. 用KanesMethod组装并求解运动方程最终组装与求解步骤 KM KanesMethod(N, q_ind[q1, q2, q3], u_ind[u1, u2, u3], kd_eqskd) (fr, frstar) KM.kanes_equations(BodyList, ForceList) MM KM.mass_matrix forcing KM.forcing rhs MM.inv() * forcing kdd KM.kindiffdict() rhs rhs.subs(kdd) rhs.simplify() mprint(rhs) Matrix([ [(4*g*sin(q2) 6*r*u2*u3 - r*u3**2*tan(q2))/(5*r)], [ -2*u1*u3/3], [ (-2*u2 u3*tan(q2))*u1]])流程解读KanesMethod(N, q_ind[q1,q2,q3], u_ind[u1,u2,u3], kd_eqskd)指定惯性系、独立广义坐标、独立广义速度与运动微分方程KM.kanes_equations(BodyList, ForceList)内部计算广义主动力Fr与广义惯性力Fr*使其满足 Kane 方程Fr Fr* 0KM.mass_matrix与KM.forcing将动力学方程整理为M * udot forcing的显式线性形式rhs MM.inv() * forcing求解udot广义速度的时间导数kdd KM.kindiffdict()得到运动微分方程的反解字典qdot → f(u)rhs.subs(kdd)用u替换掉q_dot最后simplify()化简。关于KanesMethod的构造参数从 kane.py 的类文档可以了解到其完整签名本文档只用到最核心的几个参数含义本文示例frame惯性参考系Nq_ind独立广义坐标iterable of dynamicsymbols[q1, q2, q3]u_ind独立广义速度[u1, u2, u3]kd_eqs运动微分方程线性联系广义速度与坐标时间导数kdq_dependent依赖广义坐标有约束时未使用configuration_constraints构型完整约束未使用u_dependent/velocity_constraints依赖速度 / 速度约束未使用u_auxiliary辅助广义速度求解非贡献约束力时使用未使用见第 10 节explicit_kinematics是否采用显式运动学默认True默认从源码实现看_form_fr通过partial_velocity计算偏速度并与力做点积累加得到广义主动力Fr_form_frstar将刚体拆分为平动与转动两部分分别累加质量矩阵MM与非线性的nonMM项最终得到fr_star -(MM*udot nonMM)。KM.mass_matrix与KM.forcing即对应k_d与-f_d -(fr - nonMM)。最终得到的三个方程物理含义清晰第 1 个方程含g*sin(q2)与u2*u3、u3²*tan(q2)项是倾斜角q2方向倾覆方向的动力学方程第 2、3 个方程分别对应自转与偏航方向的演化三个方程均不含质量m印证了“质量会消去”的论断。9. 结果验证测试用例与线性化基准仓库测试文件 test_kane.py 中的test_rolling_disc()完整复现了本文档的建模流程并做了两层验证方程正确性断言assert rhs.expand() Matrix([(6*u2*u3*r - u3**2*r*tan(q2) 4*g*sin(q2))/(5*r), -2*u1*u3/3, u1*(-2*u2 u3*tan(q2))]).expand()该断言与原文档mprint(rhs)的输出完全一致可作为你亲手运行代码后的对照基准。线性化临界速度基准经典力学已知结果A KM.linearize(A_and_BTrue)[0] A_upright A.subs({r: 1, g: 1, m: 1}).subs({q1: 0, q2: 0, q3: 0, u1: 0, u3: 0}) assert sympy.sympify(A_upright.subs({u2: 1 / sqrt(3)})).eigenvals() {S.Zero: 6}其物理含义是当r g m 1时直立upright圆盘的临界速度critical speed为1/sqrt(3)——在该速度下线性化状态矩阵的全部 6 个特征值均为 0。测试还通过KM.mass_matrix_full.LUsolve(KM.forcing_full)与KM.rhs()的一致性验证了全维含运动学方程的正确性。此外test_kane.py 中的test_aux()还用另一种方式手动引入辅助速度后再置零验证了辅助速度处理的内在一致性详见下一节。10. 扩展引入约束力auxiliary speeds的 Kane 方法版本姊妹文档 rollingdisc_example_kane_constraints.rst 在本文基础上增加了把非贡献约束力显式带入证据bringing the non-contributing forces into evidence的扩展原理详见 [Kane1985]Kane, T., Levinson, D.Dynamics Theory and Applications, 1985, McGraw-Hill。二者差异要点如下额外引入 3 个辅助速度u4, u5, u6与 3 个约束力分量f1, f2, f3均为动力学符号接触点C的速度不再恒为零而是分解为三个方向法向、滚动路径切向、地面内垂直方向的辅助速度线性组合 C.set_vel(N, u4 * L.x u5 * cross(Y.z, L.x) u6 * Y.z)力列表在重力基础上增加接触点处的约束力 ForceList [(Dmc, - m * g * Y.z), (C, f1 * L.x f2 * cross(Y.z, L.x) f3 * Y.z)]构造KanesMethod时通过u_auxiliary[u4, u5, u6]传入辅助速度 KM KanesMethod(N, q_ind[q1, q2, q3], u_ind[u1, u2, u3], kd_eqskd, ... u_auxiliary[u4, u5, u6])此时动力学主方程rhs与无约束版本完全一致验证了约束力不贡献广义功而约束力由KM.auxiliary_eqs给出源码 kane.py 的auxiliary_eqs属性即“用于求解非贡献力的辅助 Kane 方程组”。化简后得到 from sympy import trigsimp, signsimp, collect, factor_terms def simplify_auxiliary_eqs(w): ... return signsimp(trigsimp(collect(collect(factor_terms(w), f2), m*r))) mprint(KM.auxiliary_eqs.applyfunc(simplify_auxiliary_eqs)) Matrix([ [ -m*r*(u1*u3 u2) f1], [-m*r*u1**2*sin(q2) - m*r*u2*u3/cos(q2) m*r*cos(q2)*u1 f2], [ -g*m m*r*(u1**2*cos(q2) sin(q2)*u1) f3]])其中f1对应滚动方向约束力、f2对应地面内垂直方向约束力、f3对应法向约束力含重力平衡项-g*m。该扩展展示了KanesMethod的u_auxiliary参数在需要求接触力、摩擦力等“被隐藏”的约束反力时的用法。11. 小结与进一步阅读本文基于 rollingdisc_example_kane.rst 完整还原了 SymPy 中 Kane 方法建模滚动圆盘的每一步从接触点向上构造运动学仅用 33 个广义坐标/速度即完整描述无滑动滚动系统3-1-2Z-X-Y简单旋转序列 v2pt_theory二点理论推导质心速度inertia构造对角惯性并矢、RigidBody封装刚体、KanesMethod组装并解出M*udot forcing测试用例给出方程正确性断言与1/sqrt(3)临界速度的线性化基准扩展版本通过u_auxiliary与KM.auxiliary_eqs显式求解约束力。若想进一步对比不同建模方法可继续阅读同一教程目录下的 rollingdisc_example_lagrange.rstLagrange 方法与 rollingdisc_example_kane_constraints.rstKane 方法 约束力以及更一般的多自由度完整约束系统教程 multi_degree_freedom_holonomic_system.rst。KanesMethod的完整参数与属性说明可查阅源码 kane.py 的类文档。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考