Thrust and the rocket equation - Rocket Propulsion notes

# this line makes figures interactive in Jupyter notebooks
%matplotlib inline
from matplotlib import pyplot as plt

Thrust¶

In the last module, the expression for thrust from a rocket (or any propulsion system relying on an exhaust jet) was given. Let’s now derive that, by considering a rocket undergoing a static test fire (where it is held in place by a test stand that counteracts the thrust force).

Let’s draw a control surface just around the rocket, with the exhaust jet crossing the boundary at the exit of hte nozzle, and apply conservation of momentum in the horizontal direction:

\sum F_x = \dot{m} V_{\text{out}} - \cancel{\dot{m} V_{\text{in}}} \\ R_x + \cancel{D} + A_e (p_a - p_e) = \dot{m} V_e \;,

(1)

where the reaction force $R_x$ imparted by the test stand (to hold the rocket in place) is the same as the thrust force: $| R_x | = T$ . Then, for thrust, we have:

\begin{align} T &= \dot{m} V_e + A_e (p_e - p_a) \\ &= T_m + T_p \;, \end{align}

(2)

where $T_m$ refers to the momentum flux thrust and $T_p$ refers to the pressure thrust terms. If $p_e < p_a$ , or the nozzle is overexpanded, then the pressure thrust is negative and reduces the overall thrust. Typically, though, nozzles are designed such that $p_e \geq p_a$ .

In the vacuum of space (or near-vacuum),

T = \dot{m} V_e + p_e A_e \;.

(3)

Effective exhaust velocity¶

The effective exhaust velocity, $c$ , is defined as the actual thrust divided by the mass flow rate:

\begin{align} c &= \frac{T}{\dot{m}} = \frac{\dot{m} V_e + A_e (p_e - p_a)}{\dot{m}} \\ c &= V_e + \frac{A_e (p_e - p_a)}{\dot{m}} \;. \end{align}

(4)

This allows us to express thrust as

T = \dot{m} c

(5)

and we can relate effective exhaust velocity to specific impulse:

c = I_{\text{sp}} g_0 \;.

(6)

The rocket equation¶

We now have an expression for thrust of a rocket, considering static firing on a test stand. What about a rocket in motion, which is continuously exhausting propellant mass to generate thrust? If we apply Netwon’s second law, we can obtain the Tsiolkovsky rocket equation.

Consider a rocket in motion in some arbitrary flight direction, where we assume

thrust in the direction of flight (no thrust vectoring)
steady thrust, and
no lift generated by the rocket body.

Applying Newton’s second law:

m \frac{dV}{dt} = \sum F_F = T - D - m g \sin \theta \,;

(7)

where the $_F$ subscript indicates in the direction of flight, and

m V \frac{d\theta}{dt} = \sum F_{\perp} = -mg \cos \theta \;,

(8)

where the $_{\perp}$ subscript indicates in the direction normal to flight. Now, let’s insert our expression for thrust $T = \dot{m}c$ into Equation (7), and recall that the rate of change of the rocket’s mass is related to the mass flow rate of exhausting propellant: $-\dot{m} = \frac{dm}{dt}$ :

\begin{align} m \frac{dV}{dt} &= \dot{m} c - D - m g \sin \theta \\ \frac{dV}{dt} &= \frac{\dot{m}}{m} c - \frac{D}{m} - g \sin \theta \end{align}

(9)

\rightarrow dV = -\frac{dm}{m} c - \frac{D}{m} dt - (g \sin \theta) dt

(10)

We can integrate this from time $t = 0$ to any time $t_b$ , or the burnout time when the rocket stops thrusting and no longer exhausts propellant:

V_f - V_0 = \Delta V = -c \int_{m_0}^{m_f} \frac{dm}{m} - \int_0^{t_b} \frac{D}{m} dt - \int_0^{t_b} (g \sin \theta) dt

(11)

or, to obtain Tsiolkovsky’s rocket equation:

\boxed{\Delta V = c \ln \frac{m_0}{m_f} - \int_0^{t_b} \frac{D}{m(t)} dt - \int_0^{t_b} (g(t) \sin \theta) dt} \;,

(12)

where $m_0$ refers to the initial vehicle mass, $m_f$ refers to the final vehicle mass at time $t_b$ , and in general drag force, gravity, and flight direction may vary with time: $D = D(t)$ , $g = g(t)$ , and $\theta = \theta(t)$ .

If we neglect the gravity and drag terms, we can simplify:

\Delta V = c \ln \frac{m_0}{m_f} \;.

(13)

and, we can define the rocket mass ratio ( $R$ or MR):

R = \frac{m_f}{m_0} \;.

(14)

In addition, recalling the definition of effective exhaust velocity $c = I_{\text{sp}} g_0$ :

\Delta V = I_{\text{sp}} g_0 \ln \frac{1}{R} \;,

(15)

which is the famous form of the rocket equation. This relates specific impulse of the propulsion system (which comes from propellant choice), mass ratio of the rocket, and the resulting total change in velocity, which is driven by mission requirements.

Example: Numerical integration of vertical launch¶

Determine the burnout velocity and maximum height achieved by launching the German V-2 rocket, assuming a vertical launch. Take into account drag, the variation of atmospheric density with altitude, and the variation of gravity with altitude. Use a forward difference numerical scheme (i.e., forward Euler) to solve the system of equations.

Given data:

specific impulse: $I_{\text{sp}}$ = 250 sec
initial mass: $m_0$ = 12700 kg
propellant mass: $m_p$ = 8610 kg
burn time: $t_b$ = 60 sec
maximum body diameter: $D$ = 5 feet 4 inches (5.333 ft)
height: $H$ = 46 feet

Figure 1 shows how the coefficients of lift and drag vary with Mach number and angle of attack.

V2 rocket coefficients of lift and drag — Figure 1:Variation of coefficients of drag and lift with Mach number and angle of attack, for the German V-2 rocket. Source: Sutton & Biblarz (2016).

Let’s set up the system of ordinary differential equations for velocity ( $m$ ), altitude ( $h$ ), and mass ( $m$ ) that describes this problem:

\begin{align*} \frac{dv}{dt} &= -\frac{I_{\text{sp}} g_0}{m} \frac{dm}{dt} - \frac{1}{2} \rho (h) v^2 A C_D (h) - g(h) \\ \frac{dh}{dt} &= v \\ \frac{dm}{dt} &= -\dot{m} = -\frac{m_p}{t_b} \;, \end{align*}

(16)

where we obtained the last equation by recognizing that

\begin{align*} m_f &= m_0 - \dot{m} t_b \\ m_p &= m_0 - m_f = \dot{m} t_b \\ \dot{m} &= \frac{m_p}{t_b} \end{align*}

(17)

for steady burning. In the above, $g_0$ is the reference acceleration due to gravity (9.80065 m/s²), $\rho(h)$ is the atmospheric density, $A$ is the cross-sectional area of the body, $C_D$ is the drag coefficient, and $g(h)$ is the (varying) acceleration due to gravity.

Density and gravity: To move forward, we need to describe how to evaluate density, gravitational acceleration, and the coefficient of drag at any instant in time. The first two are known functions of altitude:

\begin{align*} \rho (h) &= \rho_0 e^{-h / H_n} \\ g(h) &= g_0 \left( \frac{R_0}{R_0 + h} \right)^2 \;, \end{align*}

(18)

where $\rho_0$ = 1.225 kg/m³ is the density at sea level, $H_n$ = 10.4 km is the height scale of the exponential fall, $R_0$ = 6378.388 km is the mean radius of Earth, and $g_0$ = 9.80665 m/s² is the sea-level acceleration due to gravity.

Coefficient of drag¶

Coefficient of drag is more difficult to handle, since it depends on angle of attack and Mach number: $C_D = C_D(M, \alpha)$ . For vertical flight, the angle of attack is zero ( $\alpha = 0$ ), but the Mach number will be changing both as the vehicle speed changes and the atmospheric conditions change with altitude:

M = \frac{v}{\sqrt{ \gamma p(h) / \rho(h) }}

(19)

Option 1: As a first approximation, we could simply assume a constant value for coefficient of drag. For example, before and after the transonic regime (where Mach number approaches 1), $C_D \approx 0.15$ .

Option 2: If we want to more-accurately capture the increase in coefficient of drag as the vehicle approaches a sonic velocity, and then the relaxation when supersonic, we need to capture the dependence of $C_D$ on Mach number.

First, we already have an expression for $\rho(h)$ , but we also need an expression for pressure. We can use a similar exponential approximation:

p(h) = p_0 e^{-h / H_p} \;,

(20)

where $p_0$ = 101325 Pa is the sea-level standard atmospheric pressure and $H_p$ = 8.4 km. With $\gamma = 1.4$ being a reasonable constant, we can now calculate Mach number based on velocity and altitude.

Now, we need to obtain a functional dependence of $C_D (M)$ ; based on Figure 1, we see that the relationship is fairly complex. The “easiest” way to create a function for $C_D (M)$ is to:

Extract data points from the plot manually, using a tool such as WebPlotDigitizer
Perform a spline interpolation on the data points, using a function such as spline in Matlab or scipy.interpolate.UnivariateSpline in Python.

Let’s see this in action, using data extracted manually from the source figure:

import numpy as np
from scipy.interpolate import UnivariateSpline

data = np.genfromtxt('coefficient-drag.csv', delimiter=',')

coefficient_drag = UnivariateSpline(data[:,0], data[:,1], s=0, ext='const')

mach = np.linspace(0.2, 6, num=200, endpoint=True)
plt.plot(data[:,0], data[:,1], 'o', mach, coefficient_drag(mach), '-')
plt.xlabel('Mach number')
plt.ylabel('Coefficient of drag')
plt.legend(['Data', 'Cubic spline'])
plt.show()

Loading...

In Matlab, the equivalent fit can be done with:

data = readmatrix('coefficient-drag.csv');
plot(data(:,1), data(:,2), 'o'); hold on

coef_drag = spline(data(:,1), data(:,2));
xx = linspace(0.2, 5, 500);
plot(xx, ppval(coef_drag,xx), '-');

Now that we have a fully defined system of equations, we need a numerical method to integrate.

We can first use a forward difference scheme, also known as the forward Euler method:

\mathbf{z}_{i+1} = \mathbf{z}_i + \Delta t \, \mathbf{f} \left( t_i, \mathbf{z}_i \right) \;,

(21)

where $\mathbf{z}_{i}$ is the vector of state variables at time $t_i$ , $\Delta t$ is the time step size, and $\mathbf{f}$ represents the vector of time derivatives.

To implement, we need to create a function to evaluate the time derivatives, and then specify initial conditions for velocity, altitude, and mass: $v(t = 0) = 0$ , $h(t = 0) = 0$ , and $m(t = 0) = m_0$ .

from scipy.integrate import solve_ivp

def vertical_launch(t, y, spec_impulse, mass_prop, time_burn, diameter, coef_drag):
    '''Evaluates system of time derivatives for velocity, altitude, and mass.
    '''
    radius_earth = 6378.388 * 1e3
    gravity_ref = 9.80665
    density_ref = 1.225
    pressure_ref = 101325
    gamma = 1.4
    height_den = 10400
    height_pres = 8400
    
    v = y[0]
    h = y[1]
    m = y[2]
    
    gravity = gravity_ref * (radius_earth / (radius_earth + h))**2
    density = density_ref * np.exp(-h / height_den)
    pressure = pressure_ref * np.exp(-h / height_pres)
    
    mach = v / np.sqrt(gamma * pressure / density)
    area = np.pi * diameter**2 / 4
    drag = 0.5 * density * v**2 * area * coef_drag(mach)
    
    dmdt = -mass_prop / time_burn
    dvdt = (-spec_impulse * gravity_ref / m) * dmdt - drag / m - gravity
    dhdt = v
    
    return [dvdt, dhdt, dmdt]

# given constants
spec_impulse = 250
mass_initial = 12700
mass_propellant = 8610
time_burn = 60
diameter = 1.626

sol = solve_ivp(
    vertical_launch, [0, time_burn], [0, 0, mass_initial], method='RK45', 
    args=(spec_impulse, mass_propellant, time_burn, diameter, coefficient_drag),
    dense_output=True
    )

time = np.linspace(0, time_burn, 100)
z = sol.sol(time)

fig, axes = plt.subplots(2, 1)

axes[0].plot(time, z[0, :])
axes[0].set_ylabel('Velocity (m/s)')
axes[0].grid(True)

axes[1].plot(time, z[1, :] / 1000)
axes[1].set_ylabel('Altitude (km)')
axes[1].set_xlabel('Time (s)')
axes[1].grid(True)

plt.show()

print(f'Burnout velocity: {z[0,-1]: 5.2f} m/s')
print(f'Burnout altitude: {z[1,-1]/1000: 5.2f} km')

Loading...

Burnout velocity:  1956.11 m/s
Burnout altitude:  44.54 km

In Matlab, the above could be done with:

clear; clc;

g0 = 9.80665;
Re = 6378.388 * 1e3;

Isp = 250;
m0 = 12700;
mp = 8610;
tb = 60;
D = 1.626;

data = readmatrix('coefficient-drag.csv');
plot(data(:,1), data(:,2), 'o'); hold on

coef_drag = spline(data(:,1), data(:,2));
xx = linspace(0.2, 5, 500);
plot(xx, ppval(coef_drag,xx), '-');

f = @(t,z) rocket(t, z, Isp, mp, tb, D, coef_drag);

[T, Z] = ode45(f, [0 tb], [0 0 m0]);

function dzdt = rocket(t, z, Isp, mp, tb, D, coef_drag)
g0 = 9.80665; Re = 6378.388 * 1e3;
mdot = mp / tb;

v = z(1);
h = z(2);
m = z(3);

g = g0 * (Re / (Re + h))^2;

rho0 = 1.225;
p0 = 101325;
Hn = 10400;
Hp = 8400;
gamma = 1.4;
rho = rho0 * exp(-h / Hn);
pres = p0 * exp(-h / Hp);

area = pi * D^2 / 4;
M = v / sqrt(gamma * pres / rho);
C_D = ppval(coef_drag, M);
drag = 0.5 * rho * v^2 * area * C_D;

dzdt = zeros(3,1);
dzdt(3) = -mdot;
dzdt(1) = (-Isp*g0/m)*dzdt(3) - drag/m - g;
dzdt(2) = v;

end