带有摆的杜芬振荡器

In this example we demonstrate the use of functionality provided in sympy.physics.mechanics for deriving the equations of motion for a system consisting of a Duffing oscillator with a pendulum. This example is inspired by the paper [P.Brzeskia2012] section 2.

系统将使用拉格朗日方程进行建模。\(M\) 是杜芬振子的质量，\(m\) 是摆的质量，\(l\) 是摆的长度。\(k_1\) 和 \(k_2\) 分别是弹簧刚度的线性和非线性部分，\(c_1\) 是杜芬振子的粘性阻尼系数。

>>> import sympy as sm
>>> import sympy.physics.mechanics as me
>>> me.init_vprinting()

定义变量

>>> 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)

  • \(h\): Duffing振荡器的高度

  • \(w\): Duffing振荡器的宽度

  • \(d\): Duffing振荡器的深度

  • \(r\): 摆锤的半径

  • \(q_1\): 广义坐标，表示Duffing振荡器的位置

  • \(q_2\): 广义坐标，表示摆的角度

定义运动学

定义所有参考框架和点。

>>> N = me.ReferenceFrame('N')
>>> B = N.orientnew('B', 'axis', (q2, N.z))

参考系中摆的角速度为：

>>> B.ang_vel_in(N)
q2'(t) n_z

Duffing 振荡器模块和摆的位置和速度是：

>>> O = me.Point('O')
>>> block_point = O.locatenew('block', q1 * N.y)
>>> pendulum_point = block_point.locatenew('pendulum', l * B.y)

O 是惯性参考系中的一个固定点。

>>> 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_x

定义惯性和刚体。这里，我们假设一个简单的单摆，它由一个质量为 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)

>>> 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))

定义力

我们计算作用在系统上的力。在这个例子中，我们在拉格朗日函数中将势能设为零，并将保守力（重力和达芬弹簧）包含在载荷中。

>>> path = me.LinearPathway(O, block_point)
>>> spring = me.DuffingSpring(k1, k2, path, 0)
>>> damper = me.LinearDamper(c1, path)

>>> loads = spring.to_loads() + damper.to_loads()

>>> bodies = [block_body, pendulum_body]

>>> for body in bodies:
...     loads.append(me.Force(body, body.mass * g * N.y))

>>> loads
            /      _____           3/2\                  /        _____           3/2\
      |     /   2       /  2\   |                  |       /   2       /  2\   |
      \k1*\/  q1   + k2*\q1 /   /*q1               \- k1*\/  q1   - k2*\q1 /   /*q1
 [(O, ------------------------------ n_y), (block, -------------------------------- n_y), (O, c1*q1'(t) n_y), (block, -c1*q1'(t) n_y), (block, M*g n_y), (pendulum, g*m n_y)]
                    _____                                         _____
                   /   2                                         /   2
                 \/  q1                                        \/  q1

拉格朗日方法

在问题设定下，可以计算拉格朗日量，并形成运动方程。

>>> L = me.Lagrangian(N, block_body, pendulum_body)
>>> L
         2      2       2     / 2       2                                     2\
 M*q1'(t)    m*r *q2'(t)    m*\l *q2'(t)  - 2*l*sin(q2)*q1'(t)*q2'(t) + q1'(t) /
 --------- + ------------ + ----------------------------------------------------
     2            5                                  2

>>> LM = me.LagrangesMethod(L, [q1, q2], bodies=bodies, forcelist=loads, frame=N)
>>> sm.simplify(LM.form_lagranges_equations())
 [                                       /                                    2          \   /          2\   ]
 [-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]
 [                                                                                                           ]
 [                     /                   2                                    2        \                   ]
 [                   m*\5*g*l*sin(q2) + 5*l *q2''(t) - 5*l*sin(q2)*q1''(t) + 2*r *q2''(t)/                   ]
 [                   ---------------------------------------------------------------------                   ]
 [                                                     5                                                     ]

运动方程在 [P.Brzeskia2012] 中：

\[ \begin{align}\begin{aligned}(M + m)y'' - ml \phi'' \sin(\phi) - ml(\phi')^2 \cos(\phi) + k_1 y + k_2 y^3 + c_1 y' = F_0 \cos(\nu t)\\ml^2 \phi'' - mly'' \sin(\phi) + mlg \sin(\phi) + c_2 \phi' = 0\end{aligned}\end{align} \]

本例中的运动方程：

\[ \begin{align}\begin{aligned}(M + m)q_1'' - mlq_2'' \sin(q_2) - ml(q_2')^2 \cos(q_2) + k_1 q_1 + k_2 q_1^3 + c_1 q_1' - (M + m)g = 0\\ml^2 q_2'' - mlq_1'' \sin(q_2) + mlg \sin(q_2) + \frac{2r^2q_2''}{5} = 0\end{aligned}\end{align} \]

运动方程的差异归因于几个因素：重力、由阻尼系数 \(c_2\) 表征的阻尼力矩，以及周期性变化的激励 \(F_0 \cos(\nu t)\)。

参考文献

[P.Brzeskia2012] (1,2)

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, https://doi.org/10.48550/arXiv.1202.5937