Building a Supersonic Nozzle CFD Simulator in Python

  

Introduction

When designing propulsion nozzles for rocket engines or supersonic aircraft, the central challenge is straightforward in concept yet demanding in execution: accelerate hot, high-pressure combustion gases as efficiently as possible to maximize thrust. This article walks through the complete process of reimplementing a MATLAB-based supersonic nozzle CFD code in Python — from nozzle geometry design and mesh generation through to numerical solution of the two-dimensional Euler equations using the MacCormack finite volume method. Along the way, the theoretical background behind each step is explained alongside the implementation details.

---

1. Nozzle Geometry Design — Method of Characteristics

The starting point of nozzle design is finding the wall geometry of a minimum-length nozzle that achieves the desired exit Mach number. This is accomplished using the Method of Characteristics (MOC), a technique that transforms partial differential equations into ordinary differential equations integrated along characteristic lines, accurately capturing wave propagation in two-dimensional supersonic flow.

 
The approach is grounded in Prandtl-Meyer expansion theory. Characteristic lines originating at the throat reflect between the axis of symmetry and the nozzle wall, and the coordinates of each intersection point are computed sequentially. At every point, the Riemann invariants $K^+ = \theta + \nu$ and $K^- = \theta - \nu$ are conserved, from which the flow angle $\theta$ and the Prandtl-Meyer function $\nu$ are determined. Extending each characteristic line to the wall then yields the nozzle wall coordinates.
 

In the Python implementation, the generate_bell_moc() function in nozzle_shape.py handles this process. The Prandtl-Meyer function is inverted at each point using Newton-Raphson iteration to recover the local Mach number, and 2×2 linear systems are solved to find characteristic intersection coordinates. Increasing the number of characteristic lines num improves design accuracy; with the default value of 15, sixteen wall coordinate points are produced.

---

2. Mesh Generation — Structured Grid

Once the nozzle geometry is established, a computational grid must be built for the CFD analysis. The build_mesh() function in noz_mesh.py generates a body-fitted structured mesh.

 

Because the MOC wall coordinates consist of only 16 unevenly spaced points, the wall curve is first resampled to 300 points using scipy's CubicSpline, producing a smooth continuous representation. The x-direction is then divided into uniform intervals, and in each x-column the y-direction is distributed linearly from the axis of symmetry (y=0) to the wall. Grid resolution is controlled by the mesh_factor parameter, ranging from 0.05 (coarse) to 0.20 (fine). At medium resolution (mesh_factor=0.10), the resulting mesh has approximately 91×72 nodes.

---

3. Governing Equations — 2D Euler Equations

The simulator solves the two-dimensional Euler equations, which neglect viscous effects. In conservation form they read:
$$\frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = 0$$
where the conserved variable vector $\mathbf{Q}$, the x-direction flux vector $\mathbf{F}$, and the y-direction flux vector $\mathbf{G}$ are defined as:
$$\mathbf{Q} = \begin{pmatrix} \rho \ \rho u \ \rho v \ e \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho u \ \rho u^2 + p \ \rho u v \ (e+p)u \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \rho v \ \rho u v \ \rho v^2 + p \ (e+p)v \end{pmatrix}$$
Pressure is recovered from the perfect gas equation of state $p = (\gamma - 1)\left(e - \frac{1}{2}\rho(u^2+v^2)\right)$, and a specific heat ratio of $\gamma = 1.25$ is used to reflect combustion gas properties.

---

4. Numerical Method — MacCormack Finite Volume Scheme

Time integration of the Euler equations uses the MacCormack scheme, a two-stage predictor-corrector method that achieves second-order accuracy in both space and time.
 
In the predictor step (dd=-1), a forward-difference approximation produces a provisional solution $\bar{\mathbf{Q}}$:
$$\bar{\mathbf{Q}}^n = \mathbf{Q}^n - \Delta t \cdot \mathcal{L}^n(\mathbf{Q}^n)$$
 
In the corrector step (dd=0), a backward-difference approximation refines the result:
$$\mathbf{Q}^{n+1} = \frac{1}{2}\left(\mathbf{Q}^n + \bar{\mathbf{Q}}^n - \Delta t \cdot \mathcal{L}(\bar{\mathbf{Q}}^n)\right)$$
 
Because the finite volume method is used, the flux divergence for each cell is assembled from contributions at its four faces — right, left, top, and bottom. Face normal vectors are computed geometrically from the grid coordinates, enabling accurate flux evaluation on non-orthogonal grids.
 
A particularly important aspect of the Python implementation is the index translation from MATLAB (1-based) to Python (0-based). The MATLAB loop variable $I = 2\ldots n_x$ (1-based) maps to Python's $i = 1\ldots n_x-1$ (0-based), and the flux array references translate as follows:

MATLAB (1-based)	Python (0-based)	dd=-1	dd=0
F(I+1+dd, J)	F[i+dd, j]	F[i-1, j]	F[i+1, j]
F(I+dd, J)	F[i+dd-1, j]	F[i, j]	F[i, j]
F(I, J+1+dd)	F[i, j+dd]	F[i, j-1]	F[i, j+1]
F(I, J+dd)	F[i, j+dd-1]	F[i, j]	F[i, j]

The initial implementation faithfully followed the MATLAB double for-loop structure. After replacing it with NumPy array slicing, a speedup of approximately 35× was achieved. The vectorized predictor-step flux computation, for example, reduces to the following compact form:

# dd=-1 (forward predictor) — vectorised
Fc = F[1:nx, 1:ny, :]    # F[i, j]   — right/top face
Fl = F[0:nx-1, 1:ny, :]  # F[i-1, j] — left face
Fb = F[1:nx, 0:ny-1, :]  # F[i, j-1] — bottom face

div = (Fc*b(sfpx) + Gc*b(sfpy) +   # right face
       Fl*b(sfmx) + Gl*b(sfmy) +   # left face
       Fc*b(sgpx) + Gc*b(sgpy) +   # top face  (Fc = F[i,j])
       Fb*b(sgmx) + Gb*b(sgmy))    # bottom face

---

5. Boundary Conditions

Correct boundary condition specification is directly tied to numerical stability. The conditions applied at each boundary are as follows.

• Inflow (i=0): A sonic inflow condition is imposed at the nozzle throat. The axial velocity is fixed at the sonic speed $u = \sqrt{\gamma R T_c}$, while the transverse velocity is extrapolated from the adjacent interior cell.
• Outflow (i=nx): Since the exit flow is fully supersonic, all flow variables are zero-order extrapolated from the interior.
• Axis of symmetry (j=0): Pressure, density, and axial velocity are mirrored from the interior cell, while the transverse velocity is sign-reversed to enforce a reflection boundary.
• Wall (j=ny): A slip-wall condition consistent with the Euler equations is applied. Pressure and density are extrapolated from the interior, and both velocity components are set to zero.

---

6. Time Integration and Stability

The time step $\Delta t$ is determined automatically at each iteration by the Courant-Friedrichs-Lewy (CFL) condition:
$$\Delta t = \frac{\text{CFL}}{\max_{i,j}\left(\frac{a_{ij} + b_{ij} + c_{ij}}{V_{ij}}\right)}$$
 
Here $a_{ij}$ and $b_{ij}$ represent the velocity contributions in each direction, $c_{ij}$ is the acoustic term, and $V_{ij}$ is the cell area computed via the shoelace formula. A default CFL number of 0.8 is used, ensuring that disturbances travel no more than one cell per time step.

---

7. Graphical User Interface

The user interface is built with Python's built-in Tkinter library combined with Matplotlib for rendering. Tkinter was chosen over PyQt5 because it requires no additional installation and runs out of the box in any standard Python environment.

 

The GUI is divided into a left-hand input panel and a right-hand results area with tabbed views. The input panel sits inside a scrollable canvas and contains controls for nozzle type selection (Cone or Bell), geometric parameters, flow conditions, and CFD settings. The Preview and Run CFD buttons are pinned to the bottom of the left column so they remain visible at all times. CFD computation runs in a background threading.Thread, keeping the interface responsive throughout the calculation.

 
Results are displayed across five tabs. Mesh Preview shows the structured grid and nozzle wall geometry. Mach Contour renders the two-dimensional Mach number distribution using the RdYlGn colormap. Axial Profiles plots the Mach number and pressure ratio along the axis of symmetry. The Characteristics tab, available for Bell nozzles, visualizes the full MOC characteristic line network and the wall intersection points used during the design process. Finally, the Summary tab presents a structured table comparing nozzle geometry, CFD results, and isentropic reference values side by side.

---

8. Validation

Running the solver under the same conditions as the original MATLAB code — throat height h_th=0.025 m, chamber pressure P_c=1.2 MPa, chamber temperature T_c=2000 K, γ=1.25, and design exit Mach number M_e=2.26 — produces a Mach number range of approximately 1.2 to 2.4, which agrees well with the MATLAB reference of 1.0 to 2.45. The minimum Mach number near the inlet does not reach exactly 1.0 because the original MATLAB boundary condition fixes only the velocity at the throat while extrapolating pressure and density from the interior; this behavior is intentionally preserved in the Python translation.

 
In terms of performance, NumPy vectorization delivers roughly a 35× speedup over the for-loop version, completing 500 time steps on a 91×72 mesh in approximately two seconds.

---

Conclusion

This project was guided throughout by a single principle: preserve the physical and numerical structure of the original MATLAB code while taking advantage of Python's ecosystem. Three core techniques — the Method of Characteristics, the MacCormack scheme, and the finite volume method — work together to simulate supersonic nozzle flow from geometry design all the way through to a converged flow field. Looking ahead, adding artificial viscosity to suppress numerical oscillations and replacing the inlet boundary condition with a fully isentropic specification would further improve solution accuracy and robustness.