Пошаговое руководство по использованию MATLAB для определения того, как космический корабль движется под действием гравитации большего тела.
Орбитальная механика (или астродинамика) включает в себя применение законов движения Ньютона и универсального закона тяготения к космическим кораблям (таким как ракеты и спутники). Он используется планировщиками миссий для прогнозирования движения космического корабля под действием гравитации, тяги и других сил. Темой данной статьи является задача двух тел. Это подкатегория астродинамики, в которой космический корабль находится под воздействием одного массивного тела без учета других сил. Масса космического корабля также считается незначительной по сравнению с большей массой. Это означает, что массивное тело будет влиять на движение космического корабля, но космический корабль не будет влиять на движение массивного тела. На самом деле интересующее тело не обязательно должно быть космическим кораблем; это может быть астероид, комета, астронавт и т. д.
Как решить задачу двух тел
Чтобы сгенерировать траекторию в системе двух тел, необходимо вывести уравнения движения космического корабля (массы интереса). Существует множество ресурсов, которые показывают этот вывод и окончательные уравнения движения, поэтому вам не нужно делать это самостоятельно. Я продемонстрировал этот процесс в статье ниже. Я призываю вас понять уравнения и то, как они были получены. Однако, если вы не хотите этого делать и просто хотите создать орбиту, я предоставил вам для справки диаграмму и необходимые уравнения движения ниже.
В этих уравнениях движения μ — гравитационный параметр первичной массы (планеты, луны, старта и т. д.). x, y и z — декартовы координаты, определяющие положение космического корабля в инерциальном пространстве. Точки представляют собой производные по времени, поэтому одна точка — это скорость, а две точки — ускорение. Мы будем использовать их в разделе кодирования статьи.
Между прочим, задача двух тел — это очень упрощенная версия сложной области орбитальной механики. В действительности при проектировании миссии приходится учитывать возмущающие силы, такие как атмосферное сопротивление, давление солнечного излучения, другие гравитационные силы и влияние несферических тел. Хотя это и не самый точный метод, задача двух тел может использоваться для начального планирования миссии космического корабля, особенно если он вращается близко к одному телу. Когда орбита находится далеко от массивного тела, другие силы будут играть большую роль в движении космического корабля. Дополнительные возмущающие силы здесь не обсуждаются, но о них полезно помнить.
Чтобы смоделировать орбиту, нам нужно будет использовать численное интегрирование. Если вы не знаете, что это такое, вы все равно сможете использовать код для создания своих собственных орбит. Чтобы численно интегрировать или решить уравнения движения космического корабля, мы будем использовать функцию ode113 в MATLAB. Это встроенная функция, которая принимает в качестве входных данных набор производных для интегрирования (определяемая пользователем функция), интересующий временной интервал, начальные условия и, возможно, допуск для численного интегрирования. Набор производных будет производным по времени от нашего вектора состояния Y, который включает компоненты положения и скорости космического корабля в векторной форме. И вектор состояния, и его производная по времени показаны ниже.
Это был краткий обзор процесса численного интегрирования уравнений движения задачи двух тел. Давайте запустим код, который будет моделировать траекторию спутника на орбите Земли!
Определяемая пользователем функция ODE
Как указывалось ранее, чтобы создать нашу траекторию, нам нужно будет численно интегрировать уравнения движения задачи двух тел. Для этого мы будем использовать ode113. Эта функция MATLAB требует, чтобы мы создали нашу собственную функцию, которая будет вызываться для выполнения численного интегрирования. Мы определим эту функцию, ODE2BP. Эта функция использует гравитационный параметр Земли mu и вектор состояния Y, полученные ранее, для создания производного по времени вектора состояния dYdt. Функция возвращает этот новый вектор, который будет использоваться ode113 для численного интегрирования. Примечание: в MATLAB определяемые пользователем функции должны располагаться в конце скрипта.
% User-Defined ODE Function
function dYdt = ODE2BP(t, Y)
mu = 3.986*10^5; % Earth's gravitational parameter [km^3/s^2]
x = Y(1); % [km]
y = Y(2); % [km]
z = Y(3); % [km]
vx = Y(4); % [km/s]
vy = Y(5); % [km/s]
vz = Y(6); % [km/s]
xddot = -mu/(x^2+y^2+z^2)^(3/2)*x; % [km/s^2]
yddot = -mu/(x^2+y^2+z^2)^(3/2)*y; % [km/s^2]
zddot = -mu/(x^2+y^2+z^2)^(3/2)*z; % [km/s^2]
dYdt = [vx;vy;vz;xddot;yddot;zddot]; % Y'
end
Численная интеграция
Далее мы заканчиваем определение входных данных для ode113. Одним из таких входных данных являются начальные условия для космического корабля. Это в виде вектора состояния из начала статьи. Обратите внимание, что компоненты положения и скорости по соглашению указаны в км и км/с. Далее у нас есть интересующий интервал времени в секундах. Наконец, есть много необязательных аргументов, которые вы можете передать ode113. Поскольку орбитальная механика требует точности, я обычно допускаю очень небольшой относительный допуск при численном интегрировании. Затем мы передаем эти входные данные с функцией, которую мы определили в предыдущем разделе, в решатель ОДУ. Функция выводит переменные t и Y , которые представляют собой временные шаги численного интегрирования и вектор состояния на каждом временном шаге. Затем мы можем извлечь историю времени для наших координат положения x, y и z из выходных данных.
% Creating Inputs for Numerical Integration
Y0 = [20000; 0; 0; 0; 2.9; 1.8]; % [x; y; z; vx; vy; vz] [km, km/s]
tspan = [0 24*60*60]; % One day [s]
options = odeset('RelTol', 1e-13); % Setting a tolerance
% Numerical Integration
[t, Y] = ode113(@ODE2BP, tspan, Y0, options);
% Pulling Position Data from Output
x = Y(:, 1); % [km]
y = Y(:, 2); % [km]
z = Y(:, 3); % [km]
Построение траектории
Наконец, мы можем использовать переменные x, y и z для создания траектории нашего космического корабля. Сначала мы инициализируем фигуру метками, сеткой и одинаковым набором осей (дает истинное представление орбиты в физическом пространстве). Затем мы можем использовать функцию sphere в MATLAB, чтобы создать представление Земли в начале нашего графика. Самое главное, траекторию можно построить с помощью функции plot3.
% Creating Figure
figure; hold on
title('Two-Body Trajectory', 'Interpreter', 'Latex')
xlabel('x', 'Interpreter', 'Latex')
ylabel('y', 'Interpreter', 'Latex')
zlabel('z', 'Interpreter', 'Latex')
axis equal
grid minor
view(30, 30)
% Creating/Plotting Spherical Earth
rm = 6378.14; % Radius of Earth [km]
[xEarth, yEarth, zEarth] = sphere(25);
surf(rm*xEarth,rm*yEarth,rm*zEarth, 'FaceColor', [0 0 1]);
% Plotting Trajectory
plot3(x, y, z, 'k')
hold off
Приведенный выше код создает следующий график:
На этом завершается процесс создания орбиты или траектории двух тел с помощью MATLAB. Поэкспериментируйте с начальными условиями и временным интервалом, чтобы создать свои собственные уникальные орбиты! Помните, что это чрезмерное упрощение большинства задач астродинамики. Если вы хотите получить лучшее представление о физике, вам нужно будет включить возмущающие силы. Это может быть довольно сложной задачей, поэтому лучше всего учиться на простой задаче двух тел и строить на ее основе. Если вы решите пойти дальше, удачи на вашем пути!
Если вы хотите создать анимацию, похожую на визуал в начале статьи, ознакомьтесь со статьей ниже. Спасибо, что нашли время прочитать. Если вы нашли это полезным, пожалуйста, подпишитесь на меня! Спасибо!
