企业官网建设流程全解析

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-2（Z、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.rst）：

1. Kane 方法（本文主题，rollingdisc_example_kane.rst）；
2. Kane 方法 + 约束力（auxiliary speeds）（rollingdisc_example_kane_constraints.rst）；
3. Lagrange 方法（rollingdisc_example_lagrange.rst）。

本案例最关键的设计决策是：运动学从接触点（接触点速度恒为零）向上构造，这样就“不需要引入广义速度”（removing the need to introduce generalized speeds）。因此，描述整个系统只需要：

• 3 个构型变量（configuration variables）：q1, q2, q3；
• 3 个速度变量（speed variables）：u1, 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_print=False)

关键点逐一说明：

• 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_print=False)关闭 pretty printing，使后续输出为常规 ASCII 矩阵形式。

从源码 kane.py 可以看到，KanesMethod内部通过dynamicsymbols._t表示时间（如self._qdot = self.q.diff(dynamicsymbols._t)），所有动力学符号的时间微分都在此基础上进行，因此在使用时必须保持符号的“动力学”属性，否则q1.diff(t)恒为 0。

3. 旋转运动学：3-1-2（Z、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直接用第二个参考系L（lean frame）的基向量分解为三个广义速度分量u1, u2, u3。这正是原文档强调“需要定义中间参考系、而非直接使用 body-three 姿态”的原因——若使用 body-three 姿态，角速度表达式的符号形式会复杂得多。而先通过w_R_N_qd = R.ang_vel_in(N)计算出“由坐标时间导数表达的角速度”，是为了后面构造运动微分方程kd时做对照。

4. 平动运动学：从接触点向上构造

平动运动学是本文档最核心的建模技巧。步骤为：

1. 创建接触点C，并显式设置其在惯性系N中的速度为零；
2. 用locatenew构造从接触点到圆盘质心Dmc的位置向量（沿L.z方向，长度为半径r）；
3. 用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.y

v2pt_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()构造惯性并矢

圆盘以L（lean 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, ixy=0, iyz=0, izx=0)定义于 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_eqs=kd) >>> (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]])

流程解读：

1. KanesMethod(N, q_ind=[q1,q2,q3], u_ind=[u1,u2,u3], kd_eqs=kd)：指定惯性系、独立广义坐标、独立广义速度与运动微分方程；
2. KM.kanes_equations(BodyList, ForceList)：内部计算广义主动力Fr与广义惯性力Fr*，使其满足 Kane 方程Fr + Fr* = 0；
3. KM.mass_matrix与KM.forcing：将动力学方程整理为M * udot = forcing的显式线性形式；
4. rhs = MM.inv() * forcing求解udot（广义速度的时间导数）；
5. kdd = KM.kindiffdict()得到运动微分方程的反解字典（qdot → f(u)），rhs.subs(kdd)用u替换掉q_dot，最后simplify()化简。

关于KanesMethod的构造参数，从 kane.py 的类文档可以了解到其完整签名，本文档只用到最核心的几个：

从源码实现看，_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()完整复现了本文档的建模流程，并做了两层验证：

1. 方程正确性断言：

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)的输出完全一致，可作为你亲手运行代码后的对照基准。

2. 线性化临界速度基准（经典力学已知结果）：

A = KM.linearize(A_and_B=True)[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）。二者差异要点如下：

1. 额外引入 3 个辅助速度u4, u5, u6与 3 个约束力分量f1, f2, f3（均为动力学符号）；
2. 接触点C的速度不再恒为零，而是分解为三个方向（法向、滚动路径切向、地面内垂直方向）的辅助速度线性组合：

>>> C.set_vel(N, u4 * L.x + u5 * cross(Y.z, L.x) + u6 * Y.z)

3. 力列表在重力基础上增加接触点处的约束力：

>>> ForceList = [(Dmc, - m * g * Y.z), (C, f1 * L.x + f2 * cross(Y.z, L.x) + f3 * Y.z)]

4. 构造KanesMethod时通过u_auxiliary=[u4, u5, u6]传入辅助速度：

>>> KM = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3], kd_eqs=kd, ... 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 方法建模滚动圆盘的每一步：

• 从接触点向上构造运动学，仅用 3+3 个广义坐标/速度即完整描述无滑动滚动系统；
• 3-1-2（Z-X-Y）简单旋转序列 +v2pt_theory二点理论推导质心速度；
• inertia构造对角惯性并矢、RigidBody封装刚体、KanesMethod组装并解出M*udot = forcing；
• 测试用例给出方程正确性断言与1/sqrt(3)临界速度的线性化基准；
• 扩展版本通过u_auxiliary与KM.auxiliary_eqs显式求解约束力。

若想进一步对比不同建模方法，可继续阅读同一教程目录下的 rollingdisc_example_lagrange.rst（Lagrange 方法）与 rollingdisc_example_kane_constraints.rst（Kane 方法 + 约束力），以及更一般的多自由度完整约束系统教程 multi_degree_freedom_holonomic_system.rst。KanesMethod的完整参数与属性说明可查阅源码 kane.py 的类文档。

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