Пример построения и численного интегрирования уравнений движения механической системы с двумя степенями свободы.
# Символные математические преобразоания
import sympy as sp
# Для конструирования функций быстрых подстановок чисел в аналитические выражения
from sympy.utilities.lambdify import lambdify
# Для печати красивых формул
sp.init_printing()
# Массивы, матрицы
import numpy as np
# Для численного интегрирования дифференциальных уранвений
from scipy.integrate import solve_ivp
# Для графиков
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.pylab as pylab
# Настройка стиля графиков
pylab.rcParams.update({'legend.fontsize': 14,
'figure.figsize': (12, 8),
'axes.labelsize': 14,
'axes.titlesize':14,
'xtick.labelsize':14,
'ytick.labelsize':14,
'axes.grid' : True,
'grid.linestyle' : ':'})
Механическая система состоит из двух тел: однородной пластины и шарика, который рассматривается как материальная точка. Полукруглая однородная пластина массой $m_1$ шарнирно закреплена и может вращаться в вертикальной плоскости вокруг точки О (точка О находится в середине диаметра полукруглой пластины). Внутри пластины в канале АВ движется тяжёлый шарик массой $m_2$. Шарик соединён с пластиной пружиной с известной жёсткостью $c$ и свободной длиной $l_0$.
Таким образом, известны следующие параметры:
Схема механической системы приведена на следующем рисунке.
# Объявляем символьные параметры, которые будут использоваться при построении уравнений движения
l0, c, R, a, b, m1, m2, yc, g, J1, t = sp.symbols('l_0 c R a b m_1 m_2 y_c g J1 t')
# Дополнительно введем словарь замен для тех символов, которые можно выразить через другие.
# Выражать зависимые символы через независимые мы не будем, чтобы не усложнять вид уравнений
# Например, расстояние b выражается через a и R:
rules = {b: sp.sqrt(R**2-a**2)}
rules
# Обобщенные координаты: положение шарика и угол поворота пластины (функции времени)
x = sp.Function('x')
phi = sp.Function('phi')
# Вектор обобщенных координат
q = [ x(t), phi(t) ]
# скоростей (производная каждого элемента из q)
dq = [ sp.diff(var,t) for var in q ]
# и ускорений
d2q = [ sp.diff(var,t) for var in dq ]
Запишем матрицу поворота:
$$ \boldsymbol A_z = \begin{bmatrix} \cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{bmatrix}, $$или матрицу преобразования координат между неподвижной системой координат $Ox_0y_0$ и подвижной системой координат, связанной с пластиной, $Ox'y'$.
Например, в системе координат, связанной с пластиной, координатный столбец радиус-вектора точки А относительно оси вращения записывается следующим образом:
$$ \boldsymbol \rho_A' = \begin{bmatrix} -b \\ -a \end{bmatrix} $$Другими словами, точка А по отношению к точке О смещена в отрицательном направлении оси x' на величину $a$ и в отрицательном направлении оси y' на величину $b$.
Координаты точки А в неподвижной системе $Ox_0y_0$:
$$ \boldsymbol \rho_A = \begin{bmatrix} a \sin \varphi - b \cos \varphi \\ -a \cos \varphi - b \sin \varphi \end{bmatrix} $$Последнее выражение можно переписать, используя матрицу поворота $\boldsymbol A_z$:
$$ \boldsymbol \rho_A = \begin{bmatrix} \cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{bmatrix} \begin{bmatrix} -b \\ -a \end{bmatrix} = \boldsymbol A_z \boldsymbol \rho_A' $$Для любого координатного столбца $\boldsymbol \rho'$ вектора в подвижной системе, его координаты $\boldsymbol \rho$ в неподвижной системе определяются выражением
$$ \boldsymbol \rho = \boldsymbol A_z \boldsymbol \rho' $$Используя матрицу $\boldsymbol A_z$, можно записать координатый столбец шарика в неподвижной системе координат, зная координатный столбец его радиус-вектора в подвижной системе $Ox'y'$. Положение шарика относительно точки О определяется координатным столбцом вектора $\boldsymbol \rho_2'$, который представляет собой сумму двух векторов, имеющих простой вид:
$$ \boldsymbol \rho_2' = \begin{bmatrix} -b \\ -a \end{bmatrix} + \begin{bmatrix} x(t) \\ 0 \end{bmatrix} = \rho_A' + \begin{bmatrix} x(t) \\ 0 \end{bmatrix} $$Координатный столбец шарика $\boldsymbol r_2$ в неподвижной системе $Ox_0y_0$ определим, умножив его на матрицу поворота:
$$ \boldsymbol r_2 = A_z \rho_2' $$Не будем вручную раскрывать это выражение, хотя оно и несложное, доверив эту работу библиотеке sympy. Зная радиус-вектор шарика в неподвижной системе координат, можно найти его скорость
$$ \boldsymbol V_2 = \frac{d \boldsymbol r_2}{dt}, $$а выражение для скорости шарика будет необходимо при формировании выражения для кинетической энергии шарика $T_2$ и механической системы в целом $T$, чтобы построить уравнения Лагранжа:
$$ T_2 = \frac{m_2 |V_2|^2}{2} $$Дифференцирование координатного столбца радиус-вектора $\boldsymbol r_2$ тоже доверим библиотеке sympy.
# Символьная матрица поворота или матрица преобразования координат
Az = sp.Matrix([ [sp.cos(phi(t)), -sp.sin(phi(t))],
[sp.sin(phi(t)), sp.cos(phi(t))] ])
# Координатный столбец точки А в подвижной системе координат Ox'y'
rhoA = sp.Matrix([ [ -b ],
[ -a ] ])
rhoA
# Координатный столбец шарика относительно точки А в подвижной системе координат Ox'y'
# или координатный столбец радиус-вектора, проведённого из точки А в точку М
rho2 = sp.Matrix([ [ x(t) ],
[ 0 ] ])
rho2
# Координатный столбец радиус-вектора шарика относительно точки О в неподвижной системе координат
r2 = Az*(rho2 + rhoA)
r2
# Координатный столбец вектора скорости шарика относительно неподвижной системы координат
v2 = sp.diff(r2, t)
v2
Момент инерции пластины относительно оси вращения $O$:
$$ J_1 = \int r^2 \mu ds = \int_0^R \int_0^{\pi} r^2 \mu r d \gamma dr = \mu \int_0^R \int_0^{\pi} r^3 d \gamma dr $$где $\mu$ - масса единицы площади пластины. Для полукруглой пластины массой $m_2$ эта величина определяется отношением массы к площади половины круга:
$$ \mu = \frac{2 m_1}{\pi R^2} $$Интеграл $J_0$ легко найти "вручную":
$$ J_1 = \mu \int_0^R \int_0^{\pi} r^3 d \gamma dr = \frac{2 m_1}{\pi R^2} \pi \frac{1}{4} R^4 = \frac{m_1 R^2}{2} $$но можно воспользоваться возможностями sympy
r, gamma = sp.symbols('r, gamma')
(2*m1/(sp.pi*R**2))*sp.integrate(sp.integrate(r**3, (r,0,R)), (gamma,0,sp.pi))
# Запомним результат в словаре замен
rules[J1] = (2*m1/(sp.pi*R**2))*sp.integrate(sp.integrate(r**3, (r,0,R)), (gamma,0,sp.pi))
rules
Положение центра масс пластины:
$$ y_c = - \frac{1}{m_1} \int_0^R \int_0^{\pi} r \sin \gamma \mu r d \gamma dr $$rules[yc] = -(2/(sp.pi*R**2))*sp.integrate(sp.integrate(r**2*sp.sin(gamma), (r,0,R)), (gamma,0,sp.pi))
rules
# v2.transpose()*v2 -- скалярное произведение скорости на саму себя
# результат - квадрат абсолютной скорости скорости шарика
T2 = m2*(v2.transpose()*v2)[0]/2
T2
T1 = J1*sp.diff(phi(t),t)**2/2
T = (T1 + T2).simplify()
T
Потенциальная энергия пластины определяется высотой её центра масс над условно нулевым уровнем потенциальной энергии. Пусть плоскость нулевого уровня энергии проходит через ось вращения пластины.
$$ \Pi_1 = m_1 g y_c \cos \varphi $$P1 = m1*g*yc*sp.cos(phi(t))
Потенциальная энергия шарика:
$$ \Pi_2 = y_2 m_2 g $$P2 = r2[1]*m2*g
Потенциальная энергия пружины
$$ \Pi_s = \frac{c (x-l_0)^2}{2} $$PS = c*(x(t)-l0)*(x(t)-l0)/2
Потенциальная энергия системы:
P = P1 + P2 + PS
P.simplify()
L = T - P
sp.simplify(L)
# Функция построяния уранения Лагранжа по заданной функции Лагранжа L и обощенной коодинате q
Lagrange_equation = lambda L, q : sp.diff(sp.diff(L,sp.diff(q)),t) - sp.diff(L,q)
# Уравнения движения (без правой части = 0)
equations = [ sp.simplify(Lagrange_equation(L, q)) for q in (x(t), phi(t)) ]
equations
Разрешим уравнения относительно вторых производных
сauchy_equations = sp.solve(equations, d2q)
сauchy_equations
# Правая часть полученных уравнений (список выражений)
сauchy_equations_right_side = list(сauchy_equations.values())
сauchy_equations_right_side
# Список всех символов, встречающихся в уравнениях движения
symbol_list = [l0, c, R, a, b, m1, m2, yc, g, J1]
symbol_list
# Объединяем список символов со списком переменных состояния, включая время
params_and_state = [*symbol_list, t, *q, *dq]
params_and_state
#
# Словарь значений независимых параметров системы
#
params = {R: 1.0, a: 0.5, m1: 5, m2: 1.0, g: 9.807, c: 25, l0: 0.7}
# для каждого символа из symbol_list выполним подстановку из rules (символ может зависеть от других символов),
# а затем подставим значения из params
params_list = [sym.subs(rules).subs(params).evalf() for sym in symbol_list]
#
for el in zip(symbol_list, params_list):
print(el[0], '= {:4.2f}'.format(el[1]))
l_0 = 0.70 c = 25.00 R = 1.00 a = 0.50 b = 0.87 m_1 = 5.00 m_2 = 1.00 y_c = -0.42 g = 9.81 J1 = 2.50
# Функция для вычисления правых частей дифференциальных уравнений
# При помощи функции lambdify определим функцию,
# которая будет вычислять значения правых частей ДУ для заданных числовых параметров.
# Первый аргумент функции lambdify -- список символов,
# которые содержаться в вычисляемом выражении (сauchy_equations_right_side),
# второй аргумент -- само символьное выражение.
f_right_side = lambdify(params_and_state, сauchy_equations_right_side)
def dydt(t, y, params):
# t - время
# y[0] - x
# y[1] - phi
# y[2] - dx/dt
# y[3] - dphi/dt
# функция возвращает
# res[0] - dx/dt
# res[1] - dphi/dt
# res[2] - d2x/dt2
# res[3] - d2phi/dt2
# вычисляем вторые производные
d2xd2phi = f_right_side( *[*params, t, *y] )
res = y.copy()
# копируем скорости
res[0:2] = y[2:4]
# записываем ускорения
res[2:4] = d2xd2phi[0:2]
return res
# Начальные условия
q0 = [0.8, 0.5, 0, 0];
# Интегрирование
solution = solve_ivp(lambda t, y: dydt(t, y, params_list), [0, 10.0], q0,
rtol = 1e-6, method="LSODA", t_eval = np.linspace(0,10,500))
plt.plot(solution.t,solution.y[0]);
plt.xlabel('t, c');
plt.ylabel('x, м');
plt.plot(solution.t,solution.y[1]);
plt.xlabel('t, c');
plt.ylabel('$\\varphi$, радиан');
# Кинетическая, потенциальная и полная энергии системы
Energy = [T, P, T + P]
f_Energy = lambdify(params_and_state, Energy)
n_Energy = [ f_Energy(*[*params_list, state[0], *(state[1])]) for state in zip(solution.t, solution.y.transpose()) ]
plt.plot(solution.t,n_Energy);
plt.legend(['Кинетическая энергия','Потенциальная энергия','Полная энергия'])
plt.xlabel('t, c');
plt.ylabel('Энергия, Дж');
from matplotlib.patches import Circle, Wedge, Polygon
from matplotlib.collections import PatchCollection
import matplotlib.lines as mlines
from matplotlib.animation import FuncAnimation
# Функция для вычисления положения шарика по заданному символьному выражению (второй аргумент функции) и
# и фактическим значениям параметров системы и ее состоянию
n_r2 = lambdify(params_and_state, r2)
# Функция создает набор координат x,y для построения зигзагообразной линии, которая изображает пружину
def get_spring_line(length, coils, diameter):
x = np.linspace(0,length,coils*2)
y = [ diameter*0.5*(-1)**i for i in range(len(x))]
return np.array([x,y])
xy = get_spring_line(10, 10, 1)
plt.plot(xy[0],xy[1]);
# матрица поворота
nAz = lambda x: np.matrix([ [np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)] ])
# создаем рисунок
fig, ax = plt.subplots();
plt.xlim(-2,2)
plt.ylim(-1.5,0.5)
ax.set_aspect(1)
# полукруг (пластина) с центром в начале координат
body1 = Wedge([0,0], 1.0, -180+np.degrees(q0[1]), np.degrees(q0[1]), fc='#AACCAA')
# шарик
body2 = Wedge([0,0], 0.1, 0, 360, fc = 'b')
# пружина
spring_xy = get_spring_line(1, 10, 0.1)
spring = mlines.Line2D(spring_xy[0], spring_xy[1], lw = 0.5, color = 'r')
# Функция init добавляет к рисунку четыре объекта:
# пластину (полукруг), шарик (круг), пружину (ломаная линия) и еще один круг,
# изображающий ось вращения.
# Далее в функции get_frame СВОЙСТВА этих объектов будут изменятся (обычно это положение объекта)
def init():
# Пластина
ax.add_patch(body1)
# Ось вращения
ax.add_patch(Wedge([0,0], 0.03, 0, 360, ec = 'black', fc = 'white'))
# Шарик
ax.add_patch(body2)
# Пружина
ax.add_line(spring)
return spring, body1, body2,
def get_frame(i):
phi = solution.y[1, i]
x = solution.y[0, i]
t = solution.t[i]
# положение шарика относительно неподвижной системы координат
r2 = n_r2(*[*params_list, t , x, phi, 0, 0]).astype(np.float).flatten()
# положение точки А относительно неподвижной системы координат
rA = nAz(phi)*(rhoA.subs(rules).subs(params).evalf())
# Повророт пластины имитируем изменением углов, обозначающих начало и конец сектора
body1.set_theta1(-180 + np.degrees(phi))
body1.set_theta2(np.degrees(phi))
body1.get_path()
# Изменяем координаты центра шарика
body2.set_center([r2[0], r2[1]])
body2.get_path()
# Точки ломаной (пружины) заданы в системе координат, связанной с пластиной
# Используем матрицу преобразования координат для определения координат ломаной в
# неподвижной системе - в системе координат рисунка
spring_xy = nAz(phi)*get_spring_line(x-0.1, 10, 0.1)
# Подменяем спискок координат ломаной, которая еще в функции init была добавлена к рисунке
# Не забываем, что пружина начинается в точке А (добавляем координаты точки А)
spring.set_xdata(spring_xy[0] + rA[0])
spring.set_ydata(spring_xy[1] + rA[1])
return spring, body1, body2
anim = FuncAnimation(fig, get_frame, init_func=init, frames = len(solution.t), interval=50, blit=True)
anim.save('animation.mp4')