Example: JP4/RFNA Rocket Motor¶

A liquid propellant rocket motor using JP-4/Red Fuming Nitric Acid (RFNA) propellant combination must produce 19.6 kN of thrust at sea level. The characteristics of the combustion products are \(T_c\) = 3000 K, \(\gamma\) = 1.2 (constant), \(W\) = 27 kg/kmol (constant), and \(p_c\) = 19.34 atm. Determine:

the throat and exit diameters, in cm, for matched nozzle exit flow at sea level assuming a nozzle efficiency of \(\eta_n\) = 95%
the characteristic velocity, \(c^*\), in m/s
the propellant mass flow rate, in kg/s
the specific impulse of the engine, in s
the thrust, in kN, at an altitude of 11.5 km where the pressure is \(p_0\) = 0.2 atm
the thrust, in kN, and exit diameter, in cm, assuming matched operation at this altitude and the same throat diameter

Solution¶

from pint import UnitRegistry
import numpy as np

units = UnitRegistry()

F_sl = 19.6 * units.kN
T_c = 3000 * units.K
gamma = 1.2 * units.dimensionless
W = 27 * units.kg / units.kmol
R = units.molar_gas_constant / W
c_p = gamma * R / (gamma - 1)
p_c = 19.34 * units.atmosphere
eta_n = 0.95 * units.dimensionless
p_sl = 1 * units.atmosphere
p_0 = 0.2 * units.atmosphere

Engineering Model¶

The nozzle is not isentropic, but is characterized by the nozzle efficiency
The specific heat ratio and molecular weight are constant
The discharge coefficient and \(\lambda\) are equal to one
The flow is adiabatic through the nozzle, even though it is not isentropic

Procedure¶

In this problem, we are initially given \(p_0\) and told the operation of the nozzle is matched. We are also given \(p_c\), so we can use the stagnation-to-static pressure ratio for isentropic flow to determine the ideal exit velocity from Eq. (46). Using the nozzle efficiency from Eq. (44), we can correct \(V_{e,i}\) to find the actual nozzle exit velocity.

Since we are given the thrust and initially the pressure is matched, we can use the actual nozzle exit velocity to calculate the mass flow rate from Eq. (36). We can also use the actual nozzle exit velocity to find the actual exit temperature from the energy equation. This determines the actual exit density via the ideal gas law. We can then find the exit area from the definition of the mass flow rate.

Then, assuming the discharge coefficient is unity, the mass flow rate coefficient, \(c_m\), is equal to the inverse of the characteristic velocity \(c^*\). The characteristic velocity is given in terms of the combustion chamber speed of sound at the specific heat ratio by Eq. (41). The mass flow rate coefficient can be used to find the throat area.

Next, the specific impulse can be determined by dividing the actual exit velocity by the gravitational constant.

At the altitude of 11.5 km, we need to correct the thrust by including the pressure thrust in Eq. (36).

Finally, we need to calculate the thrust and exit diameter at the altitude of 11.5 km assuming matched operation. This can be done by the thrust coefficient, \(c_F\), assuming that \(\lambda = c_d = 1\). Since the throat area and chamber conditions are the same, the mass flow rate is the same. This can be used to find the actual exit velocity. From there, the same process as before is used to find the exit diameter.

Work¶

pe_pc = (p_sl / p_c).to("dimensionless")
V_ei = np.sqrt(2 * gamma / (gamma - 1) * R * T_c * (1 - pe_pc**((gamma - 1) / gamma))).to("m/s")

Ve_sl = V_ei * np.sqrt(eta_n)
mdot = (F_sl / Ve_sl).to("kg/s")

Te_sl = T_c - Ve_sl**2 / (2 * c_p)
rhoe_sl = (p_sl / (R * Te_sl)).to("kg/m**3")
Ae_sl = (mdot / (rhoe_sl * Ve_sl)).to("cm**2")
de_sl = np.sqrt(4 * Ae_sl / units.pi).to("cm")

Gamma = np.sqrt(gamma) * (2 / (gamma + 1))**((gamma + 1) / (2*(gamma - 1)))
a_c = np.sqrt(gamma * R * T_c)
c_star = (a_c / (Gamma * np.sqrt(gamma))).to("m/s")
c_m = 1 / c_star
A_th = (mdot / (p_c * c_m)).to("cm**2")
d_th = np.sqrt(4 * A_th / units.pi).to("cm")

I_sp = (Ve_sl / units.gravity).to("s")

F = (mdot * Ve_sl + Ae_sl * (p_sl - p_0)).to("kN")

c_V = np.sqrt(eta_n)
c_F = c_V * Gamma * np.sqrt(2 * gamma / (gamma - 1) * (1 - (p_0/p_c)**((gamma - 1) / gamma)))
F_alt = (c_F * ((19.34 * 101325) * units.Pa) * A_th).to("kN")
Ve_alt = (F_alt / mdot).to("m/s")
Te_alt = T_c - Ve_alt**2 / (2 * c_p)
rhoe_alt = (p_0 / (R * Te_alt)).to("kg/m**3")
Ae_alt = (mdot / (rhoe_alt * Ve_alt)).to("cm**2")
de_alt = np.sqrt(4 * Ae_alt / units.pi).to("cm")

Results¶

Property	Value	Property	Value	Property	Value
\(V_{e,i}\)	2078.33 m/s	\(d_e\)	18.69 cm	\(F\) at 11.5 km, matched	22.93 kN
\(V_e\)	2025.71 m/s	\(c^*\)	1482.06 m/s	\(V_e\) at 11.5 km, matched	2369.78 m/s
\(\dot{m}\)	9.68 kg/s	\(A_{\text{th}}\)	73.18 cm²	\(T_e\) at 11.5 km, matched	1480.27 K
\(T_e\)	1889.54 K	\(d_{\text{th}}\)	9.65 cm	\(\rho_e\) at 11.5 km, matched	0.04 kg/m³
\(\rho_e\)	0.17 kg/m³	\(I_{\text{sp}}\)	206.56 s	\(A_e\) at 11.5 km, matched	918.41 cm²
\(A_e\)	274.29 cm²	\(F\) at 11.5 km	21.82 kN	\(d_e\) at 11.5 km, matched	34.20 cm

Notes for Propulsion