В этом примере показано, как смоделировать движение двойного маятника при помощи MATLAB® и Symbolic Math Toolbox™.
Решите уравнения движения двойного маятника и создайте анимацию, чтобы смоделировать двойное движение маятника.
Следующий рисунок показывает модель двойного маятника. Двойной маятник состоит из двух маятников, качается и два твердых стержня.
Опишите движение двойного маятника путем определения переменных состояния:
угловое положение первого боба
угловое положение второго боба
Опишите свойства двойного маятника путем определения переменных:
длина первого стержня
длина второго стержня
масса первого боба
масса второго боба
ускорение свободного падения
Для простоты проигнорируйте массы двух твердых стержней. Задайте все переменные при помощи 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);
Определите уравнения движения на основе законов Ньютона.
Во-первых, задайте силу первого стержня как , и сила второго стержня .
syms T_1 T_2
Затем создайте схемы свободного тела сил, которые действуют на обе массы.
Оцените силы, действующие на . Определите уравнения движения первого боба путем балансировки горизонтальных и вертикальных компонентов силы. Задайте эти два уравнения как символьные уравнения 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 =
eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g
eqy_1 =
Оцените силы, действующие на . Определите уравнения движения второго боба путем балансировки горизонтальных и вертикальных компонентов силы. Задайте эти два уравнения как символьные уравнения eqx_2
и eqy_2
.
eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t))
eqx_2 =
eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g
eqy_2 =
Четыре уравнения движения описывают кинематику двойного маятника. Оцените силы, действующие на стержни, и уменьшайте набор четырех уравнений к двум уравнениям.
Уравнения движения имеют четыре неизвестные: , , , и . Оцените эти два неизвестные и от eqx_1
и eqy_1
. Используйте solve
функционируйте, чтобы найти и .
Tension = solve([eqx_1 eqy_1],[T_1 T_2]);
Замените решениями и в 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]);
Эти два приведенных уравнения полностью описывают движение маятника.
Решите системные уравнения, чтобы описать движение маятника.
Во-первых, задайте значения для масс в , длины стержня в , и сила тяжести в (Единицы СИ). Замените этими значениями в эти два приведенных уравнения.
L_1 = 1; L_2 = 1.5; m_1 = 2; m_2 = 1; g = 9.8; eqn_1 = subs(eqRed_1)
eqn_1 =
eqn_2 = subs(eqRed_2)
eqn_2 =
Эти два уравнения являются нелинейными дифференциальными уравнениями второго порядка. Чтобы решить эти уравнения, преобразуйте их в дифференциальные уравнения первого порядка при помощи odeToVectorField
функция.
[V,S] = odeToVectorField(eqn_1,eqn_2);
Элементы векторного V
представляйте дифференциальные уравнения первого порядка, которые равны производной времени элементов S
. Элементы S
переменные состояния , , , и . Переменные состояния описывают угловые смещения и скорости двойного маятника.
S
S =
Затем преобразуйте дифференциальные уравнения первого порядка в функцию 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)')
Создайте анимацию колеблющегося двойного маятника.
Во-первых, создайте четыре функции то использование 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
проигрывать анимацию двойного маятника.