Rocket

A rocket launch simulation using a gravity turn maneuver using differential equations as described by this article in Mintoc (see below for a further description of the differential equations). The Casadi solution is based on the code given in this link, updated by Mirko Hahn and in the github repository gturn.

After setting up the environment (the program was tested in Python 3.10),

python -m venv ./venv
pip install -r requirements.txt
source ./venv/bin/activate
(or in windows)
./venv/scripts/activate

run the programs as follows:

Calculation of the thrust control:

python rocket_casadi_solution.py mintoc_20T_1.cfg

This program will create an excel trhust control file as defined in the config file. Then run the rocket launch program using this thrust control solution

python rocket_launch.py mintoc_20T_1.cfg

Where the configuration file mintoc_20T_1.cfg is

# parameters must be given in this order
rocket dry mass (kg)           m_dry  : 2_000
rocket fuel mass (kg)          m_fuel : 20_000
motor impulse (s) zero alt     Isp0   : 450
motor impulse (s) vacuum       Isp1   : 450
maximum thrust (N)             Fmax   : 600_000
rocket reference area (m^2)    A      : 3.5
initial velocity (m / s)       vel    : 1e-3
initial flight angle (radians) beta   : 0.000274
initial altitude (m)           h      : 0
gravity at zero altitude       g0     : 9.81
radius at altitude zero        r0     : 6_000_000
drag coefficient (_)           cd     : 0.75
scale height (m)               H      : 8_500
density at zero altitude       rho    : 1.2230948554874
number of shooting intervals   N      : 300
altitude objective (m)         h_obj  : 200_000
velocity objective (m/s)       v_obj  : 8_000
beta objective (degrees)       q_obj  : 90
control file (.xlsx)                  : mintoc_gravity_turn_20T_1.xlsx
time interval (s)  1.1268             : 1
status_update_step                    : 120
duration (s)                          : 20_000
speed min_max (m/s)                   : 0, 10_000
flight angle min max (degrees)        : 0, 110
altitude min max (m)                  : 0, 2_000_000
h_range min max (degrees)             : 0, 1440
acceleration min max (m*s-2)          : -100, 170
rocket sprite file                    : rocket_sprite2.png

Once the program starts it gives main parameters in the console and a display of graphs and trajectory.

Gravity turn

The gravity turn or zero lift turn is a common maneuver used to launch spacecraft into orbit from bodies that have non-negligible atmospheres. The goal of the maneuver is to minimize atmospheric drag by always orienting the vehicle along the velocity vector. In this maneuver, the vehicle's pitch is determined solely by the change of the velocity vector through gravitational acceleration and thrust. The goal is to find a launch configuration and thrust control strategy that achieves a specific orbit with minimal fuel consumption.

Physical description and model derivation

For the purposes of this model, we start with the following ODE system proposed by Culler et. al. in [Culler1957]:

where v is the speed of the vehicle, g is the gravitational acceleration at the vehicle's current altitude, F is the accelerating force and β is the angle between the vertical and the vehicle's velocity vector. In the original version of the model, the authors neglect aspects of the problem:

• Variation of g over altitude
• Decrease of vehicle mass due to fuel consumption
• Curvature of the surface
• Atmospheric drag

Changes in gravitational acceleration

To account for changes in g, we make the following substitution:

where g₀ is the gravitational acceleration at altitude zero and r₀ is the distance of altitude zero from the center of the reference body.

Decrease in vehicle mass

To account for changes in vehicle mass, we consider m a differential state with the following derivative:

where I_sp is the specific impulse of the vehicle's engine. Specific impulse is a measure of engine efficiency. For rocket engines, it directly correlates with the engine's exhaust velocity and may vary with atmospheric pressure, velocity, engine temperature and combustion dynamics. For the purposes of this model, we will assume it to be constant.

The vehicle's fuel reserve is modelled by two parameters: m₀ denotes the launch mass (with fuel) and m₁ denotes the dry mass (without fuel).

Curvature of the reference body's surface

To accomodate the reference body's curvature, we introduce an additional differential state θ which represents the change in the vehicle's polar angle with respect to the launch site. The derivative is given by

Note that the vertical changes as the vehicle moves around the reference body meaning that the derivative of β must be changed as well:

.

Atmospheric drag

To model atmospheric drag, we assume that the vehicles drag coefficient c_d is constant. The drag force is given by

where ρ is the density of the atmosphere and A is the vehicle's reference area. We assume that atmospheric density decays exponentially with altitude:

where ρ₀ is the atmospheric density at altitude zero and H is the scale height of the atmosphere. The [drag force] is introduced into the acceleration term:

.

Note that if the vehicle is axially symmetric and oriented in such a way that its symmetry axis is parallel to the velocity vector, it does not experience any lift forces. This model is simplified. It does not account for changes in temperature and atmospheric composition with altitude. Also, c_d varies with fluid viscosity and vehicle velocity. Specifically, drastic changes in c_d occur as the vehicle breaks the sound barrier. This is not accounted for in this model.