Simulate Hohmann Transfer in Satellite Scenario using Impulsive Maneuvers

This example shows how to model a Hohmann transfer between two circular coplanar orbits using impulsive maneuvers, calculate fuel requirements, and visualize the transfer in a satellite scenario viewer.

Hohmann Transfer Overview

A Hohmann transfer uses two coplanar burns to move a spacecraft between orbits on the same plane. In this example, both orbits are circular. The transfer orbit is an ellipse tangent to the initial and final orbits:

• Injection burn at the initial orbit tangent point places the spacecraft on the transfer ellipse.

• Circularization burn at the final orbit tangent point inserts the spacecraft into the target orbit.

The transfer angle is 180°. Delta‑V at each burn equals the difference between circular speed and transfer orbit speed at that point, directed tangentially. The time between burns is half the transfer orbit period. Burns are modeled as impulsive for simplicity and efficiency during early mission design.

If the circularization burn is not performed, the spacecraft follows the dashed purple trajectory back to the first tangent point.

Suppose $a_1$ and $a_2$ are the orbital radii of the initial and final circular orbits. The semimajor axis of the transfer orbit is calculated using this equation:

$$a_T=\frac{\left(a_1+a_2\right)}{2}$$

The velocity magnitude at any given point on the initial circular orbit is:

$$v_1=\sqrt{\frac{\mu }{a_1}}$$

where $$\mu$$ is the standard gravitational parameter of the planet (Earth). Immediately following the injection burn, the spacecraft is at the periapsis of the transfer orbit. The corresponding velocity is derived from this energy equation:

$$\frac{v_{T,1}^2}{2}-\frac{\mu }{a_1}=-\frac{\mu }{2a_T}$$

where $$v_{T,1}$$ is the velocity on the transfer orbit at its periapsis.

$$v_{T,1}=\sqrt{2\mu \left(\frac{1}{a_1}-\frac{1}{2a_T}\right)}$$

The delta-v for the injection burn ($\Delta v_1$) is calculated as:

$$\Delta v_1=v_{T,1}-v_1$$

The direction of $\Delta v_1$ is the tangent at the periapsis of the transfer orbit. Similarly, the delta-v for the circularization burn ($\Delta v_2$) can be calculated using this equation:

$$\begin{array}{l}{v}_2=\sqrt{\frac{\mu }{a_2}}\\ v_{T,2}=\sqrt{2\mu \left(\frac{1}{a_2}-\frac{1}{2a_T}\right)}\\ \Delta v_2=v_2-v_{T,2}\end{array}$$

The direction of $\Delta v_2$ is the tangent at the apoapsis of the transfer orbit, as illustrated above. The time between the burns is one-half the period of the transfer orbit.

$$t_T=\pi \sqrt{\frac{a_T^3}{\mu }}$$

Assume each burn achieves the required velocity change instantly, modeling them as impulsive maneuvers. Use impulsive maneuvers during early mission design for simplicity, analytical solutions in cases like Hohmann transfers, and lower computation cost. They approximate finite burns well when burn duration is short compared to the overall trajectory, which is typical for chemical propulsion missions.

Mission Overview

The mission begins on December 14, 2022, 1:04 AM UTC.

startTime = datetime( ...
    2022, ... % year
    12, ...   % month
    14, ...   % day
    1, ...    % hour
    4, ...    % minute
    0);       % second

Initial orbit parameters are set to:

• Semimajor axis = 7,000 km

• Eccentricity = 0

• Inclination = 45 degrees

• Right ascension of ascending node = 90 degrees

• Argument of periapsis = 30 degrees

• True anomaly = 30 degrees

semiMajorAxisInitial = 7000e3;      % m
eccentricity = 0;
inclination = 45;                   % deg
rightAscensionOfAscendingNode = 90; % deg
argumentofPeriapsis = 30;           % deg
trueAnomaly = 30;                   % deg

Final orbit semimajor axis is 10,000 km.

The spacecraft completes one revolution on the initial orbit, transfers to the final orbit, then completes one revolution there.

semiMajorAxisFinal = 10000e3; % m

Calculate Injection Burn Parameters

Define standard gravitational parameter for Earth.

standardGravitationalParameter = 398600.4418e9; % m^3/s^2

Calculate the circular speed on the initial orbit.

vCircularInitial = sqrt(standardGravitationalParameter/semiMajorAxisInitial);

Calculate the semimajor axis of the transfer orbit.

semiMajorAxisTransfer = (semiMajorAxisInitial + semiMajorAxisFinal)/2;

Calculate the speed at periapsis of the transfer orbit.

vPeriapsisTransfer = sqrt(2*standardGravitationalParameter* ...
    ((1/semiMajorAxisInitial) - (1/(2*semiMajorAxisTransfer))));

Calculate the delta‑V vector in the VNB frame, where the x‑axis aligns with the velocity vector, the y‑axis aligns with orbital angular momentum, and the z‑axis completes the triad. Because the delta‑V direction is tangential at the transfer orbit periapsis, which is also tangential to the initial orbit, the y and z components are zero.

injectionBurnDeltaVMagnitude = vPeriapsisTransfer - vCircularInitial;
injectionBurnDeltaV = [injectionBurnDeltaVMagnitude;0;0];
injectionBurnReferenceFrame = "vnb";

Schedule the maneuver one orbit period after the mission start because the spacecraft must complete one full revolution on the initial orbit.

initialOrbitPeriod = ...
    2*pi*sqrt((semiMajorAxisInitial^3)/standardGravitationalParameter);
injectionBurnTime = startTime + seconds(initialOrbitPeriod);

Specify a name for the maneuver.

injectionBurnName = "Injection burn";

Calculate Circularization Maneuver Parameters

Calculate the speed at apoapsis of the transfer orbit.

vApoapsisTransfer = sqrt(2*standardGravitationalParameter* ...
    ((1/semiMajorAxisFinal) - (1/(2*semiMajorAxisTransfer))));

Calculate the circular speed on the final orbit.

vCircularFinal = sqrt(standardGravitationalParameter/semiMajorAxisFinal);

Calculate the delta-V vector for the circularization burn in the VNB frame.

circularizationBurnDeltaVMagnitude = vCircularFinal - vApoapsisTransfer;
circularizationBurnDeltaV = [circularizationBurnDeltaVMagnitude;0;0];
circularizationBurnReferenceFrame = "vnb";

The circularization burn is executed at the apoapsis of the transfer orbit, which occurs one-half period of the transfer orbit after the first maneuver.

transferOrbitDuration = ...
    pi*sqrt((semiMajorAxisTransfer^3)/standardGravitationalParameter);
circularizationBurnTime = injectionBurnTime + seconds(transferOrbitDuration);

Specify a name for the maneuver.

circularizationBurnName = "Circularization burn";

Calculate Propellant Consumption

Given the initial mass $m_0$, specific impulse $I_{sp}$ of the propulsion system, and the total delta-V $\Delta v$, compute the final mass $m_f$ after applying the delta-V using the Tsiolkovsky rocket equation.

$$m_f=\frac{m_0}{e^{\frac{\Delta v}{g_0{I}_{sp}}}}$$

where $g_0$ is the standard acceleration of gravity. The difference between the two masses is the total mass of propellant consumed ($m_p$).

$$m_p=m_0-m_f$$

Assume the mass of the spacecraft prior to the injection burn to be 6,000 kg.

initialMass = 6000; % kg

Assume specific impulse of the propulsion system to be 310 s.

specificImpulse = 310; % s

Calculate the mass of the spacecraft after the execution of both maneuvers.

g0 = 9.80665; % m/s^2
deltaVMagnitudeTotal = injectionBurnDeltaVMagnitude + circularizationBurnDeltaVMagnitude;
finalMass = initialMass/exp(deltaVMagnitudeTotal/(g0*specificImpulse))

finalMass = 
4.0129e+03

Calculate the mass of the propellant consumed.

propellantConsumed = initialMass - finalMass

propellantConsumed = 
1.9871e+03

Simulate and Visualize Hohmann Transfer in Satellite Scenario

Create a satellite scenario with a start time matching the mission start, a stop time when the spacecraft completes one revolution on the final orbit after the circularization burn, and a sample time of 60 seconds.

finalOrbitPeriod = ...
    2*pi*sqrt((semiMajorAxisFinal^3)/standardGravitationalParameter);
stopTime = circularizationBurnTime + seconds(finalOrbitPeriod);
sampleTime = 60; % s
sc = satelliteScenario(startTime,stopTime,sampleTime);

Add the spacecraft to the scenario based on the orbital elements corresponding to the initial orbit. Set the orbit propagator to "two-body-keplerian".

sat = satellite(sc, ...
    semiMajorAxisInitial, ...
    eccentricity, ...
    inclination, ...
    rightAscensionOfAscendingNode, ...
    argumentofPeriapsis, ...
    trueAnomaly, ...
    OrbitPropagator="two-body-keplerian", ...
    Name="Spacecraft");

Define a timetable out of the calculated maneuvers.

Time = [injectionBurnTime;circularizationBurnTime];
DeltaV = [injectionBurnDeltaV';circularizationBurnDeltaV'];
CoordinateFrame = [injectionBurnReferenceFrame;circularizationBurnReferenceFrame];
Name = [injectionBurnName;circularizationBurnName];
maneuvers = timetable(Time,DeltaV,CoordinateFrame,Name);

Schedule the maneuvers for the spacecraft.

impulsiveManeuver(sat,Insert=maneuvers)

ans=2×3 timetable
    14-Dec-2022 02:41:08          "Injection burn"    638.7907    0    0    "vnb"
    14-Dec-2022 03:46:08    "Circularization burn"    584.0904    0    0    "vnb"

Visualize the scenario using satelliteScenarioViewer.

v = satelliteScenarioViewer(sc, ...
    CameraReferenceFrame="Inertial", ...
    PlaybackSpeedMultiplier=2000);
campos(v, ...
    43, ...    % Latitude, deg
    -97.7, ... % Longitude, deg
    22409300); % Altitude, m

Set the orbit lead and trail times to span the entire scenario duration to visualize the whole trajectory.

sat.Orbit.LeadTime = seconds(stopTime - startTime);
sat.Orbit.TrailTime = seconds(stopTime - startTime);
sat.Orbit

ans = 
  Orbit with properties:

                          LeadTime: 1.9680e+04
                         TrailTime: 1.9680e+04
                         LineWidth: 1
                         LineColor: [0.0590 1 1]
    ShowImpulsiveManeuverDirection: 1
        ShowImpulsiveManeuverLabel: 1
                    VisibilityMode: 'inherit'

Play the scenario.

play(sc);