Интегрирование уравнений движения механизма с двумя степенями свободы

Задание

1. Разработать программу для численного интегрирования уравнений движения системы с двумя степенями свободы для механизма из задания для курсовой работы по теоретической механики (2 курс).
2. Построить графики изменения обобщенных координат и скоростей (четыре графика на четырёх отдельных рисунках).
3. Построить одном рисунке графики изменения потенциальной, кинетической и полной энергии системы. Убедитесь, что полная энергия системы остается постоянной.

Пример

Схема механизма

Вывод уравнений

Вывод уравнений в Wolfram Mathematica

Первое уравнение: \[\ddot \varphi \left(m_2 \left(d^2-\left(l_0+x\right) \left(2 d-l_0-x\right)+4 h^2\right)+J_1\right)+2 m_2 \dot \varphi \dot x \left(-d+l_0+x\right)+2 h m_2 \ddot x = - g m_2 \left(\cos \varphi \left(-d+l_0+x\right)+2 h \sin \varphi \right)-g m_1 y_{c1} \sin \varphi\]

Второе уравнение: \[m_2 {\dot \varphi}^2 \left(d-l_0-x\right)+2 h m_2 \ddot \varphi +m_2 \ddot x = -c x-g m_2 \sin \varphi\]

где \(l_0\) – свободная длина пружины, \(c\) – жёсткость пружины, \(m_1\) – масса пластины, \(J_1\) – момент инерции пластины относительно оси О, \(m_2\) масса груза на пружине, \(d = AK/2\), \(y_{c1} = OC_1\).

Полученные уравнения движения можно записать в матричном виде: \[\begin{bmatrix} m_{11} & m_{12} \\ m_{21} & m_{22} \end{bmatrix} \begin{bmatrix} \ddot \varphi \\ \ddot x \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}\]

Элементы матрицы масс: \[\begin{align} & m_{11} = \left(m_2 \left(d^2-\left(l_0+x\right) \left(2 d-l_0-x\right)+4 h^2\right)+J_1\right) \\ & m_{12} = 2 h m_2 \\ & m_{21} = 2 h m_2 \\ & m_{22} = m_2 \end{align}\]

Слагаемые, не зависящие от ускорений (правые части): \[\begin{align} & b_1 = - g m_2 \left(\cos \varphi \left(-d+l_0+x\right)+2 h \sin \varphi \right)-g m_1 y_{c1} \sin \varphi - 2 m_2 \dot \varphi \dot x \left(-d+l_0+x\right) \\ & b_2 = - c x-g m_2 \sin \varphi - m_2 \left(\dot \varphi \right)^2 \left(d-l_0-x\right) \end{align}\]

MATLAB

Файл-функция правых частей dqdt.m

function dq  = dqdt(t, q, p)
    % Обобщенные координаты и скорости
    phi  = q(1);    x    = q(2);
    dphi = q(3);    dx   = q(4);

    % Элементы первой строки матрицы масс
    m11 = p.m2*(p.d^2-(p.L0+x)*(2*p.d-p.L0-x)+4*p.h^2)+p.J1;
    m12 = 2*p.h*p.m2;

    % Первая строка правой части
    b1  = -p.g*p.m2*(cos(phi)*(-p.d+p.L0+x)+2*p.h*sin(phi))-...
          p.g*p.m1*p.yc1*sin(phi) - 2*p.m2*dphi*dx*(-p.d+p.L0+x);

    % Элементы второй строки матрицы масс
    m21 = 2*p.h*p.m2;
    m22 = p.m2;
    
    % Вторая строка правой части
    b2  = -p.c*x - p.g*p.m2*sin(phi) - p.m2*(p.d-p.L0-x)*dphi^2;

    % Разрешаем систему относительно старших производных, решая 
    % систему линейных уравнений 

    % Матрица коэффициентов (симметричная)
    A   = [m11 m12; m21 m22];

    % Матрица правых частей
    B   = [b1;b2];

    % Решение системы линейных уравнений
    d2q = A\B;
    
    %  Обобщенные скорости и ускорения
    dq  = [dphi;dx;d2q];
end

Файл-скрипт для запуска процесса интегрирования main.m

clc;

p.m1  = 2;
p.m2  = 0.3;
p.L0  = 0.8;
p.c   = 5;
p.g   = 9.807;
p.R   = 1;
p.h   = p.R/2;
p.d   = sqrt(p.R^2 - p.h^2);
p.J1  = 1.258;
p.yc1 = 0.587;

phi0  = 0.3;
dphi0 = 1;
x0    = 0;
dx0   = 1;

[t, q] = ode113(@(t,q) dqdt(t,q,p), [0; 5], [phi0;x0;dphi0;dx0]);

figure;
plot(t,q(:,1),'LineWidth',2);
grid on;
xlabel('t, c'); ylabel('\phi, радиан');

figure;
plot(t,q(:,2),'LineWidth',2);
grid on;
xlabel('t, c'); ylabel('x, м');