The complete description of every model inside Project APEX — stability, drag, supersonic aerodynamics, numerical integration, and atmosphere.
Project APEX solves a one-dimensional vertical flight from launch to apogee using a physics-based model with no lookup tables for the core dynamics. Every quantity — drag coefficient, stability margin, air density, speed of sound — is computed analytically at each timestep from the current state of the rocket and the atmosphere.
The model has three main layers. The aerodynamics layer computes the centre of pressure (CP), the drag coefficient (Cd), and the lift-curve slope from the rocket's geometry. The flight integration layer uses those forces together with the motor thrust curve to propagate position and velocity forward in time. The atmosphere layer supplies air density, dynamic viscosity, and speed of sound at each altitude, updated at every step.
All internal calculations use SI units throughout. Inputs entered in imperial units are converted before the simulation runs and converted back for display.
| Layer | Output | Method |
|---|---|---|
| Stability | Xcp, CNα, static margin | Barrowman + Rogers Kbf (1966/2011) + DATCOM 3D Jones supersonic (M > 1.2) |
| Drag | Total Cd (7 components) | Analytical, Re-based skin friction, power-on base drag correction |
| Compressibility | Fin CNα Mach correction | PG subsonic (M < 0.8), DATCOM 3D Jones supersonic (M ≥ 1.2), linear blend transonic |
| Integration | Altitude, velocity vs time | 4th-order Runge-Kutta, dt = 0.05 s |
| Atmosphere | ρ, T, P, a, μ | ISA 1976, 3 layers to 32 km |
The centre of pressure is calculated using the Barrowman method, which treats the rocket as a body of revolution and sums the normal-force contributions of the nosecone and fins separately. The CP is the weighted average of those contributions, weighted by their respective normal-force slope coefficients (CNα).
Xcp = (CNα_nose · Xcp_nose + CNα_fins · Xcp_fins) / (CNα_nose + CNα_fins)
All distances are measured from the nose tip. A larger Xcp means the CP is further aft, which increases the stability margin.
For a slender body of revolution, the nosecone normal-force slope is 2.0 regardless of shape. The shape affects the CP location of the nose contribution:
| Nose shape | Xcp_nose (from tip) |
|---|---|
| Cone | 2/3 · Ln |
| Ogive | 0.466 · Ln |
| Von Kármán | 0.437 · Ln |
| Parabolic | 0.500 · Ln |
APEX uses the Rogers Modified Barrowman method, which corrects two omissions in Barrowman's original 1966 formulation. Standard Barrowman accounts for how the body changes fin lift (factor Kfb), but omits the lift induced on the body surface adjacent to the fin roots by the fin's shed vortex (factor Kbf). Both are included:
Kfb = 1 + r/(s+r) // fins-in-presence-of-body (standard Barrowman) Kbf = Kfb · (r/(s+r))² // body-in-presence-of-fins (Rogers addition) K = Kfb + Kbf // combined interference factor used in APEX lm = √(s² + sb²) // mid-chord sweep length sb = sweepLen + (CT − CR)/2 // Barrowman mid-chord sweep distance CNα_fins = K · (4·N·(s/d)²) / (1 + √(1 + (2·lm / (CR + CT))²))
where N = number of fins, s = semi-span (body surface to tip), d = body diameter, CR and CT = root and tip chord. For typical HPR geometry (r/(s+r) ≈ 0.33), the Kbf body-carryover term adds approximately 11% to fin CNα compared to standard Barrowman, shifting CP slightly aft and increasing stability margin by 0.3–0.5 calibers.
Xcp_fins = Xf + (lm/3)·(CR + 2·CT)/(CR + CT) + (1/6)·(CR + CT − CR·CT/(CR+CT))
where Xf is the distance from nose tip to the fin root leading edge.
APEX recomputes CP at every recorded timestep as Mach changes through the flight. The correction to fin normal-force slope uses three regimes:
M < 0.8 : CNα_fin_eff = CNα_fin · 1/√(1−M²) // Prandtl-Glauert subsonic boost 0.8–1.2 : linear blend between PG and DATCOM anchors M ≥ 1.2 : CNα_fin_sup = K · N · [4·AR/(AR·β+2)] · (A_fin/A_ref) // DATCOM 3D Jones AR = s²/A_fin (panel aspect ratio), β = √(M²−1), A_fin = single-fin planform area
Why DATCOM rather than the simpler Ackeret 1/√(M²−1) multiplier? The Ackeret formula is a 2D thin-wing result. Applying it as a multiplier on Barrowman's already-3D fin CNα double-penalises fin authority at high Mach. At M = 2.2, the 2D approach models fins at only 51% of subsonic effectiveness, pushing CP far forward and producing an unrealistically low stability margin. The DATCOM 3D Jones formula correctly accounts for finite span: for low-AR fins (AR ≈ 0.5–1.5, typical in HPR) it retains 77–90% of subsonic effectiveness at M = 2.2 — consistent with RASAero's DATCOM wind-tunnel-validated results at M = 1.5–3.
Transonic stability margin dip. As the rocket passes through M = 0.8–1.2, fin CNα transitions from the PG-boosted value (KPG ≈ 1.67 at M = 0.8, fins 67% stronger) to the DATCOM supersonic value (fins at ~85–90% of subsonic effectiveness at M = 2+). This creates a real, physical dip in the Stability chart through the transonic band. For a well-designed HPR rocket the minimum margin stays above 2 calibers.
The static margin is expressed in calibers (body diameters):
Static Margin = (Xcp − Xcg) / d
A margin of 1–2 calibers is the conventional HPR target. APEX applies an increasing drag penalty for unstable rockets — ramping from 1× at 0 calibers to up to 10× at −2 calibers — to reflect the tumbling behaviour of an unstable flight.
During the burn, the CG moves forward as propellant is consumed from the motor location:
Xcg(t) = (m_launch · Xcg0 − m_burned(t) · Xcg_motor) / m_current(t)
This means stability margin typically increases during the burn as the CG moves forward of its initial loaded position.
The total drag coefficient is the sum of seven analytically computed terms, each recalculated at every simulation step from the current Mach number and Reynolds number:
F_drag = 0.5 · ρ · v² · Cd_total · A_ref
where A_ref = π · r² (body cross-sectional area).
| Component | Symbol | Regime |
|---|---|---|
| Skin friction (body) | Cd_body | All speeds |
| Nose pressure drag | Cd_nose | All speeds |
| Base drag | Cd_base | All speeds (reduced during burn) |
| Fin friction + profile | Cd_fin_friction | All speeds |
| Fin LE/TE pressure | Cd_fin_LE/TE | All speeds |
| Fin interference | Cd_fin_int | All speeds |
| Fin wave drag | Cd_fin_wave | M > 0.75 |
| Body shoulder wave drag | Cd_wave | 0.82 < M < 1.25 |
Body skin friction uses the Prandtl-Schlichting turbulent flat-plate formula applied to the wetted body surface area:
Re = ρ · v · L / μ // Reynolds number (rocket length) Cf_smooth = 0.455 / (log₁₀(Re))²·⁵⁸ // Prandtl-Schlichting turbulent Cf_rough = 0.032 · (ε/L)^0.2 // Schlichting fully-rough (ε = roughness) Cf = max(Cf_smooth, Cf_rough) Cd_body = Cf · (1 + 60/(L/d)³ + 0.0025·(L/d)) · (π·d·L) / A_ref
The form factor term accounts for body fineness ratio. Surface roughness ε sets the lower limit — once fully rough, increasing Re no longer reduces Cf.
Fin skin friction is computed separately at the fin chord Reynolds number:
Re_fin = Re · (c_mac / L) // Re at mean aerodynamic chord Cf_fin = max(Cf_smooth_fin, Cf_rough_fin) // uses fin roughness ε_fin and c_mac Cd_fin_fric = Cf_fin · finFF · (2·N·A_fin) / A_ref
The Fin Roughness parameter is independent of body roughness — a painted carbon body (~20 µm) and polished aluminium fins (~0.5 µm) have very different Cf values and should be modelled separately.
The nosecone generates shape-dependent pressure drag that grows with Mach number. Below M = 0.8 it follows the Prandtl-Glauert correction; above M = 1.2 APEX uses shape-specific DATCOM coefficients anchored to the ogive via benchmark flight validation:
θ = atan(r / Ln) // nose half-angle (radians) β = √(M² − 1) // supersonic factor M ≤ 1.2: Cd_nose = Cd_nose_0(shape) · PG_correction(Mach) M > 1.2: Cd_nose = Kwave_n · θ² · √0.44 / β // DATCOM shape-specific wave drag
| Nose shape | Kwave_n | Relative to ogive | Physical reason |
|---|---|---|---|
| Cone | 3.0 | 1.43× more drag | Constant slope — steepest area gradient |
| Parabolic | 1.8 | 0.86× less drag | Gentler area distribution |
| Ogive | 2.1 | baseline | Anchored to M2.26 benchmark flight |
| Von Kármán | 1.4 | 0.67× less drag | LD-Haack series — minimum wave drag |
The √0.44 factor ensures continuity at M = 1.2. Switching from ogive to Von Kármán on a M2.2 rocket gains approximately 0.6% altitude; switching to a cone loses approximately 2.5%.
The rocket's blunt base creates a low-pressure wake. The formula is continuous at Mach 1:
M < 1.0: Cd_base = 0.025 + 0.165·M² // subsonic: 0.025 at rest → 0.190 at M=1 M ≥ 1.0: Cd_base = 0.19 / M^1.25 // supersonic: continuous at M=1, ~0.08 at M=2
Power-on plume correction. During motor burn, the exhaust plume fills the base wake and suppresses base drag. APEX applies this correction when a nozzle exit diameter is entered:
Cd_base_powered = Cd_base · max(1 − Ae/Ab, 0) Ae = π·(exitDia/2)² // nozzle exit area (wide end, not the throat) Ab = π·r² // rocket base area (= A_ref for flat-base rocket)
Set the Nozzle Exit Diameter input (Aerodynamics section, 0–100 mm) to activate this correction. With a typical Aerotech 75mm motor (exit ≈ 42 mm) in a 3" airframe, Ae/Ab ≈ 0.28 — base drag is reduced by 28% throughout the entire burn phase. Leave at zero to match OpenRocket's behaviour (no plume correction).
Fins contribute three terms beyond skin friction. The cross-section profile sets both the LE/TE pressure drag and the Ackeret wave drag coefficient:
| Cross-section | LE | TE | Kwave (supersonic) |
|---|---|---|---|
| Flat Plate | Blunt | Blunt | 4·t_c (linear) |
| Rounded LE | Rounded | Blunt | 4·t_c (linear) |
| LE Bevel (only) | Chamfered | Blunt | 2·t_c²/f_le |
| Quarter-Chord Double Wedge | Knife | Knife | 8·t_c² |
| Half-Chord Wedge (Diamond) | Knife | Knife | 4·t_c² |
| Airfoil (NACA symmetric) | Rounded | Cusped | 16/3·t_c² |
where t_c = fin thickness / root chord. Flat plate and rounded LE have Kwave ∝ t_c (not t_c²) because their blunt edges create a detached bow shock — wave drag scales with frontal area rather than deflection angle. All knife-edge profiles have Kwave ∝ t_c², which falls off much faster for thin fins.
Cd_fin_int = 2 · Cf · N · (t_fin · CR) / A_ref
Fin wave drag (Ackeret):
M_eff = M · cos(sweep_angle) // effective Mach at fin LE Cd_fin_wave = Kwave / √(max(M_eff²−1, 1.0)) // M_eff > 1.0 Cd_fin_wave = Kwave · sin(π·(M_eff−0.75)/0.5) // 0.75 < M_eff < 1.0 (build-up)
Swept fins maintain a lower M_eff than freestream Mach, delaying wave drag onset.
Body shoulder transonic wave drag — exists only in the transonic band where flow recompresses at the nosecone-to-body junction:
Cd_wave = 0.18 · sin(π·(M−0.82)/0.43) · (d/Ln)^0.4 // 0.82 < M < 1.25 only
Fin normal-force generation changes significantly with Mach number. APEX uses a three-regime model for the fin CNα Mach correction:
M < 0.8 : KPG = 1/√(1−M²) // Prandtl-Glauert — boosts fin CNα toward M=0.8 0.8–1.2 : linear blend // anchored at PG(0.8) ≈ 1.667 and DATCOM(1.2) M ≥ 1.2 : DATCOM 3D Jones formula // replaces the scalar PG multiplier entirely
Subsonic (M < 0.8): PG correction increases fin lift-curve slope. At M = 0.8, KPG ≈ 1.67 — fins are 67% more effective, CP is further aft, and stability margin is higher than at launch.
Supersonic (M ≥ 1.2): Previous APEX versions applied Ackeret's 2D result (1/√(M²−1)) as a multiplier on Barrowman's already-3D fin CNα. This double-penalised fin authority — at M = 2.2, fins dropped to only 51% of subsonic effectiveness and the stability margin chart showed an unrealistic dip to ~0.7 calibers. APEX now uses the DATCOM 3D Jones formula (see Stability section), which retains 77–90% of subsonic fin effectiveness at M = 2.2, producing a transonic margin dip of ~4 calibers — in agreement with RASAero's DATCOM-validated predictions.
Transonic blend (M = 0.8–1.2): Neither PG nor DATCOM is physically valid in the transonic band. APEX linearly blends from the M = 0.8 PG anchor (KPG ≈ 1.667) to the M = 1.2 DATCOM anchor. The transonic phase typically lasts less than 0.5 seconds for a typical HPR flight.
The flight is solved as a system of two first-order ODEs — altitude (dh/dt = v) and velocity (dv/dt = net acceleration) — using the fourth-order Runge-Kutta method (RK4) with a fixed 50 ms timestep.
k1 = deriv(h, v, t ) k2 = deriv(h+k1·dt/2, v+k1·dt/2, t+dt/2 ) k3 = deriv(h+k2·dt/2, v+k2·dt/2, t+dt/2 ) k4 = deriv(h+k3·dt, v+k3·dt, t+dt ) dh = (k1.dh + 2·k2.dh + 2·k3.dh + k4.dh) / 6 dv = (k1.dv + 2·k2.dv + 2·k3.dv + k4.dv) / 6
Each derivative call performs a full atmosphere lookup, computes Mach and Reynolds number, evaluates all seven drag components (including the power-on plume correction when the motor is burning), reads thrust from the interpolated thrust curve, and returns net acceleration.
a_net = (Thrust(t) − F_drag(h,v,t)) / m(t) − g
Mass decreases linearly through the burn as propellant is depleted. After burnout, mass is constant at the dry mass. Gravity uses g = 9.80665 m/s².
The integrator runs at 50 ms internally. Output is recorded every 100 ms (every second step), producing ~300 data points per channel for a typical 30-second flight.
Air density, temperature, pressure, speed of sound, and dynamic viscosity are computed from altitude using the ISA 1976. The model covers three layers up to 32 km — well above the ceiling of any HPR rocket.
| Layer | Altitude range | Lapse rate | Base temperature |
|---|---|---|---|
| Troposphere | 0 – 11 km | −6.5 K/km | 288.15 K (15 °C) |
| Tropopause | 11 – 20 km | 0 K/km (isothermal) | 216.65 K (−56.5 °C) |
| Lower stratosphere | 20 – 32 km | 0 K/km (isothermal) | 216.65 K |
T(h) = T0 − L·h // T0=288.15 K, L=0.0065 K/m (troposphere) P(h) = P0 · (T(h)/T0)^(g/L·R) // P0=101325 Pa ρ(h) = P(h) / (R·T(h)) // R=287.05 J/kg·K a(h) = √(γ · R · T) // γ=1.4 (speed of sound) μ(h) = 1.458×10⁻⁶ · T^1.5 / (T + 110.4) // Sutherland's law (viscosity)
All five quantities update at every integration step, correctly accounting for the large air density change between launch altitude and apogee. The launch altitude input allows the simulation to start at field elevation — important for high-desert launch sites where sea-level assumptions would overestimate drag.
With correct geometry, surface roughness, and fin profile inputs, APEX's analytical model matches real flight data without any correction factor. Two benchmark flights spanning the subsonic and supersonic regimes:
Surface roughness is the most influential calibration input. Body Roughness and Fin Roughness should match the actual finish of each surface:
| Surface | Roughness (µm) |
|---|---|
| Polished metal / bare carbon fibre | 0.5 – 2 |
| Smooth paint (glossy finish) | 5 – 20 |
| Matte / spray paint | 30 – 60 |
| Rough fibreglass | 100 – 300 |
If simulated apogee is consistently above measured flights after confirming correct geometry, increase surface roughness. Rougher surfaces add drag and lower apogee.
When a Featherweight Blue Raven altimeter CSV is loaded, the Aerodynamics tab back-calculates actual Cd vs Mach from the coast phase — the interval between burnout and apogee when thrust is zero and drag is the only unknown:
F_drag(t) = m · (g + a_measured(t)) // coast phase: thrust = 0
Cd(t) = F_drag(t) / (0.5 · ρ(h) · v² · A_ref)
This produces a measured Cd vs Mach curve overlaid on the model prediction. A persistent offset usually means the surface roughness or fin profile inputs are wrong — not the model itself.
APEX can also import an OpenRocket Cd CSV as a visual reference overlay on the Aerodynamics chart. This does not replace the APEX drag model — the simulation always uses APEX's computed Cd — but it allows direct visual comparison to identify where the two models diverge.