Моделирование подпружиненного маятника на вращающемся диске

Как только вы устраните ограничение маятника фиксированной длины, вероятно, будет проще решить эту проблему в декартовых координатах, используя 2-й закон Ньютона (F = ma), а не с помощью уравнений движения Лагранжа.

Рассмотрим пружину естественной длины L и жесткости k, один конец которой прикреплен к точке на краю диска радиуса R, вращающегося с угловой скоростью ω, а другой удерживает шарик массы m.

Точка привязки на диске имеет координаты, зависящие от времени.

У Боба есть координаты

На боб действуют две силы: сила упругости пружины величиной

и его вес

Уравнение движения просто

или

где n — единичный вектор направления троса, а g — ускорение свободного падения.

Обратите внимание, что здесь нет механизма рассеивания энергии, поэтому избегайте резонанса (т. е. собственная круговая частота пружины sqrt(k/m) не должна быть близка к круговой частоте омега точки вращения). .

Чтобы получить стабильную, а не очень «прыгучую» орбиту, см. основной код ниже.

Код:

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D

g = 9.81                                                             # acceleration due to gravity (m/s^2)
gravity = np.array( [ 0.0, 0.0, -g ] )

def plot_animation( t, qx, qy, qz, x, y, z ):
   plotInterval = 1
   fig = plt.figure( figsize=(4,4) )
   ax = fig.add_subplot(111, projection='3d' )
   ax.view_init(30, 45)
   ax.set_xlim( -2.0, 2.0 );   ax.set_ylim( -2.0, 2.0 );   ax.set_aspect('equal')

   ax.plot( qx, qy, qz, 'k' )                                        # tether point
   a = ax.plot( [qx[0],x[0]], [qy[0],y[0]], [qz[0],z[0]], 'g' )      # pendulum string
   b = ax.plot( [x[0]], [y[0]], [z[0]], 'ro' )                       # pendulum bob
   def animate( i ):                                                 # update anything that has changed
      a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] )   
      b[0].set_data_3d( [x[i]], [y[i]], [z[i]] )
   ani = animation.FuncAnimation( fig, animate, interval=4, frames=len( t ) )
   plt.show()
#  ani.save( "demo.gif", fps=50 )


def deriv( t, Y, R, omega, L, k, m ):
    x, v = Y[0:3], Y[3:6]                                                      # bob position and velocity (vectors)
    q = np.array( [ R*math.cos(omega*t), R*math.sin(omega*t), 0.0 ] )          # position vector of tether
    spring = q - x                                                             # vector from bob to tether
    length = np.linalg.norm( spring )                                          # stretched length of spring                                             
    unitVector = spring / length                                               # direction of elastic force
    extension = length - L                                                     # extension of spring
    a = ( k / m ) * extension * unitVector + gravity                           # F combines elastic force and gravity
    return np.hstack( ( v, a ) )                                               # return velocity and acceleration


R = 0.5            # radius of disk (m)
omega = 2.0        # angular velocity of disk (rad/s)
L = 1.0            # natural length of spring (m)
k = 20.0           # stiffness of spring (N/m)             # warning: avoid resonance!
m = 2.0            # mass (kg)
Y0 = np.array( [ R, 0.0, -L, 0.0, 0.0, 0.0 ] )                                 # initial position and velocity

period = 2 * np.pi / omega
tmax = 5 * period
solution = solve_ivp( deriv, [0, tmax], Y0, args=(R,omega,L,k,m), rtol=1.0e-6, dense_output=True )

t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )

# Position of bob
x = Y[0,:]
y = Y[1,:]
z = Y[2,:]

# Position of tether on disk
qx = R * np.cos( omega * t )
qy = R * np.sin( omega * t )
qz = np.zeros_like( qx )

plot_animation( t, qx, qy, qz, x, y, z )

Можно придать качанию маятника начальное положение и скорость, которые позволят ему перейти на устойчивую орбиту (постоянное движение по кругу, синхронно с тросом на диске). Если вы установите центростремительную силу, равную горизонтальной составляющей силы упругости, а вес, равный вертикальной составляющей силы упругости, то вы сможете одновременно решить пару уравнений для тангенса угла к вертикали и растяжения пружины и тем самым обеспечьте соответствующие начальные условия для Y0. Попробуйте заменить инициализацию Y0 следующими строками. Вероятно, существуют также синхронные решения, в которых боб отстает от точки привязки, но я их не исследовал.

# "Equilibrium motion" - moves in a circle, synchronised with the tether
t, e = 0, 0                  # tan(theta) and extension
w2g = omega ** 2 / g
mgk = m * g / k
for _ in range( 100 ):
     c = 1 / math.sqrt( 1 + t*t );   s = t * c        # cosine and sine
     t = w2g * ( R + (L+e) * s )
     e = mgk / c
     # print( t, e )
x = R + (L+e) * s
z =   - (L+e) * c
Y0 = np.array( [ x, 0.0, z, 0.0, x * omega, 0.0 ] )