SymPy 多体动力学实战基于拉格朗日方法推导受迫 Duffing 振子-摆系统的运动方程【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy导读本文基于 SymPy 的sympy.physics.mechanics模块完整演示如何用拉格朗日方法Lagranges Method为一个Duffing 振子悬挂单摆的耦合系统推导运动方程。该系统同时包含非线性弹簧Duffing 弹簧、线性阻尼器、刚体惯量与重力载荷是理解 SymPy 多体动力学建模标准流程——定义符号 → 运动学参考系与点→ 惯量与刚体 → 载荷作动器→ 拉格朗日方程——的理想范例。读完本文你将掌握DuffingSpring、LinearDamper、LinearPathway、RigidBody与LagrangesMethod的组合用法并能自行复现并核对论文中的动力学方程。原文档位于 doc/src/tutorials/physics/mechanics/duffing-example.rst其示例参考了 P. Brzeskia 等人的论文The dynamics of the pendulum suspended on the forced Duffing oscillatorJournal of Sound and Vibration, 2012。系统描述与建模思路本示例建模一个由 Duffing 振子非线性弹簧-质量块与悬挂在其上的单摆组成的系统MDuffing 振子振荡器的质量m摆的质量l摆长质量不计的刚性轻杆k_1、k_2弹簧刚度的线性部分与非线性部分Duffing 弹簧c_1Duffing 振子的黏性阻尼系数。系统采用拉格朗日方程建模。文档特意说明本例在拉格朗日量中将势能置零而把保守力重力与 Duffing 弹簧力显式加入载荷列表loads中这是sympy.physics.mechanics中两种可互换的建模习惯之一。建模整体流程如下定义常量与动力学符号定义参考系惯性系N、摆系B与关键点O、block_point、pendulum_point计算刚体惯量并构建两个RigidBody通过LinearPathwayDuffingSpringLinearDamper生成弹簧与阻尼器载荷并叠加重力载荷计算拉格朗日量L用LagrangesMethod.form_lagranges_equations()得到运动方程。准备工作导入模块示例开头导入 SymPy 与力学模块并开启漂亮打印 import sympy as sm import sympy.physics.mechanics as me me.init_vprinting()sm提供symbols、simplify等符号运算工具me提供ReferenceFrame、Point、RigidBody、LinearPathway、DuffingSpring、LinearDamper、Lagrangian、LagrangesMethod等力学建模类me.init_vprinting()将动力学符号如q1(t)、q1(t)以导数记号漂亮显示。定义变量与参数 M, m, l, k1, k2, c1, g, h, w, d, r sm.symbols(M, m, l, k1, k2, c1, g, h, w, d, r) q1, q2 me.dynamicsymbols(q1 q2) q1d me.dynamicsymbols(q1, 1)各符号的物理含义符号含义hDuffing 振子块的高度用于计算块体惯量wDuffing 振子块的宽度dDuffing 振子块的深度r摆的质量球半径q_1广义坐标表示 Duffing 振子沿N.y的位移q_2广义坐标表示摆相对竖直方向的转角q_1广义速度即q1.diff(t)这里q1, q2 me.dynamicsymbols(q1 q2)等价于sm.symbols(q1 q2, clssm.Function)得到的是关于时间t的符号函数q1(t)、q2(t)而me.dynamicsymbols(q1, 1)返回其一阶导数q1(t)。这正是动力学建模中广义坐标随时间变化的标准化表达。定义运动学参考系与点首先创建惯性参考系N并让摆参考系B绕N.z轴旋转角度q2 N me.ReferenceFrame(N) B N.orientnew(B, axis, (q2, N.z))orientnew(..., axis, (q2, N.z))表示B相对N绕公共轴N.z旋转q2角度。立即可以验证摆的角速度 B.ang_vel_in(N) q2(t) n_z接下来定义系统关键点。O为惯性系中的固定点地面振子块位于O沿N.y方向q1处摆球位于振子块沿B.y方向l处 O me.Point(O) block_point O.locatenew(block, q1 * N.y) pendulum_point block_point.locatenew(pendulum, l * B.y)locatenew(name, pos)在父点基础上相对定位创建新点。然后为各点设置速度 O.set_vel(N, 0) block_point.set_vel(N, q1d * N.y) pendulum_point.v2pt_theory(block_point, N, B) q1(t) n_y -l*q2(t) b_xO.set_vel(N, 0)固定点在惯性系中速度为 0block_point.set_vel(N, q1d * N.y)振子块只有竖直平动速度pendulum_point.v2pt_theory(block_point, N, B)利用同一刚体上两点速度关系v2pt_theory velocity of point 2 from point 1 theory计算摆球速度——结果为平移速度q1(t) n_y与旋转引起的切向速度-l*q2(t) b_x之和这正是摆球相对惯性系的实际速度符号上自动带出B系基向量。定义惯量与刚体摆被视为质量球 无质量轻杆的简单摆模型球质量为m杆长l铰接点固定在 Duffing 振子块上。两个刚体的惯量张量分别为 I_block M / 12 * me.inertia(N, h**2 d**2, w**2 d**2, w**2 h**2) I_pendulum 2*m*r**2/5*me.inertia(B, 1, 0, 1)I_block把振子块视为长方体长h、宽w、深d绕质心的惯量M/12*(h²d²)等项即长方体惯量公式(1/12)M(边长组合²)I_pendulum把摆球视为实心球半径r2mr²/5是实心球绕直径的转动惯量me.inertia(B, 1, 0, 1)表示沿B系三个主轴方向惯量分量分别为1、0、1绕B.y轴的惯量为 0 对应点质量理想化。注意惯量张量以摆坐标系B表达保证随摆转动。构建刚体时指定质心、连体参考系、质量与惯量 block_body me.RigidBody(block, block_point, N, M, (I_block, block_point)) pendulum_body me.RigidBody(pendulum, pendulum_point, B, m, (I_pendulum, pendulum_point))RigidBody(name, masscenter, frame, mass, (inertia, point))的第五个参数是一个二元组其中惯量张量可相对某个非质心点给出此处恰为质心。block_body的连体坐标系取惯性系Npendulum_body的连体坐标系取摆系B。定义力与载荷本例通过**作动器actuator**体系生成力载荷。先在O与block_point之间建立直线通路pathway path me.LinearPathway(O, block_point) spring me.DuffingSpring(k1, k2, path, 0) damper me.LinearDamper(c1, path)LinearPathway是连接一对附着点attachments的最简单通路两点之间沿直线其length表示两点间欧氏距离sqrt(q1**2)保证恒为正extension_velocity为其对时间的导数。通路决定了力的作用线与符号约定正力沿两点分开的方向expansileto_loads()会在两端点产生一对大小相等、方向相反的力。Duffing 弹簧的力-位移关系源码位于 sympy/physics/mechanics/actuator.pyF -linear_stiffness * displacement - nonlinear_stiffness * displacement**3 displacement pathway.length - equilibrium_length即弹簧力由线性项k1·x与立方非线性项k2·x³组成x为相对平衡长度 0 的伸长量本例equilibrium_length0。DuffingSpring的构造参数为(linear_stiffness, nonlinear_stiffness, pathway, equilibrium_length0)其中equilibrium_length为可选参数默认S.Zero所有系数经sympify(strictTrue)严格符号化。测试用例 test_actuator.py 验证了to_loads()产生的载荷恰为[Force(pA, F*N.x), Force(pB, -F*N.x)]。阻尼器LinearDamper(c1, path)产生与相对速度成比例的黏性阻尼力F -c1·extension_velocity。将两者载荷合并 loads spring.to_loads() damper.to_loads()然后为每个刚体添加重力保守力沿N.y方向g为重力加速度 bodies [block_body, pendulum_body] for body in bodies: ... loads.append(me.Force(body, body.mass * g * N.y))最终载荷列表文档输出为[(O, (k1*sqrt(q1²) k2*(q1²)^(3/2))*q1/sqrt(q1²) n_y), (block, -(k1*sqrt(q1²) k2*(q1²)^(3/2))*q1/sqrt(q1²) n_y), (O, c1*q1(t) n_y), (block, -c1*q1(t) n_y), (block, M*g n_y), (pendulum, g*m n_y)]可见弹簧与阻尼器在O与block两端成对出现牛顿第三定律重力则只作用在各自刚体的质心上。这里sqrt(q1²)保留了LinearPathway.length的长度恒正表达F化简后即-(k1·q1 k2·q1³)。拉格朗日方法求运动方程先计算系统的拉格朗日量L T - V本例势能已置零故只剩动能 L me.Lagrangian(N, block_body, pendulum_body) L M*q1(t)²/2 m*r²*q2(t)²/5 m*(l²*q2(t)² - 2*l*sin(q2)*q1(t)*q2(t) q1(t)²)/2三项分别对应振子块的平动动能Mq1²/2、摆球的旋转动能mr²q2²/5、以及摆球平动动能m/2·(l²q2² - 2l·sin(q2)·q1·q2 q1²)交叉项来自振子与摆的相对运动耦合。me.Lagrangian(frame, *bodies)会遍历各刚体RigidBody与Particle自动累加T mv²/2 ω·I·ω/2形式的动能。随后建立拉格朗日方法对象并求运动方程 LM me.LagrangesMethod(L, [q1, q2], bodiesbodies, forcelistloads, frameN) sm.simplify(LM.form_lagranges_equations()) [−M·g M·q1(t) c1·q1(t) − g·m − m·(l·sin(q2)·q2(t) l·cos(q2)·q2(t)² − q1(t)) (k1 k2·q1²)·q1] [m·(5·g·l·sin(q2) 5·l²·q2(t) − 5·l·sin(q2)·q1(t) 2·r²·q2(t))/5]LagrangesMethod(L, qs, ...)的bodies、forcelist、frame参数用于处理非保守力阻尼与显式加入的保守力重力、弹簧力。form_lagranges_equations()内部执行d/dt(∂L/∂q̇) − ∂L/∂q Q广义力得到两个耦合的二阶常微分方程。与论文方程的对照论文 [P.Brzeskia2012] 中的运动方程为(M m)y − ml·φ·sin(φ) − ml·(φ)²·cos(φ) k1·y k2·y³ c1·y F0·cos(νt) ml²·φ − ml·y·sin(φ) ml·g·sin(φ) c2·φ 0本示例得到的方程为以q1 ↔ y、q2 ↔ φ对照(M m)q1 − ml·q2·sin(q2) − ml·(q2)²·cos(q2) k1·q1 k2·q1³ c1·q1 − (M m)g 0 ml²·q2 − ml·q1·sin(q2) ml·g·sin(q2) (2r²/5)·q2 0两式差异来源清晰可辨重力项论文将重力并入等效外激励本例把(Mm)g作为常数重力项留在方程左侧第二式多了mlg·sin(q2)与论文一致第一式多出−(Mm)g外激励论文含周期性强迫力F0·cos(νt)受迫 Duffing 振子本例为自由振动未引入摆端阻尼论文第二式含阻尼力矩c2·φ本例未定义摆的铰点阻尼c2未出现在变量列表中摆惯量本例把摆球视为实心球显式带出转动惯量项(2r²/5)·q2论文采用无质量杆-质点摆理想化故无此项。进阶Duffing 弹簧作动器的扩展用法DuffingSpring位于 sympy/physics/mechanics/actuator.py继承自ForceActuator进而继承ActuatorBase类层次为ActuatorBase → ForceActuator → DuffingSpring。其全部属性在构造后不可修改setter 中通过hasattr(self, _attr)抛出AttributeError保证建模对象不可变、可安全重复用于不同求解器。参数速查参数类型含义默认值linear_stiffnessExpr线性刚度系数k1即论文中的beta必填nonlinear_stiffnessExpr立方项系数k2即论文中的alpha必填pathwayPathwayBase弹簧作用通路决定作用点与方向必填equilibrium_lengthExpr弹簧平衡长度位移按length − equilibrium_length计0to_loads()返回(Point, Vector)载荷元组列表与KanesMethod凯恩方法和LagrangesMethod拉格朗日方法的loads/forcelist参数格式兼容可以直接拼接其他载荷后交给任意求解方法。对pathway参数有严格类型检查必须是PathwayBase实例系数则通过sympify(..., strictTrue)统一符号化保证非符号输入也能被正确转换。若希望把弹簧的平衡位置设在某非零长度如预紧弹簧只需传入第 4 个位置参数例如me.DuffingSpring(k1, k2, path, l0)。该类的 repr 与to_loads行为均有对应单元测试覆盖见 test_actuator.py可作为自定义作动器如CoulombKineticFriction实现时的参考模板。运行与验证上述代码可在任意支持 SymPy 的 Python 环境中交互式执行如python -c ...或 Jupyter。注意文中为 doctest 提示符复制时可去掉需要 SymPy 版本中包含DuffingSpring该作动器属于较新 API请确认当前安装版本的sympy.physics.mechanics.actuator中存在该类本仓库的sympy.physics.mechanics源码与文档位于 sympy/physics/mechanics 与 doc/src/tutorials/physics/mechanics同目录下还有rollingdisc_example.rst、bicycle_example.rst、atwoods_machine_example.rst等更多建模示例可供对照学习。验证要点form_lagranges_equations()输出的第一式系数展开后与(Mm)q1 − ml·sin(q2)·q2 − ml·cos(q2)·q2² (k1k2q1²)q1 c1q1 − (Mm)g完全一致第二式经simplify后与ml²q2 − ml·sin(q2)·q1 mlg·sin(q2) (2r²/5)q2一致。你可以自行替换参数例如把equilibrium_length改为l0或增加摆端阻尼力矩来观察方程如何变化从而深化对拉格朗日建模各环节的理解。参考资料P. Brzeskia, P. Perlikowskia, S. Yanchukb, T. Kapitaniaka,The dynamics of the pendulum suspended on the forced Duffing oscillator, Journal of Sound and Vibration, 2012本示例方程对照的原始文献。duffing.svg本示例的系统结构示意图随文档发布的 SVG 矢量图。相关 API 文档可参阅 doc/src/modules/physics/mechanics 目录下的 mechanics 模块说明。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
