Построение годографа единичного вектора продольной оси тела при движении в случае Лагранжа
Рассматривается движение твёрдого тела в случае Лагранжа. Необходимо построить годограф продольной оси тела (\(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]);