Recently, I’ve been working on a playground for the Glenn.jl package. The goal is to create a repository of ideas and usage examples for the package to inspire users. The repository will also serve as a way to aggregate ideas from everyone who wants and is able to contribute.
Feel free to fork it, work on brilliant ideas, and submit a pull request when you’re ready — just don’t forget to place your project in a separate directory and document it very well.
Below I’m leaving the simulation report of a KNSB Motor.
Playground: GitHub - ProfLeao/PlaygroundGlenn.jl · GitHub
Glenn.jl: GitHub - ProfLeao/Glenn.jl: A Julia Lang package to access the NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species · GitHub
KNSB Motor Simulation Report
Original run: 2026-08-11 11:20:00 · Revised: 2026-08-15
Simulator: MiniRocket.jl — Solid Rocket Combustion CFD
Model: 0D internal ballistics (uniform chamber)
NAR/Tripoli Classification: Class B (2.51–5.00 N·s) — original run · Class I after burn-rate correction
1. Physical Model
1.1 Governing Equations (0D Model)
The zero-dimensional model treats the entire combustion chamber as a single uniform control volume. Three coupled equations govern the system:
Conservation of mass in the chamber:
\frac{d}{dt}(\rho_c V_c) = \underbrace{\rho_p A_b r_b}_{\text{gas generation}} - \underbrace{\dot{m}_{nozzle}}_{\text{nozzle exit}}
where \rho_c is the chamber gas density, V_c the chamber free volume, \rho_p the propellant density, A_b the burning surface area, r_b the linear burning rate, and \dot{m}_{nozzle} the mass flow through the nozzle.
Saint Robert’s burning rate law:
r_b = a \cdot P_c^{\,n}
This empirical power-law relates the linear regression rate of the propellant surface to chamber pressure P_c. The coefficients a (burning rate coefficient) and n (pressure exponent, typically 0.3–0.5 for sugar propellants) are determined experimentally. The exponent n < 1 is critical for stability — if n \geq 1, the motor is inherently unstable and prone to catastrophic overpressure.
Equation of state (ideal gas):
P_c = \frac{m_{gas} \cdot R_{gas} \cdot T_c}{V_c}
The chamber gas is treated as an ideal gas at the adiabatic flame temperature T_c, with gas constant R_{gas} = R_u / MW_{products}. This assumption is reasonable for the low-to-moderate pressures (~1–20 bar) typical of amateur solid motors.
Nozzle mass flow (choked):
\dot{m}_{nozzle} = \frac{P_c \cdot A_t}{c^*}
where A_t is the throat area and c^* is the characteristic velocity:
c^* = \frac{\sqrt{\gamma \cdot R_{gas} \cdot T_c}}{\gamma \cdot \left(\frac{2}{\gamma+1}\right)^{\frac{\gamma+1}{2(\gamma-1)}}}
Choked flow (M = 1 at the throat) is assumed whenever P_c / P_{amb} \gtrsim 1.8, which holds for all meaningful operating conditions.
Thrust equation:
F = \eta \cdot \left[\dot{m}_{nozzle} \cdot v_e + (P_e - P_{amb}) \cdot A_e\right]
where the exit velocity v_e and exit pressure P_e are computed via isentropic expansion from the throat to the exit area ratio A_e/A_t, and \eta accounts for nozzle losses (friction, divergence, two-phase flow).
1.2 The Kn Parameter — Key to Motor Design
The dimensionless parameter Kn = A_b / A_t (burning area to throat area ratio) is the single most important design variable for solid rocket motors:
- The equilibrium chamber pressure is approximately P_c \propto Kn^{1/(1-n)}
- Higher Kn → higher pressure → higher thrust → faster burn
- Lower Kn → lower pressure → risk of extinction
- Typical amateur KNSB motors operate at $Kn \approx 200$–$450$ (Nakka), with \approx 300 a common design point
For this motor, the initial Kn = 250. Because the grain is BATES with uninhibited ends, the burning area changes over time. For this geometry (L/D = 3.75, web = 15 mm) the area is progressive: A_b grows from 70.7 cm² toward ~150 cm² at burnout, so Kn increases during the burn (250 → ~360 within the simulated 10 s) and the equilibrium pressure slowly rises with it.
2. Motor Configuration
2.1 Propellant: KNSB
| Property |
Value |
Description |
| Name |
KNSB (65% KNO₃ + 35% C₆H₁₄O₆) |
Potassium Nitrate + Sorbitol |
| Density |
1750 kg/m³ |
Compacted granular propellant |
| T_adiabatic (estimated) |
1600 K |
From literature (Nakka) |
| γ (estimated) |
1.220 |
Specific heat ratio of product gases |
| MW products |
28.5 g/mol |
Average molecular weight of exhaust |
| Heat of explosion |
3.8 MJ/kg |
Chemical energy content |
| Burning rate: a |
8.26×10⁻⁶ m/s/Paⁿ |
Saint Robert coefficient |
| Burning rate: n |
0.319 |
Pressure exponent (< 1 → stable) |
KNSB is the most widely used propellant in Brazilian model rocketry due to its low cost, stable burning characteristics, and ease of manufacture. The combustion reaction (simplified) is:
10\,\text{KNO}_3 + 3\,\text{C}_6\text{H}_{14}\text{O}_6 \rightarrow 5\,\text{K}_2\text{CO}_3 + 6\,\text{CO}_2 + 7\,\text{CO} + 14\,\text{H}_2\text{O} + 7\,\text{H}_2 + 5\,\text{N}_2
This is a balanced, representative form (10 : 3 molar ratio = 65.0 : 35.0 by weight, matching the 65/35 KNSB mix). Because KNSB is fuel-rich, the products contain substantial CO and H₂ in addition to CO₂; the exact CO₂/CO/H₂/H₂O split is set by chemical equilibrium (the water–gas shift), so this equation is representative rather than unique. The original report’s equation (5 KNO₃ + C₆H₁₄O₆ → …) was not atom-balanced — K, N, C and O did not match.
2.2 Grain Geometry: BATES
| Property |
Value |
Significance |
| Length |
150 mm |
Total grain length |
| Outer diameter |
40 mm |
Fits standard 40 mm PVC casing |
| Inner diameter |
10 mm |
Initial port for flame propagation |
| Outer inhibited |
Yes |
PVC casing prevents outer surface burn |
| Ends inhibited |
No |
Both ends burn → progressive thrust profile |
| Initial burning area |
70.7 cm² |
Sum of inner cylinder + both end faces |
| Propellant volume |
176.7 cm³ |
Total solid propellant |
| Propellant mass |
309.3 g |
≈ 0.3 kg fuel load |
The BATES (Ballistic Test and Evaluation System) geometry is a hollow cylinder that burns on all non-inhibited surfaces. For this grain the inner cylindrical surface dominates: as the port enlarges (r_i: 5 → 20 mm) the inner area grows faster than the end faces shrink, so the total burning area is progressive (70.7 → ~150 cm² at burnout), not regressive. The original report described this geometry as regressive, which is incorrect — regressive behaviour occurs only in short, end-face-dominated grains.
2.3 Motor Assembly
| Property |
Value |
Significance |
| Chamber (L × D) |
180 × 40 mm |
30 mm extra length for combustion zone |
| Throat diameter (D_t) |
6.0 mm |
Critical flow control |
| Exit diameter (D_e) |
12.0 mm |
4:1 area expansion ratio |
| Nozzle efficiency |
92% |
Accounts for losses |
| Initial Kn (A_b/A_t) |
250 |
Design point for pressure balance |
The nozzle is the convergent-divergent type (de Laval). At 6.0 mm throat, A_t = 28.3\text{ mm}^2 and Kn = 250 — within the typical KNSB operating range, so the throat is not abnormally large. The real mismatch is the expansion ratio: \varepsilon = 4 (D_e = 12 mm) is matched to a chamber pressure of ~20 bar, whereas this motor only reaches ~3 bar, so the nozzle is severely over-expanded (see §3.2).
3. Simulation Results
3.1 Performance Summary
| Metric |
Value |
Notes |
| Burn time |
10.000 s |
Truncated at T_max — only ~21% of propellant consumed |
| Max chamber pressure |
3.08 bar |
Still slowly rising at cutoff (no plateau) |
| Avg chamber pressure |
2.39 bar |
Corrected (original said ~0.91 bar) |
| Max thrust |
2.28 N |
Occurs at t = 10 s (cutoff), not at ignition |
| Avg thrust |
0.46 N |
≈ 0 for the first ~6 s (over-expanded nozzle) |
| Total impulse |
4.64 N·s |
Low for a 309 g propellant load |
| Specific impulse (Isp) |
74.9 s |
≈ 76% of the ideal ~99 s at 3 bar |
| Propellant consumed |
65 g (21%) |
244 g remains |
| NAR/Tripoli Class |
B |
2.51–5.00 N·s (corrected; original said 1/2A) |
3.2 Thrust and Pressure Curves
The combined plot reveals several key features (several differ from the original report):
- Gradual pressurization, no sharp ignition spike — pressure rises from 1.01 bar to ~1.8 bar over the first ~1–2 s, then slowly climbs to 3.08 bar, tracking the increasing Kn. The model starts at ambient pressure with a 50% “incipient” burn rate below 1.5\,P_{amb}, so no explosive ignition transient is captured.
- No pressure plateau — at t = 10 s the pressure is still rising; it never levels off because the equilibrium pressure itself rises as the port grows.
- Thrust ≈ 0 for the first ~6 s despite rising pressure — thrust does not simply follow pressure here.
- Over-expanded nozzle is the cause: at \varepsilon = 4 and \gamma = 1.22, M_e \approx 2.65 and P_e/P_c \approx 0.042. At P_c = 3.08 bar this gives P_e \approx 0.13 bar, so the pressure term (P_e - P_{amb})\,A_e \approx -10 N nearly cancels the \approx 12.5 N momentum term. Thrust only turns positive once P_c \gtrsim 2.5 bar. The ideal expansion ratio at 3 bar would be \varepsilon \approx 1.14; the 4:1 nozzle is sized for ~20 bar.
3.3 Grain Evolution
The grain evolution plots show (corrected values from knsb_results.csv):
- Burning area: increases from 70.7 to 102.2 cm² (+44%) over the 10 s — the grain is progressive, not regressive.
- Burned web: 4.28 mm of the 15 mm web (28.5%). The burn rate is ~0.35–0.46 mm/s (r_b \propto P_c^{0.319} at 1–3 bar).
- Remaining mass: falls from 309 g to 244 g — 65 g (21%) consumed, not “nearly flat”.
3.4 Why Did This Motor Underperform? (Root-Cause Analysis)
The original report attributed the failure to “oversized throat / low Kn”. The review below shows the dominant cause is a ~12× error in the burning-rate coefficient, not the throat.
1. Dominant error — the Saint Robert coefficient is ~12× too low. The model uses a = 8.26\times10^{-6}\ \text{m/(s·Pa}^n), which yields a burn rate of only \approx 0.68 mm/s at 10 bar and \approx 1.25 mm/s at 68 bar. The widely cited Nakka law for KNSB is r[\text{mm/s}] = 8.26\cdot P[\text{MPa}]^{0.319} — i.e. \approx 8.3 mm/s at 10 bar and \approx 15 mm/s at 68 bar, ~12× faster. This is consistent with a unit-conversion error (mm/s per MPaⁿ vs m/s per Paⁿ); the correct SI coefficient is a \approx 1.0\times10^{-4}\ \text{m/(s·Pa}^n), not 8.26\times10^{-6}.
2. Consequence. The equilibrium pressure P_{eq} = \left(\rho_p a c^* Kn\right)^{1/(1-n)} scales as a^{1.47}. With the wrong a the motor settles near ~1.8 bar (rising with Kn to ~3 bar). After correcting a and re-running, this same 40 mm motor reaches ~200 bar peak, 822 N and Class I (466.8 N·s) — far beyond what a 40 mm PVC casing can withstand. The 6 mm throat is not the problem; with the correct burn rate this throat is already aggressive for such a casing.
3. Secondary error — over-expanded nozzle. Even at the low pressures reached, the 4:1 nozzle cancels ~80% of the momentum thrust (see §3.2), which is why thrust ≈ 0 for the first 6 s and Isp is reduced.
4. The parametric study’s “fix” is unsafe. Its recommendation (reduce D_t to 4–5 mm) is an artifact of the wrong a. With a corrected coefficient, D_t = 4 mm (Kn \approx 700) would imply hundreds of bar — a CATO (catastrophic failure) risk. The parametric study has since been re-run with the corrected coefficient — see parametric_report.md: all 25 configurations are Class I, 85–601 bar, confirming this risk.
4. Glenn.jl Thermochemical Verification
Glenn.jl provides high-fidelity thermochemical properties from the NASA-7 polynomial database (~2030 species). This section cross-validates the propellant’s estimated properties against Glenn.jl calculations.
4.1 Product Species at Adiabatic Temperature
| Species |
Cp at 1600 K (J/mol·K) |
S° at 1600 K (J/mol·K) |
Phase |
| CO₂ |
58.87 |
295.98 |
Gas |
| CO |
35.47 |
250.71 |
Gas |
| H₂O |
48.34 |
253.74 |
Gas |
| N₂ |
35.13 |
244.14 |
Gas |
| H₂ |
32.73 |
180.94 |
Gas |
| CH₄ |
93.37 |
287.69 |
Gas |
The NASA-7 polynomials capture the temperature dependence of Cp(T) with high accuracy (typically < 1% error in the 200–6000 K range). The large Cp of CH₄ (93.37 J/mol·K) reflects its polyatomic structure with many vibrational degrees of freedom, while the diatomic species (CO, N₂, H₂) cluster around ~33–35 J/mol·K — close to the theoretical \frac{7}{2}R = 29.1 J/mol·K for diatomic ideal gases.
4.2 Adiabatic Temperature and Impulse Comparison
| Comparison |
Value |
Source |
| T_adiabatic (Glenn.jl) |
4118 K |
NASA-7 + simplified equilibrium |
| T_adiabatic (estimated) |
1600 K |
Literature (Nakka) |
| Isp theoretical (Glenn.jl + isentropic) |
141.0 s |
Frozen isentropic expansion |
| Isp simulated (0D) |
74.9 s |
0D ballistics at ~3 bar |
Critical observation: The Glenn.jl adiabatic temperature (4118 K) overpredicts the literature value (~1600 K). The root causes are more specific than “simplified equilibrium”:
- Zero reactants enthalpy.
knsb_rocket.jl calls adiabatic_flame_temperature(calc, 0.0, products_glenn). Passing 0.0 for the reactants enthalpy makes the energy balance H_{products}(T) = 0, which searches for the temperature where the products’ absolute enthalpy (including their strongly negative heats of formation) crosses zero. This is not a meaningful flame-temperature condition — the heat of explosion (3.8 MJ/kg) is never supplied to the calculation.
- Condensed species treated as gas. K₂CO₃ and KOH are condensed (solid/liquid) products;
estimate_combustion_products deliberately selects their gas-phase entries, removing their latent heat of vaporization and inflating T.
- No Gibbs minimization. The product set is a fixed guess; a full equilibrium calculation (NASA CEA) would redistribute species and lower T.
Isp consequence. The “theoretical” 141.0 s is inflated by the bogus T_{ad} (I_{sp} \propto \sqrt{T_c}). Using the correct T_c = 1600 K, the ideal sea-level frozen Isp at P_c = 3.08 bar is \approx 99 s; the simulated 74.9 s is \approx 76\% of that (consistent with \eta = 0.92 plus over-expansion). The original claim that “141 s aligns with literature” was a coincidence of two compensating errors. After the energy-balance correction, the re-run gives T_{ad} \approx 2602 K and Isp$_{theo} \approx 224$ s, with simulated Isp \approx 154 s (\approx 68\% of theoretical).
Recommendation: pass the heat of explosion as the reactants enthalpy, treat K₂CO₃/KOH as condensed, and cross-validate against NASA CEA.
5. Conclusions
| Aspect |
Finding |
| Motor Class |
B (2.51–5.00 N·s) — corrected from “1/2A” |
| Primary Issue |
Burning-rate coefficient ~12× too low (unit-conversion error) |
| Secondary Issue |
4:1 nozzle over-expanded at ~3 bar (ideal \varepsilon \approx 1.14) |
| Grain behaviour |
Progressive (A_b: 70.7 → 102 cm² in 10 s), not regressive |
| Glenn.jl Validation |
Species Cp/S° accurate; T_{ad} = 4118 K is an artifact of zero reactants enthalpy + gas-phase condensed species |
| Fix |
Correct a \approx 1.0\times10^{-4} m/(s·Paⁿ) and re-run before trusting any throat recommendation |
| Model Limitations |
0D, fixed T_c, no heat losses, no ignition transient |
6. Model Issues Identified During Review
| # |
Issue |
Severity |
Location |
Suggested Fix |
| 1 |
Saint Robert coefficient a ~12× too low (unit error) |
Critical |
src/propellant.jl (KNSB) |
a \approx 1.0\times10^{-4} m/(s·Paⁿ) |
| 2 |
NAR/Tripoli class thresholds deviate from the standard table (4.64 N·s is reported as “1/2A” but is Class B) |
High |
src/ballistics.jl (print_results), README.md |
Use the standard table (B = 2.51–5.00 N·s) |
| 3 |
adiabatic_flame_temperature called with zero reactants enthalpy |
High |
examples/knsb_rocket.jl |
Pass the heat of explosion (3.8 MJ/kg) |
| 4 |
Condensed products (K₂CO₃, KOH) looked up as gas phase |
Medium |
src/equilibrium.jl |
Keep condensed-phase entries |
| 5 |
No ignition transient / fixed T_c / no heat losses |
Medium |
src/ballistics.jl |
Add ignition model, T_c(P), or quasi-1D |
| 6 |
Fixed 4:1 expansion assumed appropriate at all pressures |
Low |
examples/*.jl |
Size \varepsilon to the design P_c |
Safety note. Before flying or static-testing any configuration, re-run the simulation with the corrected burning-rate coefficient. With the original (too-low) coefficient the tool systematically underpredicts pressure and thrust: after the correction, this same 40 mm motor re-runs at ~200 bar / 822 N / Class I (vs. the original 3 bar / 2.3 N / Class B), confirming that throat sizes optimized against the wrong coefficient (e.g., D_t → 4–5 mm) would be catastrophic.
Report based on execution of examples/knsb_rocket.jl using MiniRocket.jl; revised 2026-08-15 after a full review of the code, the raw results (examples/data/knsb_results.csv) and the execution logs. After correcting the burning-rate coefficient and re-running (2026-08-15), the motor predicts ~200 bar / Class I — see §3.4. The parametric study (parametric_report.md) has also been re-run with the corrected coefficient.
