Simulating a Rocket Nozzle in 150 Lines of Python

Introduction
Point a garden hose at your car, press your thumb over the end, and the water leaves faster. Narrow the opening, and the water exits faster. That's the intuition almost everyone carries about nozzles, and for an ordinary garden hose, it's a useful one.

A rocket nozzle works differently. It narrows to a waist in the middle, the throat, and then flares back open toward the exit. And it is precisely in that widening, flaring section that the exhaust does most of its accelerating, screaming past it faster and faster the more room it is given.

Today I want to explore this using Python.

Particles accelerating through a de Laval nozzle, dark and slow in the chamber, bright and fast at the supersonic exit. — Figure 1. Tracer particles moving through the solved velocity field. They crawl in the chamber, funnel through the throat, and accelerate as the nozzle widens. Colour is speed: dark is slow, bright is fast.

Trading Heat for Speed

The shape in Figure 1 is a de Laval nozzle: a converging section, a throat, then a diverging section. To understand why it is built this way, we need one idea and one equation.

The idea is a trade. Combustion fills the chamber with hot, high-pressure, nearly motionless gas. The nozzle's only job is to convert that stored thermal energy and pressure into directed kinetic energy.

For a steady, adiabatic nozzle with no shaft work, the energy equation is

c_p T_0 = c_p T + \frac{u^2}{2}

where:

$T_0$ is the stagnation temperature, the temperature the gas would have if brought to rest adiabatically;
$T$ is the local static temperature of the flowing gas;
$c_p$ is the specific heat at constant pressure;
$u$ is the local flow velocity.

The equation shows the energy trade directly: as

u

increases,

T

must decrease, assuming

T_0

remains constant.

For an ideal gas with constant specific heat,

h = c_p T

where:

$h$ is the local static enthalpy of the gas;
$c_p$ is the specific heat at constant pressure;
$T$ is the local static temperature.

For an ideal gas with constant

c_p

, the equation expresses the gas's enthalpy directly in terms of its temperature.

The consequence is simple: as the gas accelerates, its static temperature falls. Thermal energy is being converted into kinetic energy¹.

But that still does not explain the nozzle's geometry.

The following equation connects two things: how the nozzle's width is changing at a given point, and how the gas's speed is changing there. In other words, it answers the question "if I make the pipe a little wider or narrower here, does the gas speed up or slow down?"

\frac{dA}{A} = (M^2 - 1)\,\frac{du}{u}

where:

$\frac{dA}{A}$ is the fractional change in the nozzle's cross-sectional area $A$ , positive where the nozzle widens and negative where it narrows;
$\frac{du}{u}$ is the fractional change in the local flow speed $u$ ;
$M = u/a$ is the Mach number, the flow speed expressed as a multiple of $a$ , the local speed of sound ( $M < 1$ is subsonic, $M > 1$ is supersonic).

That little (M² − 1) term is where everything interesting hides. Here is how to read the equation without doing any algebra.

Suppose we want the gas to go faster, so its change in speed is positive. The equation then tells us what the area has to do to make that happen, and the answer flips completely depending on whether the gas is moving slower or faster than sound.

Take the slow case first. When the gas is subsonic (M < 1), the term (M² − 1) is negative, and working through the equation, that means the area has to get smaller for the gas to speed up. This is the behavior everyone already knows; put your thumb over the end of a garden hose and the water shoots out faster. Squeeze the flow and it accelerates. It is exactly why the nozzle narrows on the way in.

Now the fast case, which is where it gets strange. Once the gas is supersonic (M > 1), that same term flips to positive, and the whole rule inverts. To keep the gas accelerating, the area now has to get bigger. Widen the channel and supersonic gas speeds up; narrow it and it would slow back down. This is the exact opposite of the hose, and it is why a rocket nozzle flares open after its waist.

And right in between sits the turning point. At exactly M = 1, the coefficient (M² − 1) vanishes. For a smooth nozzle that passes continuously from subsonic to supersonic flow, the sonic point therefore occurs at the minimum-area section: the throat. In the idealized choked de Laval flow, the throat is the critical point where the flow reaches Mach 1.

Why does supersonic gas behave so strangely? Because mass must be conserved:

\dot{m} = \rho u A = \text{constant}

where

\rho

is density.

As the gas accelerates through the nozzle, its pressure and temperature fall. Its density falls too.

In supersonic flow, the density decrease associated with acceleration changes the continuity balance: the area must increase for the velocity to continue rising while mass flow remains constant.

The widening is not itself the source of the acceleration. The pressure gradient is. The geometry provides the area change required for the supersonic flow to keep accelerating while conserving mass and energy.

That distinction matters. It is tempting to think of the diverging wall as somehow pushing the gas faster. It is not. The nozzle shape constrains the flow, and the pressure field produced by that constraint accelerates the gas.

Mach number and velocity by position along nozzle. — Figure 2. Both Mach number and velocity keep rising after the throat, in the diverging section. The Mach curve crosses 1 exactly at the throat.

Choking

Once the throat reaches Mach 1, something worth pausing on happens: the flow chokes. The mass passing through the throat reaches a maximum fixed by the throat area and the chamber's stagnation conditions, and the choked throat becomes insensitive to further downstream pressure changes as far as the upstream mass flow is concerned.

The reason is about information. A pressure disturbance travels through a gas at the local speed of sound. At the sonic throat, the gas is moving at the local speed of sound, so an upstream-propagating acoustic disturbance has zero velocity relative to the throat. It cannot propagate upstream through the choked section.

Lower the pressure outside the nozzle further, and the downstream flow can respond, but the information cannot propagate back through the sonic throat and increase the upstream mass flux. The throat has reached the maximum mass flow allowed by the upstream state and its area.

So dropping the pressure outside the nozzle cannot coax more gas through a choked throat. The throat becomes a bottleneck, and the mass flow is fixed by the upstream stagnation conditions, the throat area, and the gas properties.

For an ideal gas, the choked mass flow can be written directly as

\dot{m} = A^* \frac{P_0}{\sqrt{T_0}} \sqrt{\frac{\gamma}{R}} \left( \frac{2}{\gamma+1} \right)^{\frac{\gamma+1}{2(\gamma-1)}}

where:

$\dot{m}$ is the mass flow rate through the choked nozzle;
$A^*$ is the throat area, the minimum cross-sectional area where the flow reaches $M = 1$ ;
$P_0$ is the chamber stagnation pressure;
$T_0$ is the chamber stagnation temperature;
$\gamma$ is the ratio of specific heats, $\gamma = c_p/c_v$ ;
$R$ is the specific gas constant.

Once the throat is choked,

\dot{m}

is determined by these quantities. Increasing

A^*

P_0

increases the mass flow rate, while increasing

T_0

decreases it.

The diverging section then takes that sonic flow and accelerates it to supersonic speed.

Solving the Flow, Simulating the Particles

I do not simulate particles that push on one another and bounce around pretending to be a fluid. That path is both more difficult to achieve and less accurate. Instead I solve the steady-state flow field first, under the assumptions below, and then release massless tracer particles to drift through it. The particles are a visualisation of the physics (and not a substitute).

Two assumptions make that field solvable.

The first is that the flow is quasi-one-dimensional. Velocity, pressure, temperature, and density are treated as uniform across a given slice of the nozzle, changing only along its length. A real nozzle has a faster core and slower gas near the walls; quasi-1D flattens each cross-section into a single set of numbers. It is the assumption that turns a three-dimensional flow into a curve I can solve one point at a time.

The second is that the flow is isentropic, which here means both adiabatic and reversible. Adiabatic means no heat passes through the nozzle walls. Reversible means there is no entropy generation from effects such as friction or shocks. Hold both, and the gas' entropy stays constant from chamber to exit. That is what makes the thermodynamic properties become clean functions of Mach number².

The field comes from inverting the area-Mach relation. Given the nozzle's area at each point, I need the Mach number there. The relation itself is

\frac{A}{A^\ast} = \frac{1}{M} \left[ \frac{2}{\gamma+1} \left( 1 + \frac{\gamma-1}{2}M^2 \right) \right]^{\frac{\gamma+1}{2(\gamma-1)}}

where:

$A$ is the local cross-sectional area of the nozzle;
$A^\ast$ is the throat area, the area at which $M = 1$ ;
$A/A^\ast$ is the local nozzle area measured relative to the throat area;
$M$ is the local Mach number, the ratio of flow speed to the local speed of sound;
$\gamma$ is the ratio of specific heats of the gas, $\gamma = c_p/c_v$ .

The equation relates the nozzle's local area to the Mach number of the flow. Given

M

and

\gamma

, it determines the required area ratio

A/A^\ast

. The throat is the reference point because

A = A^\ast

when

M = 1

For any area ratio greater than 1, the area-Mach relation has two valid solutions, one subsonic and one supersonic. The converging section takes the subsonic root, and the diverging section takes the supersonic root³.

It is a transcendental relation, so rather than hand-roll a fragile root-finder for every point, I sweep Mach number, tabulate the area ratio it produces, and invert by lookup, one table for each branch:

def area_mach(M):
    g = GAMMA
    return (1/M) * ((2/(g+1)) * (1 + 0.5*(g-1)*M**2))**((g+1)/(2*(g-1)))

# Build the inversion tables once: one subsonic branch, one supersonic.
M_sub, M_sup = np.linspace(1e-4, 1.0, 20000), np.linspace(1.0, 6.0, 20000)
AR_sub, AR_sup = area_mach(M_sub), area_mach(M_sup)

def mach_from_area(ar, supersonic):
    if supersonic:
        return np.interp(ar, AR_sup, M_sup)
    return np.interp(ar, AR_sub[::-1], M_sub[::-1])

From Mach number, the isentropic relations hand over temperature, the local speed of sound, and finally the velocity field u(x). The code is almost a literal translation:

def solve_field(T0):
    x = np.linspace(0, 1, 400)
    ar = area_ratio(x)

    M = np.where(
        x < X_THROAT,
        mach_from_area(ar, supersonic=False),
        mach_from_area(ar, supersonic=True)
    )

    T = T0 / (
        1 + 0.5*(GAMMA - 1)*M**2
    )

    a = np.sqrt(GAMMA * R_S * T)
    u = M * a

    return x, ar, M, T, u

The flow solver is essentially structured as a series of yields:

geometry → Mach number → temperature → speed of sound → velocity

This total is our flow solver and the particles can be added later.

Each particle carries a position and moves by the real local velocity, then recycles back to the chamber when it leaves:

def advect():
    px += np.interp(px, x, u) / u_max * STEP   # move by the real local velocity
    out = px > 1.0
    px[out] = 0.0                              # recycle a particle at the chamber
    pf[out] = rng.uniform(-1, 1, out.sum())

Everything else about the particles is presentation: their count, brightness, size, frame rate, and the way they are spread vertically.

The horizontal motion follows the computed velocity field after normalization for animation, while the vertical spread and other visual properties are illustrative.

Checking the Work

It is easy to write a solver that produces a plausible-looking picture. Although a resulting animation may look reasonable, we still want to test it for correctness:

First, the throat should be sonic:

print(
    f"throat Mach = "
    f"{M[np.argmin(ar)]:.4f}"
)

The numerical grid gives M = 1.003, which is close to the target of 1.0. The small offset comes from the finite sampling of the nozzle and the interpolation used for the area-Mach inversion.

Second, the exit velocity from the solver can be checked against the closed-form energy equation:

u_e = \sqrt{ 2c_pT_0 \left( 1 - \left( \frac{P_e}{P_0} \right)^{\frac{\gamma-1}{\gamma}} \right) }

where:

$u_e$ is the flow velocity at the nozzle exit;
$c_p$ is the specific heat at constant pressure;
$T_0$ is the chamber stagnation temperature;
$P_e$ is the static pressure at the nozzle exit;
$P_0$ is the chamber stagnation pressure;
$\gamma$ is the ratio of specific heats of the gas, $\gamma = c_p/c_v$ .

The equation gives the ideal exit velocity resulting from the pressure expansion through the nozzle. For fixed

c_p

T_0

, and

P_0

, a lower exit pressure

P_e

corresponds to greater expansion and therefore a higher exit velocity

u_e

The independently calculated exit velocity matches the velocity obtained from the Mach field to four significant figures.

Finally, conservation of mass gives us another useful field-level check:

\dot{m} = \rho u A = \text{constant}

where:

$\dot{m}$ is the mass flow rate through the nozzle;
$\rho$ is the local gas density;
$u$ is the local flow velocity;
$A$ is the local cross-sectional area of the nozzle.

The equation expresses conservation of mass: for steady flow, the same mass must pass through every cross-section of the nozzle per unit time. Therefore, if the area

A

changes, the velocity

u

and/or density

\rho

must change so that the mass flow rate

\dot{m}

remains constant.

The code evaluates

P = P0 * (T/T0)**(
    GAMMA/(GAMMA - 1)
)

rho = P / (R_S * T)

flux = rho * u * ar

and measures the spread of flux across the nozzle.

For this run, the mass-flux variation is only about 2.7 × 10⁻⁶ percent.

That is a useful consistency check. The area-Mach relation and the thermodynamic relations come from the same conservation laws, so this does not prove the model independently. What it does show is that the numerical inversion, interpolation, and property reconstruction remain internally consistent across the nozzle.

The animation can be wrong while still looking convincing. A conservation check gives us something much harder to fool.

The Combustion Knob

You will notice the nozzle burns no fuel. There is no flame in the model at all, and that is correct, because combustion never reaches into the nozzle as a process. It reaches in only through a set of conditions in the chamber: temperature, pressure, and the properties of the gas that combustion produced.

Change what you burn, and those chamber conditions change. In this simplified model, the main knob I expose is the chamber temperature T₀⁴.

So T₀ is the one dial that stands in for combustion. For fixed nozzle geometry and gas properties, increasing T₀ increases the exhaust velocity.

Figure 3. Exit velocity against chamber temperature. For a choked, quasi-one-dimensional, isentropic flow with fixed gas properties, the Mach number profile is determined by the nozzle's area ratio. For fixed geometry and gas properties, the Mach-number profile is determined by the area ratio. Changing $T_0$ changes the temperature and speed-of-sound fields, and therefore scales the resulting velocity field.

For the same nozzle geometry and gas properties, the Mach number at every point is fixed by the area ratio. Temperature does not change that dimensionless profile.

Velocity, however, scales with the local speed of sound:

u = Ma = M\sqrt{\gamma R T}

where:

$u$ is the local flow velocity;
$M$ is the local Mach number;
$a$ is the local speed of sound;
$\gamma$ is the ratio of specific heats of the gas;
$R$ is the specific gas constant;
$T$ is the local static temperature.

The equation relates flow velocity to Mach number and the local speed of sound. For an ideal gas, the speed of sound is

a = \sqrt{\gamma R T}

, so the flow velocity is

u = M\sqrt{\gamma R T}

and because the temperature field scales with $T_0$ ,

u \propto \sqrt{T_0}

where:

$u$ is the local flow velocity;
$T_0$ is the stagnation temperature.

The proportionality means that, for otherwise comparable flow conditions, the characteristic flow velocity increases with the square root of the stagnation temperature. Doubling

T_0

therefore increases

u

by a factor of

\sqrt{2}

Geometry determines the dimensionless flow structure, while thermodynamic conditions determine the velocity scale.

There is an important simplification here. In a real rocket, you do not turn a T₀ knob directly. You choose a propellant combination and mixture ratio, and those determine the combustion products, temperature, molecular weight, γ, and therefore the performance.

So T₀ is a convenient stand-in for combustion, not a complete combustion model.

Limitations

Because the flow is quasi-1D, there is no wall boundary layer and no radial profile; the model has no dimension to put them in. A real nozzle has viscous effects near the walls, and the velocity there is not the same as it is along the centreline.

Because it is isentropic, it contains no shocks. In a rocket plume you may see "Mach diamonds", shock-cell structures produced as the exhaust adjusts to the surrounding pressure. This model stops at the exit plane, so it does not capture them.

The flow field is solved from an idealized model, and the particles are simulated from that field. The vertical spread of the particles is decorative, because quasi-1D flow gives us no radial velocity field to solve.

There is one more simplification worth keeping in mind. γ and R are treated as constants. Real combustion products have temperature-dependent thermodynamic properties and can change composition through the nozzle. A higher-fidelity model would need to account for those effects.

Closing Thoughts

Below the speed of sound, gas accelerates when the channel narrows whereas above the speed of sound, gas accelerates when the channel widens. This transition is made possible by the throat, a critical point where the mach number is exactly 1 (equal to the speed of sound).

To validate the solution we confirm that the throat approaches Mach 1, the independently calculated exit velocity agrees with the solver, and mass flux stays effectively constant across the nozzle. Within the assumptions of the quasi-one-dimensional, isentropic model, the picture and the numbers agree.

That is the whole program, and it fits in about 150 lines.

I hope you found this read enjoyable and valuable. I invite you to share your own thoughts and ideas through the comments below.

You can find the code for this piece on my GitHub

You might assume that since the converging section narrows, the gas must be compressed. However, this is not the case. In the idealized nozzle, the converging section accelerates the gas, converting enthalpy into kinetic energy, while its temperature and pressure fall. The chamber is where the gas is pressurized and heated before entering the nozzle. ↩
The chamber values T₀ and P₀ are stagnation conditions, the temperature and pressure the gas would have if brought to rest isentropically. The assumption is that the combustion chamber is a large reservoir where the gas is effectively motionless, so its conditions are the stagnation conditions. The nozzle inlet is not quite at rest, though: by the time the gas reaches the start of the converging section it already carries a small subsonic velocity. True zero velocity lives in the reservoir upstream, not at the inlet plane. ↩
The two branches are the one genuine subtlety. For any area ratio greater than 1, the area-Mach relation has two valid solutions, one subsonic and one supersonic. The converging section takes the subsonic root, the diverging section the supersonic root. Picking the wrong branch is the classic way this solver can produce a smooth-looking but physically wrong result. ↩
T₀ is a useful stand-in for combustion, but it is not combustion itself. In a real rocket, the propellant combination and mixture ratio determine the combustion products and their thermodynamic properties. Those, in turn, determine the chamber temperature and the gas constants used by the nozzle. ↩