Построение годографа единичного вектора продольной оси тела при движении в случае Лагранжа

Динамика твёрдого тела и систем тел

Рассматривается движение твёрдого тела в случае Лагранжа. Необходимо построить годограф продольной оси тела (\(Oz_2\)).

clear all;
% Параметры тела
m   = 4;    % Масса тела
Jxy = 0.5;  % Поперечный момент инерции
Jz  = 1;    % Продольный момент инерции
s   = 0.5;  % Расстояние от точки опоры до центра масс
g   = 9.807;% Ускорение свободного падения

% Начальные условия
theta0 = 0.7;   dtheta0 = 0.0;
psi0   = 0;     dpsi0   = 0.0;
phi0   = 0;     dphi0   = 10;

% Угловые скорости
wx      =   dtheta0;
wy      =   dpsi0*sin(theta0);
wz      =   dphi0+dpsi0*cos(theta0);

% Интегралы энергии и кинетического момента
E       = 0.5*(Jxy*(wx^2+wy^2)+Jz*wz^2+2*m*g*s*cos(theta0));
% Проекция веткора кинетического момента на направление вертикали
L       = Jxy*dpsi0*sin(theta0)^2+Jz*wz*cos(theta0);

Коэффициенты полинома относительно u (13)

coeff_dotu2  = [2*m*g*s/Jxy, -(2*E-Jz*wz^2)/Jxy-Jz^2*wz^2/Jxy^2,...
    -2*m*g*s/Jxy+2*Jz*wz*L/Jxy^2, (2*E-Jz*wz^2)/Jxy-L^2/Jxy^2];  

Корни полинома, отсортированные по возрастанию

u = sort(roots(coeff_dotu2));

\(k^2\) в выражении (16)

k2 = (u(2)-u(1))/(u(3)-u(1));

Коэффициент масштаба времени \(\tau\) (18)

tau_k = sqrt((u(3)-u(1))*m*g*s/Jxy);

Решения для углов и их производных. Угол нутации \(\theta\):

theta = @(t) acos(u(1)+(u(2)-u(1))*ellipj(tau_k*t,k2).^2);

Угловая скорость прецессии \(\dot \psi\) и угол прецессии \(\psi\):

dpsi  = @(t) (L-Jz*wz*cos(theta(t)))./(Jxy*(1-cos(theta(t)).^2));
psi   = @(t) integral(dpsi,0,t);

Угловая скорость собственного вращения \(\dot \varphi\)

dphi  = @(t) wz-dpsi(t)*cos(theta(t));

График изменения угла нутации

fplot(@(t) theta(t),[0, 1]);
xlabel('t, c'); 
ylabel('\theta, градус');

Формируем массив значений моментов времени для построения графика.

time_array = (0:0.01:3)';

При помощи функции arrayfun вызываем функцию psi(t) для каждого элемента массива time_array. Такой способ формирования массива значений угла \(\psi\) обусловлен тем, что функция integral, используемая для определения функции psi(t), не принимает только скалярные значение пределов интегрирования.

psi_array = arrayfun(psi,time_array);

Значения углов \(\theta\) может

theta_array = theta(time_array);

Координаты точек годографа единичного вектора, направленного по продольной оси тела:

godograph = [sin(theta_array).*sin(psi_array), ...
            -sin(theta_array).*cos(psi_array), ...
            cos(theta_array)];

Строим график:

plot3(godograph(:,1), godograph(:,2), godograph(:,3),'r-','LineWidth',2);

hold on;
[xs,ys,zs] = sphere(64);
sf = surf(xs,ys,zs); 
sf.EdgeColor = 'none';
sf.FaceColor =[0.8 0.8 0.8];
sf.FaceAlpha = 0.8;
hold off;

xlim([-1 1]); ylim([-1 1]); zlim([-1 1]);
box on; 
daspect([1 1 1]);


© 2023. All rights reserved.

Powered by Hydejack v9.1.6