Шаг 1. Определение смещения, скорости и ускорения двойных маятниковых масс

На следующем рисунке показана модель двойного маятника. Двойной маятник состоит из двух маятниковых бобин и двух жестких стержней.

Опишите движение двойного маятника, определив переменные состояния:

Опишите свойства двойного маятника, определив переменные:

Для простоты игнорируйте массы двух жестких стержней. Укажите все переменные с помощью syms.

syms theta_1(t) theta_2(t) L_1 L_2 m_1 m_2 g

Определите перемещения двойного маятника в декартовых координатах.

x_1 = L_1*sin(theta_1);
y_1 = -L_1*cos(theta_1);
x_2 = x_1 + L_2*sin(theta_2);
y_2 = y_1 - L_2*cos(theta_2);

Определение скоростей путем дифференциации смещений по времени с помощью diff функция.

vx_1 = diff(x_1);
vy_1 = diff(y_1);
vx_2 = diff(x_2);
vy_2 = diff(y_2);

Найдите ускорения, дифференцируя скорости по времени.

ax_1 = diff(vx_1);
ay_1 = diff(vy_1);
ax_2 = diff(vx_2);
ay_2 = diff(vy_2);

Шаг 2: Определение уравнений движения

Определите уравнения движения на основе законов Ньютона.

Сначала задайте натяжение первого стержня, как ${\mathit{}}_{\mathrm{T1}}$, и натяжение второго стержня ${\mathit{}}_{\mathrm{T2}}$.

syms T_1 T_2

Затем создайте диаграммы свободного тела сил, которые действуют на обе массы.

Оцените силы, действующие на ${\mathit{}}_{\mathrm{m1}}$. Определите уравнения движения первого боба, уравновешивая горизонтальную и вертикальную составляющие силы. Укажите эти два уравнения как символические eqx_1 и eqy_1.

eqx_1 = m_1*ax_1(t) == -T_1*sin(theta_1(t)) + T_2*sin(theta_2(t))

eqx_1 = 
$$-m_1\hspace{0.17em}\left(L_1\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2-L_1\hspace{0.17em}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)\right)=T_2\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)-T_1\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)$$

eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g

eqy_1 = 
$$m_1\hspace{0.17em}\left(L_1\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)+L_1\hspace{0.17em}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2\right)=T_1\hspace{0.17em}\cos \left({\theta }_1\left(t\right)\right)-g\hspace{0.17em}{m}_1-T_2\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)$$

Оцените силы, действующие на ${\mathit{}}_{\mathrm{м2}}$. Определите уравнения движения второго боба, уравновешивая горизонтальную и вертикальную составляющие силы. Укажите эти два уравнения как символические eqx_2 и eqy_2.

eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t))

eqx_2 = 
$$-m_2\hspace{0.17em}\left(L_1\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2+L_2\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_2\left(t\right)\right)}^2-L_1\hspace{0.17em}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)-L_2\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_2\left(t\right)\right)=-T_2\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)$$

eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g

eqy_2 = 
$$m_2\hspace{0.17em}\left(L_1\hspace{0.17em}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2+L_2\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_2\left(t\right)\right)}^2+L_1\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)+L_2\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_2\left(t\right)\right)=T_2\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)-g\hspace{0.17em}{m}_2$$

Шаг 3: Оценка сил и уменьшение уравнений системы

Четыре уравнения движения описывают кинематику двойного маятника. Оцените силы, действующие на стержни, и уменьшите набор из четырех уравнений до двух уравнений.

Уравнение движения имеет четыре неизвестных: ${}_{\mathrm{start1}}$, ${}_{\mathrm{start2}}$, ${\mathit{}}_{\mathrm{T1}}$и ${\mathit{}}_{\mathrm{T2}}$. Оцените двух неизвестных ${\mathit{}}_{\mathrm{T1}}$ и ${\mathit{}}_{\mathrm{T2}}$ из eqx_1 и eqy_1. Использовать solve функция для поиска ${\mathit{}}_{\mathrm{T1}}$ и ${\mathit{}}_{\mathrm{T2}}$.

Tension = solve([eqx_1 eqy_1],[T_1 T_2]);

Заменить решения для ${\mathit{}}_{\mathrm{T1}}$ и ${\mathit{}}_{\mathrm{T2}}$ на eqx_2 и eqy_2.

eqRed_1 = subs(eqx_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);
eqRed_2 = subs(eqy_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);

Два приведенных уравнения полностью описывают движение маятника.

Шаг 4: Решение системных уравнений

Решите системные уравнения, чтобы описать движение маятника.

Сначала определите значения для масс в $\mathrm{кг}$, длины стержней в $\mathrm{м}$и силы тяжести в $\mathrm{}{\mathrm{}}^{\mathrm{м/с2}}$ (единицы СИ). Замените эти значения двумя приведенными уравнениями.

L_1 = 1;
L_2 = 1.5;
m_1 = 2;
m_2 = 1;
g = 9.8;
eqn_1 = subs(eqRed_1)

eqn_1 = 
$$\begin{array}{l}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\sigma }_1-\frac{3\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_2\left(t\right)\right)}^2}{2}-\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2+\frac{3\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_2\left(t\right)}{2}=-\frac{2\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\left({\cos \left({\theta }_1\left(t\right)\right)}^2\hspace{0.17em}{\sigma }_1+{\sin \left({\theta }_1\left(t\right)\right)}^2\hspace{0.17em}{\sigma }_1+\frac{49\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)}{5}\right)}{\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)-\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)}\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)\end{array}$$

eqn_2 = subs(eqRed_2)

eqn_2 = 
$$\begin{array}{l}\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_1\left(t\right)\right)}^2+\frac{3\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}{\left(\frac{\partial }{\partial t}\mathrm{ }{\theta }_2\left(t\right)\right)}^2}{2}+\sin \left({\theta }_1\left(t\right)\right)\hspace{0.17em}{\sigma }_1+\frac{3\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_2\left(t\right)}{2}=\frac{2\hspace{0.17em}\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\left({\cos \left({\theta }_1\left(t\right)\right)}^2\hspace{0.17em}{\sigma }_1+{\sin \left({\theta }_1\left(t\right)\right)}^2\hspace{0.17em}{\sigma }_1+\frac{49\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)}{5}\right)}{\cos \left({\theta }_1\left(t\right)\right)\hspace{0.17em}\sin \left({\theta }_2\left(t\right)\right)-\cos \left({\theta }_2\left(t\right)\right)\hspace{0.17em}\sin \left({\theta }_1\left(t\right)\right)}-\frac{49}{5}\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\frac{{\partial }^2}{\partial t^2}\mathrm{ }{\theta }_1\left(t\right)\end{array}$$

Эти два уравнения являются нелинейными дифференциальными уравнениями второго порядка. Чтобы решить эти уравнения, преобразуйте их в дифференциальные уравнения первого порядка с помощью odeToVectorField функция.

[V,S] = odeToVectorField(eqn_1,eqn_2);

Элементы вектора V представляют дифференциальные уравнения первого порядка, которые равны временной производной элементов S. Элементы S являются переменными состояния ${}_{\mathrm{start2}}$, $\mathit{}{}_{\mathrm{}}\mathrm{dü\; 2/dt}$, ${}_{\mathrm{start1}}$и $\mathit{}{}_{\mathrm{}}\mathrm{dstart1/dt}$. Переменные состояния описывают угловые смещения и скорости двойного маятника.

S

S = 
$$\left(\begin{array}{c}{\theta }_2\\ {\mathrm{Dtheta}}_2\\ {\theta }_1\\ {\mathrm{Dtheta}}_1\end{array}\right)$$

Затем преобразуйте уравнения первого порядка в функцию MATLAB с помощью маркера M.

M = matlabFunction(V,'vars',{'t','Y'});

Определите начальные условия переменных состояния как [pi/4 0 pi/6 0]. Используйте ode45 функция для решения переменных состояния. Решения являются функцией времени в интервале [0 10].

initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);

Постройте график решений переменных состояния.

plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')

Шаг 5: Создание анимации колеблющегося двойного маятника

Создайте анимацию колебательного двойного маятника.

Сначала создайте четыре функции, которые используют deval для оценки координат обоих маятников из решений sols.

x_1 = @(t) L_1*sin(deval(sols,t,3));
y_1 = @(t) -L_1*cos(deval(sols,t,3));
x_2 = @(t) L_1*sin(deval(sols,t,3))+L_2*sin(deval(sols,t,1));
y_2 = @(t) -L_1*cos(deval(sols,t,3))-L_2*cos(deval(sols,t,1));

Затем создайте объект анимации стоп-движения первого маятникового узла с помощью fanimator функция. По умолчанию fanimator создает анимационный объект с 10 сгенерированными кадрами в единицу времени в диапазоне t от 0 до 10. Постройте график координат с помощью plot функция. Задайте для осей X и Y одинаковую длину.

fanimator(@(t) plot(x_1(t),y_1(t),'ro','MarkerSize',m_1*10,'MarkerFaceColor','r'));
axis equal;

Затем добавьте объекты анимации первого жесткого стержня, второго маятника и второго жесткого стержня.

hold on;
fanimator(@(t) plot([0 x_1(t)],[0 y_1(t)],'r-'));
fanimator(@(t) plot(x_2(t),y_2(t),'go','MarkerSize',m_2*10,'MarkerFaceColor','g'));
fanimator(@(t) plot([x_1(t) x_2(t)],[y_1(t) y_2(t)],'g-'));

Добавление фрагмента текста для подсчета прошедшего времени с помощью text функция. Использовать num2str для преобразования параметра времени в строку.

fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)));
hold off;

Используйте команду playAnimation для воспроизведения анимации двойного маятника.