Contents

Example: LR91 Rocket Motor

The LR91 rocket motor uses the liquid bipropellants nitrogen tetroxide (N2O4) and Aerozine 50 (a 50/50 mixture by weight of hydrazine and unsymmetrical dimethyl hydrazine). Assume that the LR91 is an ideal rocket motor with nozzle expansion ratio of \(\varepsilon\) = 49.2, chamber pressure of \(p_c\) = 5.6 MPa, chamber temperature of \(T_c\) = 3,400 K. The exhaust gas products of combustion have a ratio of specific heats \(\gamma\) = 1.3 and a gas constant of \(R\) = 390.4 J/(kg K). The nozzle exit diameter is \(d_e\) = 1.6 m. Determine:

a. the altitude at which the nozzle is designed for optimal expansion (\(p_e = p_0\)), in km and ft b. the nozzle throat diameter (\(d_{\text{th}}\)), in cm and in. c. the mass flow rate (\(\dot{m}\)), in kg/s and lb/s d. the ideal thrust at the optimal altitude, in kN and lbf

Solution

from pint import UnitRegistry
import numpy as np
from scipy import optimize
from ambiance import Atmosphere

units = UnitRegistry()

T_c = 3400 * units.K
p_c = 5.6 * units.MPa
d_e = 1.6 * units.m
epsilon = 49.2 * units.dimensionless
gamma = 1.3 * units.dimensionless
W = 21.30 * units.kg / units.kmol
R = units.molar_gas_constant / W

Engineering Model

  1. Ideal rocket motor, all processes are isentropic

  2. Expansion in the nozzle is optimized, so \(p_e = p_0\)

  3. The specific heat ratio and molecular weight are constant

Procedure

To determine the altitude, we need to determine the ambient pressure around the rocket motor. With the ambient pressure, we can look up or find the altitude in standard atmospheric tables.

Since we are looking for the optimal expansion, we can find the ambient pressure by finding the exit pressure of the nozzle, since they are equal. The exit pressure of the nozzle is going to be found from the isentropic equations for compressible flow, since we are assuming this is an ideal nozzle.

The isentropic equations for compressible flow require the exit Mach number, which we will find from the expansion ratio. Likewise, the nozzle throat diameter is determined by the expansion ratio and the exit diameter.

Finally, the mass flow rate can be found from the maximum flow rate equations, the exit velocity from the exit Mach number, and the ideal thrust by multiplying these two quantities.

Work

We have two ways of solving for the exit Mach number from the expansion ratio. One is to use the area-Mach-number relation, Eq. (94), which we have seen before. The other is to introduce a quantity called the X-function, defined as:

(55)\[X = M \sqrt{\gamma} \left[1 + \frac{\gamma - 1}{2} M^2\right]^{\frac{-\left(\gamma + 1\right)}{2\left(\gamma - 1\right)}}\]

Notice that the \(X\)-function is equivalent to the \(\Gamma\) (Gamma) function we defined in Eq. (40) when \(M = 1\). Whether you use the area-Mach-number relation, or the \(X\)-function, you’ll get the same answer for \(M\) at the end.

In terms of the \(X\)-function, the ideal maximum mass flow rate can be written as:

\[\dot{m} = \frac{p_t A X}{\sqrt{R T_t}}\]

Since the mass flow rate in the nozzle must be constant, we can plug in the throat quantities and the exit quantities and equate them. The total pressure and temperature are constant, since the nozzle process is isentropic. In addition, at the throat, the Mach number is one, so \(X_{\text{th}} = \Gamma\).

\[\frac{p_{t} A_{\text{th}} \Gamma}{\sqrt{R T_t}} = \frac{p_t A_e X_e}{\sqrt{R T_t}}\]

From this equation, the gas constant, total pressure, and total temperature will cancel, leaving:

\[X_e = \frac{1}{\varepsilon}\Gamma = \frac{1}{\varepsilon}\left[\sqrt{\gamma}\left(\frac{2}{\gamma + 1}\right)^{\frac{\gamma + 1}{2\left(\gamma - 1\right)}}\right]\]
X_e = 1 / epsilon * (np.sqrt(gamma) * (2 / (gamma + 1))**((gamma + 1)/(2*(gamma - 1))))

This gives \(X_e\) = 0.0136. Then, the exit Mach number is found from Eq. (55). This function is transcendental in \(M_e\), so we use the scipy.optimize.brentq() root finding method.

def X_function(M, X, gamma):
    power = -(gamma + 1) / (2 * (gamma - 1))
    return M * np.sqrt(gamma) * (1 + (gamma - 1) / 2 * M**2)**(power) - X
M_e = optimize.brentq(
    X_function, 
    1.1,
    10.0,
    args=(X_e.magnitude, gamma.magnitude)
) * units.dimensionless

This gives an exit Mach number of \(M_e\) = 5.07. Next, we can determine the ratio of the stagnation to the static pressure at the exit from Eq. (77):

\[\frac{p_t}{p_e} = \frac{p_c}{p_0} = \left(1 + \frac{\gamma - 1}{2} M_e^2\right)^{\frac{\gamma}{\gamma - 1}}\]
pc_p0 = (1 + (gamma - 1) / 2 * M_e**2)**(gamma / (gamma - 1))
p_0 = (p_c / pc_p0).to("Pa")

This gives a pressure ratio of 938.4786 and an ambient pressure of 5967.10 Pa. Then, we can use ambiance to determine the altitude:

atmos = Atmosphere.from_pressure(p_0.magnitude)
h = atmos.h[0] * units.m

Finally, the altitude is \(h\) = 19.51 km = 64021 ft.

The throat diameter is found from the expansion ratio and the known exit diameter:

\[\frac{A_e}{A_{\text{th}}} = \frac{d_e^2}{d_{\text{th}}^2} \Rightarrow d_{\text{th}} = d_e\sqrt{\frac{1}{\varepsilon}}\]
d_th = d_e * np.sqrt(1 / epsilon)

The diameter is \(d_{\text{th}}\) = 22.8 cm = 9.0 in. Then, the mass flow rate can be determined from the exit conditions:

A_e = np.pi * d_e**2 / 4
mdot = p_c * A_e * X_e / np.sqrt(R * T_c)

The mass flow rate is \(\dot{m}\) = 132.6 kg/s = 292.2 lb/s. Finally, we can find the exit velocity in a few ways. Let’s compute the speed of sound from the static temperature at the exit, and then use the Mach number to find the exit velocity.

Tt_Te = 1 + (gamma - 1) / 2 * M_e**2
T_e = T_c/Tt_Te
a_e = np.sqrt(gamma * R * T_e)
V_e = M_e * a_e

The exit static temperature is \(T_e\) = 700.7 K, the speed of sound is \(a_e\) = 596.30 m/s, and the exit velocity is \(V_e\) = 3021.89 m/s. Finally, the thrust is computed from Eq. (36), using the fact that the exit pressure is matched to the ambient.

F = mdot * V_e

The thrust is \(F\) = 400.55 kN = 90048.15 lbf. Alternately, we can calculate the thrust using the ideal thrust coefficient, given by Eq. (48).

Gamma = np.sqrt(gamma) * (2 / (gamma + 1))**((gamma + 1)/(2*(gamma - 1)))
p_e = p_0
pe_pc = (p_e / p_c).to("dimensionless")
c_fi = Gamma * np.sqrt(2 * gamma / (gamma - 1) * (1 - pe_pc**((gamma - 1) / gamma)))
A_th = np.pi * d_th**2 / 4
F_2 = A_th * p_c * c_fi

With this approach, the thrust is \(F\) = 400.55 kN, equivalent to the first approach. This provides a good check that our work is correct and consistent.

Design of Rocket Nozzles Example: JP4/RFNA Rocket Motor