Changing Orbits: Hohmann Transfer

Satellites often manoeuvre to different orbits for various reasons and, while simple, most coplanar manoeuvres can still be modelled with Hohmann Transfer, where the satellite changes semimajor axis and eccentricity, for example correcting post-launch injection errors, executing a Collision Avoidance Manoeuvre or moving to the mission orbit.

For this example, imagine a satellite on an elliptical orbit after the injection. The aim will be to increase the altitude and circularise the orbit. As the orbits are not much different from rails that trains move on, if two orbits do not intersect, it is mandatory to use an intermediate orbit that connects the initial orbit, and the target orbit. This intermediate or transfer orbit will start from the periapsis of the initial orbit and will have its apoapsis on the target orbit. There will be instantaneous velocity changes (or thruster firings) at these intersection points to jump from one orbit to the other one.

It is possible to compute the result analytically, without using the tools offered by SatMAD, but we will also demonstrate how we can put together a propagator, add an instantaneous velocity change (like a control input, as computed via analytical means) and use the result as the initial condition for the new propagation, putting together a complete mission profile. This procedure, computing the correct time of firing and the correct velocity change (or thrust vector) direction is the key to orbit control. The steps given in this example are thus identical to real orbit change operations.

For our example, we will start with an elliptic orbit at a semimajor axis of 7056 km, and we will end up at a circular orbit at a semimajor axis of 7500 km.

We start by defining the initial orbit:

import numpy as np
from astropy.time import Time
from astropy import units as u

from satmad.coordinates.frames import init_pvt
from satmad.core.celestial_bodies_lib import EARTH
from satmad.propagation.classical_orb_elems import OsculatingKeplerianOrbElems
from satmad.propagation.numerical_propagators import NumericalPropagator
from satmad.utils.discrete_time_events import DiscreteTimeEvents
from satmad.utils.timeinterval import TimeInterval

time = Time("2020-01-11T11:00:00.000", scale="utc")
central_body = EARTH

sm_axis = 7056.0 * u.km
ecc = 0.02 * u.dimensionless_unscaled
incl = 0 * u.deg
raan = 0 * u.deg
arg_perigee = 90 * u.deg
true_an = 20 * u.deg

init_orb_elems = OsculatingKeplerianOrbElems(
    time, sm_axis, ecc, incl, raan, arg_perigee, true_an, central_body
)

The first step is to find the time of first periapsis where the firing will be executed. For this, we will propagate the orbit for one period and look for the minimum radius value (which is the apoapsis point).

# generate cartesian initial conditions
pvt0 = init_orb_elems.to_cartesian()

# Set up propagation config
stepsize = 10 * u.s

prop_start = pvt0.obstime
prop_duration = init_orb_elems.period

# init propagator with defaults - run propagation and get trajectory
prop = NumericalPropagator(stepsize)
trajectory_init = prop.gen_trajectory(pvt0, TimeInterval(prop_start, prop_duration))

# Extract search range
time_list = trajectory_init.coord_list.obstime
r_list = trajectory_init.coord_list.cartesian.without_differentials().norm()

# Find time events
events = DiscreteTimeEvents(time_list, r_list)

# Min / Max Event times
print(events.max_min_table)

# Cross-check with orbital elems
print("\nCheck with initial conditions:")
print(f"Apoapsis r : {init_orb_elems.apoapsis}")
print(f"Periapsis r: {init_orb_elems.periapsis}")

          time          type       value
                                     km
----------------------- ---- ------------------
2020-01-11T11:43:54.264  max  7197.120000000502
2020-01-11T12:33:03.563  min 6914.8800000043475

Check with initial conditions:
Apoapsis r : 7197.12 km
Periapsis r: 6914.88 km

As can be seen, we have found the times corresponding to minimum and maximum radii, and the values match the periapsis and apoapsis of the initial orbit.

We can generate the transfer orbit analytically, using the properties of the initial and target orbits:

# semimajor axis and eccentricity of the transfer orbit
a_tgt = 7500 * u.km
a_transfer = 0.5*(init_orb_elems.periapsis + a_tgt)
ecc_transfer = (a_tgt - init_orb_elems.periapsis) /(2*a_transfer)

print("Transfer orbit properties:")
print(f"Transfer sma : {a_transfer}")
print(f"Transfer ecc : {ecc_transfer}")

Transfer orbit properties:
Transfer sma : 7207.4400000000005 km
Transfer ecc : 0.040591388898138576

The problem is to calculate how much velocity change is required to jump from the initial orbit periapsis to the transfer orbit (\(\Delta V_1\)), and then how much velocity change is required to jump from the transfer orbit apoapsis to the target orbit (\(\Delta V_2\)).

We will use the classical analytical method and go through the energy computed at the intersection points. The energy equation in terms of cartesian coordinates and semimajor axis is given as:

\(\varepsilon = \dfrac{v^2}{2} - \dfrac{\mu}{r} = - \dfrac{\mu}{2 a}\)

We can evaluate the equation at the intersection point of initial orbit and transfer orbit and solve for the velocity - this is essentially the velocity after the first thruster firing (\(v_{init@peri} + \Delta V_1\)). Note that the applied \(\Delta V\) is along the initial velocity vector.

\(\varepsilon_{tr} = \dfrac{v_{tr@peri}^2}{2} - \dfrac{\mu}{r_{init@peri}} = - \dfrac{\mu}{2 a_{tr}}\)

\(v_{tr@peri}^2 = \dfrac{2\mu}{r_{init@peri}} - \dfrac{\mu}{a_{tr}}\)

\(v_{tr@peri} = |\vec{v}_{init@peri}| + \Delta v_1\)

# Analytical method: velocity of the transfer orbit at periapsis
v_tr_peri_comp = np.sqrt(EARTH.mu * (2/init_orb_elems.periapsis - 1/a_transfer))

# Using satmad: find the next periapsis
# velocity of the init orbit at periapsis
t_init_peri = events.max_min_table[1]["time"]
rv_init_peri = trajectory_init(t_init_peri)

# required delta_v_1
delta_v_1_comp = v_tr_peri_comp - rv_init_peri.velocity.norm()
print("Compare velocities at init orbit periapsis and transfer orbit periapsis")
print(f"Velocity init orbit @ periapsis     : {rv_init_peri.velocity.norm()}")
print(f"Velocity transfer orbit @ periapsis : {v_tr_peri_comp}")
print(f"Required delta V                    : {delta_v_1_comp.to(u.m/u.s)}")

Compare velocities at init orbit periapsis and transfer orbit periapsis
Velocity init orbit @ periapsis     : 7.667903697041672 km / s
Velocity transfer orbit @ periapsis : 7.744915393079457 km / s
Required delta V                    : 77.01169603778446 m / s

This shows that a \(\Delta V\) of 77 m/s is required to make the jump from the initial orbit to the transfer orbit. As a cross-check, we will apply this \(\Delta V\) to the initial orbit at periapsis and see whether we can actually get the transfer orbit properties.

# Unit velocity and delta V vectors
v_unit_1 = rv_init_peri.velocity / rv_init_peri.velocity.norm()
delta_v_1 = delta_v_1_comp * v_unit_1

# Generate the coordinates of the transfer orbit
v_tr_peri = rv_init_peri.velocity + delta_v_1
r_tr_peri = rv_init_peri.cartesian.without_differentials()

rv_tr_peri = init_pvt("gcrs", rv_init_peri.obstime, r_tr_peri, v_tr_peri)

# Compute orbital elements of the transfer orbit
tr_orb_elems = OsculatingKeplerianOrbElems.from_cartesian(rv_tr_peri)

print("Computed Transfer orbit properties:")
print(f"Transfer sma : {tr_orb_elems.sm_axis}")
print(f"Transfer ecc : {tr_orb_elems.eccentricity}")

Computed Transfer orbit properties:
Transfer sma : 7207.439999999193 km
Transfer ecc : 0.040591388898082684

This confirms that the additional \(\Delta V\) really changes the initial orbit to the transfer orbit.

Before we go on to compute \(\Delta V_2\), we have to obtain the velocity at the apoapsis of the transfer orbit. This is where we are touching the target orbit but, because of the lack of kinetic energy, we cannot jump to the target orbit. We can either compute the velocity analytically from the energy at the transfer orbit evaluated at the apoapsis, or, as above, we can run a propagation, compute the apoapsis location and evaluate the velocity there. We have already seen that the results are the same.

\(\varepsilon_{tr} = \dfrac{v_{tr@apo}^2}{2} - \dfrac{\mu}{r_{tr@apo}} = - \dfrac{\mu}{2 a_{tr}}\)

\(\dfrac{v_{tr@apo}^2}{2} = \dfrac{\mu}{r_{tr@apo}} - \dfrac{\mu}{2 a_{tr}}\)

\(v_{tr@apo}^2 = \dfrac{2\mu}{r_{tr@apo}} - \dfrac{\mu}{a_{tr}}\)

# analytical method: velocity of the transfer orbit at periapsis
v_tr_apo_comp = np.sqrt(EARTH.mu * (2/tr_orb_elems.apoapsis - 1/a_transfer))

# Using satmad: propagate the orbit, find the apoapsis
# Set up propagation config and run the propagation
prop_start = rv_tr_peri.obstime
prop_duration = tr_orb_elems.period

trajectory_tr = prop.gen_trajectory(rv_tr_peri, TimeInterval(prop_start, prop_duration))

# Extract search range
time_list = trajectory_tr.coord_list.obstime
r_list = trajectory_tr.coord_list.cartesian.without_differentials().norm()

# Find time events
events = DiscreteTimeEvents(time_list, r_list)

# Min / Max Event times
print(events.max_min_table)

          time          type       value
                                     km
----------------------- ---- -----------------
2020-01-11T12:33:03.563  min 6914.879999999618
2020-01-11T13:23:48.319  max 7500.000000000302
2020-01-11T14:14:33.076  min  6914.87999996385

This shows where the next apoapsis will be, at the target radius of 7500 km.

Now that we are on the transfer orbit, we can calculate the \(\Delta V_2\), corresponding to the second manoeuvre, to jump from the transfer orbit apoapsis to the target orbit. We will use the same technique as above, computing the energy at the transfer orbit and getting the required additional \(\Delta V\). Note that, at this stage \(r_{tr@apo}\) is equal to \(a_{tgt}\).

\(\varepsilon_{tgt} = \dfrac{v_{tgt}^2}{2} - \dfrac{\mu}{r_{tr@apo}} = - \dfrac{\mu}{2 a_{tgt}}\)

\(\dfrac{v_{tgt}^2}{2} = \dfrac{\mu}{r_{tr@apo}} - \dfrac{\mu}{2 a_{tgt}} = \dfrac{\mu}{a_{tgt}} - \dfrac{\mu}{2 a_{tgt}}\)

\(v_{tgt}^2 = \dfrac{\mu}{a_{tgt}}\)

\(v_{tgt} = v_{tr@apo} + \Delta V_2\)

Finally, the total \(\Delta V\) is the sum of \(\Delta V_1\) and \(\Delta V_2\) values.

# velocity of the target orbit
v_tgt = np.sqrt(EARTH.mu / a_tgt)

# required delta_v_2
delta_v_2_comp = v_tgt - v_tr_apo_comp

print("Compare velocities at  transfer orbit apoapsis and target orbit")
print(f"Velocity transfer orbit @ apoapsis : {v_tr_apo_comp}")
print(f"Velocity target orbit              : {v_tgt}")
print(f"Required delta V                   : {delta_v_2_comp.to(u.m/u.s)}")
print("\n--------------------------------------------------------------")
print(f"Required total delta V             : {(delta_v_1_comp +delta_v_2_comp).to(u.m/u.s)}")

Compare velocities at  transfer orbit apoapsis and target orbit
Velocity transfer orbit @ apoapsis : 7.140688073774203 km / s
Velocity target orbit              : 7.290180078251383 km / s
Required delta V                   : 149.49200447717993 m / s

--------------------------------------------------------------
Required total delta V             : 226.50370051496438 m / s

The final step is to extract the coordinates of the transfer orbit at its apoapsis and add \(\Delta V_2\) to find the final orbit numerically - confirming that the numerical method agrees with the analytical method.

# coordinates of the transfer orbit at apoapsis
t_tr_apo = events.max_min_table[1]["time"]
rv_tr_apo = trajectory_tr(t_tr_apo)

# Unit velocity and delta V vectors
v_unit_2 = rv_tr_apo.velocity / rv_tr_apo.velocity.norm()
delta_v_2 = delta_v_2_comp * v_unit_2

# Generate the coordinates of the transfer orbit
v_tgt = rv_tr_apo.velocity + delta_v_2
r_tgt = rv_tr_apo.cartesian.without_differentials()

rv_tgt = init_pvt("gcrs", rv_tr_apo.obstime, r_tgt, v_tgt)

# Compute orbital elements of the target orbit
tgt_orb_elems = OsculatingKeplerianOrbElems.from_cartesian(rv_tgt)

print("Computed Target orbit properties:")
print(f"Target sma : {tgt_orb_elems.sm_axis}")
print(f"Target ecc : {tgt_orb_elems.eccentricity}")

Computed Target orbit properties:
Target sma : 7500.000000004719 km
Target ecc : 6.437006045365392e-11

This confirms that we were able to reach the target orbit via adding \(\Delta V\) vectors to the velocities at the correct times executing a Hohmann Transfer numerically.