Анимация и решение двойного движения маятника

В этом примере показано, как смоделировать движение двойного маятника при помощи MATLAB® и Symbolic Math Toolbox™.

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

Шаг 1: задайте смещение, скорость и ускорение двойных масс маятника

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

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

  • угловое положение первого боба θ1(t)

  • угловое положение второго боба θ2(t)

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

  • длина первого стержня L1

  • длина второго стержня L2

  • масса первого боба m1

  • масса второго боба m2

  • ускорение свободного падения g

Для простоты проигнорируйте массы двух твердых стержней. Задайте все переменные при помощи 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: определите уравнения движения

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

Во-первых, задайте силу первого стержня как T1, и сила второго стержня T2.

syms T_1 T_2

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

Оцените силы, действующие на 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 = 

-m1L1sin(θ1(t))t θ1(t)2-L1cos(θ1(t))2t2 θ1(t)=T2sin(θ2(t))-T1sin(θ1(t))- m_1* (L_1*sin (theta_1 (t)) * (diff (theta_1 (t), t)) ^2 - L_1*cos (theta_1 (t)) *diff (theta_1 (t), t, 2)) == T_2*sin (theta_2 (t)) - T_1*sin (theta_1 (t))

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

m1L1sin(θ1(t))2t2 θ1(t)+L1cos(θ1(t))t θ1(t)2=T1cos(θ1(t))-gm1-T2cos(θ2(t))m_1* (L_1*sin (theta_1 (t)) *diff (theta_1 (t), t, 2) + L_1*cos (theta_1 (t)) * (diff (theta_1 (t), t)) ^2) == T_1*cos (theta_1 (t)) - g*m_1 - T_2*cos (theta_2 (t))

Оцените силы, действующие на m2. Определите уравнения движения второго боба путем балансировки горизонтальных и вертикальных компонентов силы. Задайте эти два уравнения как символьные уравнения eqx_2 и eqy_2.

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

-m2L1sin(θ1(t))t θ1(t)2+L2sin(θ2(t))t θ2(t)2-L1cos(θ1(t))2t2 θ1(t)-L2cos(θ2(t))2t2 θ2(t)=-T2sin(θ2(t))- m_2* (L_1*sin (theta_1 (t)) * (diff (theta_1 (t), t)) ^2 + L_2*sin (theta_2 (t)) * (diff (theta_2 (t), t)) ^2 - L_1*cos (theta_1 (t)) *diff (theta_1 (t), t, 2) - L_2*cos (theta_2 (t)) *diff (theta_2 (t), t, 2)) ==-T_2*sin (theta_2 (t))

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

m2L1cos(θ1(t))t θ1(t)2+L2cos(θ2(t))t θ2(t)2+L1sin(θ1(t))2t2 θ1(t)+L2sin(θ2(t))2t2 θ2(t)=T2cos(θ2(t))-gm2m_2* (L_1*cos (theta_1 (t)) * (diff (theta_1 (t), t)) ^2 + L_2*cos (theta_2 (t)) * (diff (theta_2 (t), t)) ^2 + L_1*sin (theta_1 (t)) *diff (theta_1 (t), t, 2) + L_2*sin (theta_2 (t)) *diff (theta_2 (t), t, 2)) == T_2*cos (theta_2 (t)) - g*m_2

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

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

Уравнения движения имеют четыре неизвестные: θ1, θ2, T1, и T2. Оцените эти два неизвестные T1 и T2 от eqx_1 и eqy_1. Используйте solve функционируйте, чтобы найти T1 и T2.

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

Замените решениями T1 и 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: решите системные уравнения

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

Во-первых, задайте значения для масс в kg, длины стержня в m, и сила тяжести в m/s2 (Единицы СИ). Замените этими значениями в эти два приведенных уравнения.

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

cos(θ1(t))σ1-3sin(θ2(t))t θ2(t)22-sin(θ1(t))t θ1(t)2+3cos(θ2(t))2t2 θ2(t)2=-2sin(θ2(t))cos(θ1(t))2σ1+sin(θ1(t))2σ1+49sin(θ1(t))5cos(θ1(t))sin(θ2(t))-cos(θ2(t))sin(θ1(t))where  σ1=2t2 θ1(t)because(theta_1 (t)) *diff (theta_1 (t), t, 2) - (3*sin (theta_2 (t)) * (diff (theta_2 (t), t)) ^2)/2 - sin (theta_1 (t)) * (diff (theta_1 (t), t)) ^2 + (3*cos (theta_2 (t)) *diff (theta_2 (t), t, 2))/2 == - (2*sin (theta_2 (t)) * (cos (theta_1 (t)) ^2*diff (theta_1 (t), t, 2) + sin (theta_1 (t)) ^2*diff (theta_1 (t), t, 2) + (49*sin (theta_1 (t)))/5)) / (cos (theta_1 (t)) *sin (theta_2 (t)) - because(theta_2 (t)) *sin (theta_1 (t)))

eqn_2 = subs(eqRed_2)
eqn_2 = 

cos(θ1(t))t θ1(t)2+3cos(θ2(t))t θ2(t)22+sin(θ1(t))σ1+3sin(θ2(t))2t2 θ2(t)2=2cos(θ2(t))cos(θ1(t))2σ1+sin(θ1(t))2σ1+49sin(θ1(t))5cos(θ1(t))sin(θ2(t))-cos(θ2(t))sin(θ1(t))-495where  σ1=2t2 θ1(t)because(theta_1 (t)) * (diff (theta_1 (t), t)) ^2 + (3*cos (theta_2 (t)) * (diff (theta_2 (t), t)) ^2)/2 + sin (theta_1 (t)) *diff (theta_1 (t), t, 2) + (3*sin (theta_2 (t)) *diff (theta_2 (t), t, 2))/2 == (2*cos (theta_2 (t)) * (cos (theta_1 (t)) ^2*diff (theta_1 (t), t, 2) + sin (theta_1 (t)) ^2*diff (theta_1 (t), t, 2) + (49*sin (theta_1 (t)))/5)) / (cos (theta_1 (t)) *sin (theta_2 (t)) - because(theta_2 (t)) *sin (theta_1 (t))) - sym (49/5)

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

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

Элементы векторного V представляйте дифференциальные уравнения первого порядка, которые равны производной времени элементов S. Элементы S переменные состояния θ2, dθ2/dt, θ1, и dθ1/dt. Переменные состояния описывают угловые смещения и скорости двойного маятника.

S
S = 

(θ2Dtheta2θ1Dtheta1)[theta_2; Dtheta_2; theta_1; Dtheta_1]

Затем преобразуйте дифференциальные уравнения первого порядка в функцию 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));

Затем создайте объект Animation движения остановки первого боба маятника при помощи fanimator функция. По умолчанию, fanimator создает объект Animation с 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;

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

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 проигрывать анимацию двойного маятника.

Для просмотра документации необходимо авторизоваться на сайте