Определение момента сопротивления в шарнире физического маятника

Оценка постоянного момента сопротивления в оси вращения физического маятника по продолжительности его колебаний из некоторого начального положения до полной остановки. Численное интегрирование движения физического маятника в среде Python.

Постановка задачи

Рассмотрим движение физического маятника с известной массой \(m\), моментом инерции относительно оси вращения \(J\) и расстоянием от оси вращения до центра масс \(l\). Маятник совершает плоские колебания в вертикальной плоскости. В оси вращения маятника «О» действует малый постоянный момент сопротивления \(M\).

Предположим, что нам неизвестен момент сопротивления \(M\) и этот момент необходимо определить из эксперимента, оценив продолжительность затухания колебаний маятника при движении и заданного начального положения, определяемого углом \(\varphi_0\).

Уравнения движения

Запишем уравнение малых колебаний маятника: \[J \ddot{\varphi} = - m g l \varphi - M \text{sgn} (\varphi)\]

Рассмотрим движение маятника в течение половины периода из начального положения \(\varphi = \varphi_0 < 0\) при \(\dot \varphi > 0\). На этом интервале направление действия момента сопротивления не изменяется, поэтому уравнение движения можно записать в следующем виде: \[\ddot \varphi = -k^2 \varphi - f,\]

где \[k^2 = \frac{mgl}{J}, \quad f = \frac{M}{J}\]

Решение этого линейного неоднородного уравнения определяется как: \[\varphi = \left(\frac{f}{k^2} + \varphi_0\right) \cos ⁡kt - \frac{f}{k^2}\]

Зная период колебаний маятника: \[T = \frac{2 \pi}{k} = 2 \pi \sqrt{\frac{J}{mgl}},\]

определим угол поворота маятника через половину периода: \[\varphi(T/2) = \varphi(\pi/k) = - \varphi_0 - 2 \frac{f}{k^2}\]

Учитывая, что в начальный момент \(\varphi_0 < 0\), амплитуда колебаний маятника за половину периода уменьшится на \(\Delta A(T/2)= 2f/k^2\), а за весь период – на \(\Delta A(T) = 4f/k^2\). Средняя скорость уменьшения амплитуды колебаний равна: \[\frac{\Delta A(T)}{T} = \frac{4f}{k^2 T}\]

Теперь можно оценить продолжительность затухания колебаний маятника: \[t_d = \frac{\varphi_0}{\Delta A(T)} T = \frac{\varphi_0 k^2}{4f} T = \frac{\sqrt{mglJ}}{2M} \pi \varphi_0\]

Постоянный момент сопротивления \(M\) в шарнире при известной продолжительности затухания колебаний: \[M = \frac{\sqrt{mglJ}}{2t_d} \pi \varphi_0.\]

Численное интегрирование и “проверка” формулы

Численно проинтегрируем в среде Python уравнение движения нелинейного маятника (без допущения о малости его колебаний) для проверки полученной формулы.

Подключаемые библиотеки:

import numpy as np
from scipy.integrate import solve_ivp
from dataclasses import dataclass
import matplotlib.pyplot as plt

Объявляем гладкую функцию для аппроксимации функции модуля, которая в малой окрестности нуля (“малость” этого интервала зависит от параметра n) плавно изменяется от минус 1 до 1:

def sigmoid(x):
    if not hasattr(x, '__iter__'):
        x = (x,)
    n = 500
    return [np.sign(xi) if np.abs(xi) >= 1 else 2.0/(1.0+np.exp(-2*n*xi/(1-xi*xi)))-1 for xi in x]

Функция правых частей дифференциальных уравнений

def dqdt(t, q, p):  
    phi  = q[0]
    dphi = q[1]
    d2phi = - (p.m*p.g*p.l/p.J)*np.sin(phi) - p.M*np.sign(dphi)/p.J
    return np.hstack([dphi,d2phi])

Объявляем тип на основе dataclass для хранения параметров системы:

@dataclass
class Data:
    # Масса маятника
    m : float = 100;
    # Расстояние от точки подвески до центра масс (от оси вращения)
    l : float = 0.5;
    # Момент инерции маятника относительно оси вращения
    J : float =  100*0.5**2/3    
    # Момент сопротивления
    M : float = 3
    # Ускорение свободного падения
    g : float = 9.807

# Создаем экземпляр класса Data (объект)
p = Data()

Численно интегрируем уравнение движения маятника методом Рунге-Кутты с начальным углом отклонения маятника 20 градусов:

q0 = [np.deg2rad(20),0]
sol = solve_ivp(lambda t,q: dqdt(t,q,p), [0,13], q0, method='RK45', rtol = 1e-8)

Построим график изменения угла поворота маятника:

plt.plot(sol.t,np.rad2deg(sol.y[0]))
plt.grid(ls=':')
plt.xlabel('t, c')
plt.ylabel('Угол поворота, градус')

Из графика следует, что колебания затухают приблизительно через 11,5 секунд. Оценим момент трения по выведенной ранее формуле:

np.sqrt(p.m*p.g*p.l*p.J)*np.pi*q0[0]/(2*11.5)

> 3.047837910536908

Получилось 3.05 Нм. Эта оценка отличается от значения, заданного при интегрировании (3 Нм) менее чем на 2 %.

Литература

1. Гладков Сергей Октябринович, Богданова Софья Борисовна К вопросу учета силы сопротивления в шарнирной точке крепления физического маятника и ее влияние на динамику движения // Известия вузов. ПНД. 2019. №1.
2. Mungan C. E., Lipscombe T. C. Simple pendulum with speed-independent bearing friction //European Journal of Physics. – 2022. – Т. 43. – №. 4. – С. 045001. DOI 10.1088/1361-6404/ac646a.