Post-processing¶

Aeolus ships a layered post-processing toolkit covering the outputs that GT-POST and AVL Cruise M present by default: indicated-engine performance, frequency-domain / NVH analysis, and system-level reports (turbo operating-point overlay, energy balance, cumulative emissions, drive-cycle integration, steady-state convergence detection). All of it sits on top of the standard tracer / HDF5 output — every function is a pure NumPy reduction over time series, so it composes cleanly with custom recorders and external data.

This page describes what each module computes, the underlying formulas, and a representative usage pattern.

Module map¶

Module	What it gives you
`aeolus1d.engine_post`	Per-cycle IMEP / PMEP / IMEP_gross / p-max / heat release / CA10/50/90 / polytropic exponents / volumetric efficiency / COV(IMEP)
`aeolus1d.analysis.engine_metrics`	Live `CycleAccumulator` + `MultiCycleAccumulator` (multi-cylinder fleet rollup) used during a run via `on_step=`
`aeolus1d.analysis.acoustics`	SPL / dB(A) / dB(C) / dB(Z) / spectrum / 1·1·octave + 1·3·octave bands (brick-wall or Butterworth roll-off) / engine-order analysis / Campbell diagram / spectrogram / muffler transmission (naive + 4-mic decomposition) and insertion loss / monopole free-field radiation / combustion-noise level
`aeolus1d.analysis.system_reports`	Turbo operating-point overlay + map envelope / Sankey-ready energy balance / cumulative emissions / drive-cycle integrator / steady-state convergence detector
`aeolus1d.audio`	WAV-file export of a probe pressure trace (sub-cell interpolation → polyphase anti-alias decimation → single-pole DC high-pass → 16-bit PCM)

Engine cycle diagnostics¶

The aeolus1d.engine_post module reduces a (t, p, V) trace from a cylinder probe into the engine-textbook performance metrics. The trace can come from any source: a Tracer HDF5, a custom on-step recorder, or even measured pressure data fed in directly.

IMEP, gross / net / PMEP¶

from aeolus1d.engine_post import engine_metrics

em = engine_metrics(
    t, p, V,
    V_disp=5.0e-4,
    omega_crank_rad_s=engine_omega,
    ivc_deg=200.0, evo_deg=540.0,           # Phase 7g.2: gross/PMEP partition
    m_air_trapped_kg=4.5e-4,                # Phase 7g.2: η_vol
    rho_intake_ref_kg_m3=1.0,
    m_fuel_per_cycle=5e-5, LHV=4.35e7,      # ISFC + η_thermal
)

print(f"IMEP_net    : {em.imep_mean_pa/1e5:.2f} bar")
print(f"IMEP_gross  : {em.imep_gross_mean_pa/1e5:.2f} bar")
print(f"PMEP        : {em.pmep_mean_pa/1e5:.3f} bar")
print(f"η_vol       : {em.eta_volumetric:.3f}")
print(f"η_thermal   : {em.eta_thermal:.3f}")
print(f"ISFC        : {em.isfc_g_kwh:.0f} g/kWh")
print(f"COV(IMEP)   : {em.cov_imep*100:.2f} %")

Mean effective pressure is the cycle's net work normalised to displacement volume:

\[ \mathrm{IMEP}_\text{net} = \frac{\oint p\,\mathrm{d}V}{V_\text{disp}} \]

The gross / net split depends on where the cycle is integrated. With IVC (intake-valve close) and EVO (exhaust-valve open) in the 720° four-stroke convention, cycle_metrics(...) partitions each cycle into:

Closed loop (IVC → EVO, compression + power) → IMEP_gross
Gas-exchange loop (the rest of the 720°) → contributes to IMEP_net via the pumping work; PMEP = IMEP_gross − IMEP_net

This is the Heywood §2.4 convention. Without ivc_deg/evo_deg the IMEP_gross and PMEP fields stay at their backward-compatible defaults (gross = net, PMEP = 0).

Heat release / mass-fraction-burned¶

from aeolus1d.engine_post import heat_release_analysis

hr = heat_release_analysis(
    p_cycle, V_cycle, theta_deg_cycle,
    gamma=1.35,
    ivc_deg=200.0, evo_deg=480.0,
)

print(f"CA10 / CA50 / CA90 : {hr.ca10_deg:.1f} / {hr.ca50_deg:.1f} / {hr.ca90_deg:.1f} deg")
print(f"Burn duration       : {hr.burn_duration_deg:.1f} deg")
print(f"Q_total             : {hr.Q_total_j:.1f} J")
print(f"Polytropic n (comp) : {hr.n_compression:.3f}")
print(f"Polytropic n (exp)  : {hr.n_expansion:.3f}")

Inverts the first-law balance on the closed cycle:

\[ \frac{\mathrm{d}Q_n}{\mathrm{d}\theta} = \frac{\gamma}{\gamma - 1}\,p\,\frac{\mathrm{d}V}{\mathrm{d}\theta} + \frac{1}{\gamma - 1}\,V\,\frac{\mathrm{d}p}{\mathrm{d}\theta} \]

The cumulative integral normalised to its asymptote is the mass-fraction-burned curve x_b(θ) ∈ [0, 1]. CA10 / CA50 / CA90 are the crank angles where x_b crosses 0.10 / 0.50 / 0.90 — the canonical combustion-phasing indicators. Polytropic exponents are fitted on the early-compression and late-expansion log p / log V segments; values close to γ indicate negligible heat loss / blowby.

Volumetric efficiency¶

from aeolus1d.engine_post import volumetric_efficiency

eta_vol = volumetric_efficiency(
    m_air_trapped_kg=4.5e-4,
    rho_ref_kg_m3=intake_manifold_rho,
    V_disp_m3=0.5e-3,
)

\[ \eta_\text{vol} = \frac{m_\text{air,trapped}}{\rho_\text{ref}\,V_\text{disp}} \]

1 indicates tuned-intake / forced-induction; < 1 reflects throttling or heating losses.

COV of IMEP¶

from aeolus1d.engine_post import cov_imep, cycle_metrics

cycles = cycle_metrics(t, p, V, V_disp=0.5e-3, omega_crank_rad_s=...)
cov = cov_imep(cycles)

\[ \mathrm{COV}(\mathrm{IMEP}) = \frac{\sigma(\mathrm{IMEP})}{|\bar{\mathrm{IMEP}}|} \]

The combustion-stability indicator. ≈ 5 % is the rule-of-thumb upper bound for a well-running gasoline engine; > 10 % flags misfire / unstable combustion.

Acoustic / NVH analysis¶

The aeolus1d.analysis.acoustics module ships the frequency-domain post-processing that GT-POST and AVL Cruise M ship for muffler / intake-runner tuning workflows. Inputs are plain probe pressure arrays, so the module composes with the Tracer HDF5 output and the WAV exporter in aeolus1d.audio.

SPL — sound pressure level¶

from aeolus1d.analysis.acoustics import spl_db, spl_db_a, spl_db_c, spl_db_z

linear_db   = spl_db(probe_p, p_ref=20e-6)         # dB (un-weighted)
weighted_db = spl_db_a(probe_p, sample_rate_hz)    # dB(A) per IEC 61672
c_db        = spl_db_c(probe_p, sample_rate_hz)    # dB(C) — peaky / LF-rich noise
z_db        = spl_db_z(probe_p, sample_rate_hz)    # dB(Z) — flat (≡ spl_db)

Subtracts the static-pressure DC component before computing RMS (only the AC component radiates). The default p_ref = 20 µPa is the human-hearing reference in air (use 1 µPa for underwater). dB(A) applies the standard IEC 61672 A-weighting curve — the regulatory noise-pollution metric that attenuates the low band and boosts 1–4 kHz. dB(C) (IEC 61672 C-curve) is nearly flat across the audible band and is the weighting for peak / low-frequency-rich noise such as exhaust boom. dB(Z) is the zero (flat) weighting — identical to the un-weighted spl_db, provided for API symmetry (its sample_rate_hz argument is unused).

Power spectrum + octave bands¶

from aeolus1d.analysis.acoustics import power_spectrum, octave_band_levels

freqs, amp = power_spectrum(probe_p, fs, window="hann")
bands_one_third = octave_band_levels(freqs, amp, n_per_octave=3)
for b in bands_one_third:
    print(f"{b.centre_hz:>7.1f} Hz  {b.spl_db:5.1f} dB")

Returns RMS-normalised one-sided amplitude spectrum (Pa per bin) and IEC 61260-preferred 1/1 or 1/3 octave-band reductions. A pure sine of peak amplitude A shows up at A/√2 in the right bin.

By default each band is an ideal brick wall — the quadratic sum of the bins strictly inside the base-2 band edges. Pass rolloff_order=<order> to instead weight every bin by a Butterworth band-pass analyzer filter of that order (−3 dB at the nominal edges, asymptotic skirt slope 6·order dB/octave per side), modelling a real 1/1- or 1/3-octave analyzer's finite selectivity — adjacent bands overlap and out-of-band energy leaks in. Higher orders approach the brick-wall limit.

Order analysis¶

For a varying-RPM trace (engine run-up / run-down), the canonical NVH plot is amplitude vs frequency vs engine order:

from aeolus1d.analysis.acoustics import order_analysis

track = order_analysis(
    probe_p, sample_rate_hz=fs,
    omega_crank_t_s=omega_t_axis,
    omega_crank_rad_s=omega_history,
    orders=[0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 8],
    window_s=0.05, overlap=0.5,
)
# track.amplitude_pa[i, j] is the amplitude of orders[j] at rpm[i]

Engine "orders" are integer / half-integer multiples of crank rotation frequency. Order 2 is the dominant gas-pressure harmonic in a 4-cylinder 4-stroke; order 0.5 is the firing frequency of a single-cylinder 4-stroke; order 1 is once-per-rev (reciprocating- mass imbalance harmonics). The function slides a Hann-windowed FFT across the trace and samples each window's spectrum at every requested order (default orders 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 8).

Campbell diagram¶

from aeolus1d.analysis.acoustics import campbell_diagram

cd = campbell_diagram(track, resonances_hz=[220.0, 550.0])
for x in cd.crossings:
    print(f"order {x.order:g} crosses {x.resonance_hz:.0f} Hz "
          f"at {x.crossing_rpm:.0f} rpm  (amp {x.amplitude_pa:.1f} Pa)")

Builds the interference (Campbell) diagram from an OrderTrack. The order lines are f_order(rpm) = (rpm/60)·order; a structural / acoustic resonance f_r is excited by order o at the critical speed rpm = 60·f_r/o. Each supplied resonance × order whose crossing speed falls inside the swept RPM range yields a CampbellCrossing with the order amplitude interpolated at the crossing — a direct read of how strongly that resonance is driven. resonances_hz may be omitted (returns the order_freqs_hz grid only, for plotting).

Spectrogram¶

from aeolus1d.analysis.acoustics import spectrogram

times, freqs, mag = spectrogram(probe_p, fs, window_size=1024, overlap=0.5)
# mag[i, j] is the RMS-normalised amplitude of freqs[i] at times[j]

Manual STFT (Hann-windowed, coherent-gain normalised, np.fft only — no scipy) for run-up / waterfall plots. Uses the same tonal normalisation as power_spectrum / order_analysis (a pure tone reads at its true amplitude, not 1.76 dB low).

Muffler design — transmission loss / insertion loss¶

from aeolus1d.analysis.acoustics import transmission_loss, insertion_loss

# Two probes straddling a muffler element:
freqs, tl_db = transmission_loss(p_upstream, p_downstream, fs)

# Two runs of the same case, with and without the muffler:
freqs, il_db = insertion_loss(p_baseline, p_with_muffler, fs)

TL is a property of the element alone: TL(f) = 20·log10(|P_in(f)| / |P_out(f)|). Positive means the element attenuates that frequency.
IL is the in-situ noise reduction the element actually delivers — TL plus the impedance-match interaction with the rest of the system. The metric that matters for end-user noise rating.

Naive TL is anechoic-only

transmission_loss is the naive total-pressure amplitude ratio at each probe — correct only for a reflection-free (anechoic-terminated, no upstream standing wave) setup. Wherever a reflection exists the upstream probe sees incident + reflected superimposed and the metric over/under-states TL.

For reflective terminations use the four-microphone wave-decomposition TL (Munjal §4; Chung & Blaser 1980):

from aeolus1d.analysis.acoustics import (
    transmission_loss_decomposed, decompose_waves,
)

# Two probe pairs straddling the element, each spaced `spacing_m`:
freqs, tl_db = transmission_loss_decomposed(
    p_up1, p_up2, p_dn1, p_dn2,
    spacing_m=0.05, sample_rate_hz=fs, c_sound=343.0,
)

decompose_waves(p1, p2, spacing_m=..., sample_rate_hz=..., c_sound=...) is the underlying two-microphone forward/backward plane-wave split (returns complex (freqs, A, B)); transmission_loss_decomposed isolates the incident forward wave upstream and the transmitted forward wave downstream and takes their amplitude ratio. The decomposition is singular at ks = nπ (spacing = n·λ/2), so the valid band is roughly f < c/(2·spacing).

Radiated & combustion noise¶

The in-pipe or in-cylinder pressure is not the listener SPL — two helpers close that gap:

from aeolus1d.analysis.acoustics import (
    freefield_monopole_pressure, combustion_noise_db,
)

# Orifice → far-field mic as an acoustic monopole (compact source):
p_mic = freefield_monopole_pressure(
    u_orifice_trace, area_m2=A, distance_m=1.0,
    sample_rate_hz=fs, rho0=1.2, c_sound=343.0,
)   # feed to spl_db() / power_spectrum() for the listener spectrum

# In-cylinder pressure → radiated combustion-noise level [dB]:
level_db = combustion_noise_db(p_cyl_trace, fs)

freefield_monopole_pressure applies the compact-monopole far field p(r,t) = ρ₀·S/(4π·r)·du/dt|_(t − r/c) (volume velocity Q = S·u), with an integer-sample retarded-time shift for the propagation delay r/c. combustion_noise_db filters the cylinder-pressure spectrum by a frequency-dependent structure attenuation A(f) (cylinder pressure → radiated surface noise; defaults to generic_structure_attenuation(), an illustrative block low-pass — pass your own (freqs, db) table for a real engine) and collapses to an overall RMS level. Sharper pressure rise (knock / premixed Diesel burn) carries more high-frequency content, which the block transmits to the surface.

System-level reports¶

The aeolus1d.analysis.system_reports module ships the system-level plots GT-POST / Cruise M ship as default reports.

Turbo operating-point overlay¶

from aeolus1d.analysis.system_reports import (
    turbo_map_envelope, turbo_operating_points,
)

# Envelope sampled from the map's PR-grid edges along each speed line:
env = turbo_map_envelope(compressor.map)

# Per-step trajectory (your run records dicts or component snapshots):
op = turbo_operating_points(times, comp_snapshots)

# Plot env.{pr,mdot}_{surge,choke}_per_line as the static envelope
# and (op.mdot_corr, op.pr) as the trajectory through it.
# op.regime carries the Phase 7c-surge regime flag per sample.

Energy balance (Sankey)¶

from aeolus1d.analysis.system_reports import energy_balance

eb = energy_balance(
    m_fuel_per_cycle_kg=5e-5,
    LHV_J_per_kg=4.35e7,
    indicated_work_j=800.0,
    pumping_work_j=50.0,
    exhaust_enthalpy_j=900.0,
    wall_heat_j=300.0,
)
print(eb.fractions)   # → {'indicated_work': 0.37, 'pumping_work': 0.02, ...}

Per-cycle first-law partition: fuel chemical energy → indicated work + pumping loss + exhaust sensible enthalpy + wall heat + residual. eb.fractions returns the Sankey-ready percentages summing to 1.0. A large residual flags unmodelled losses or integration errors.

Cumulative emissions¶

from aeolus1d.analysis.system_reports import cumulative_emissions

totals = cumulative_emissions(
    times_s,
    mdot_kg_s=tail_pipe_mdot,
    Y_species={"CO": Y_CO_history, "NOx": Y_NOx_history},
)
print(f"Total CO emitted : {totals['CO']*1000:.3f} g")
print(f"Total NOx emitted: {totals['NOx']*1000:.3f} g")

Trapezoidal integral of ∫ ṁ · Y · dt for each species — used for post-test-cycle emissions reporting (g/test-cycle, equivalent to the regulatory metrics from FTP / WLTP / NEDC test runs).

Drive-cycle integration¶

from aeolus1d.analysis.system_reports import (
    DriveCyclePoint, drive_cycle_runner,
)

# 1. Build the steady-state operating-point map once.
def my_engine_op(rpm, torque):
    # Lookup / regression / re-run aeolus1d at the operating point.
    # Must return at minimum:
    return {
        "fuel_rate_kg_s": ...,
        "power_w": rpm * torque * 2 * np.pi / 60,
        "CO_rate_kg_s": ...,
        "NOx_rate_kg_s": ...,
    }

# 2. Author the drive cycle (NEDC / WLTP / FTP / your own).
schedule = [
    DriveCyclePoint(t_s=0.0,    rpm=800,   torque_nm=0),
    DriveCyclePoint(t_s=10.0,   rpm=1500,  torque_nm=60),
    # ... full cycle ...
]

# 3. Integrate.
res = drive_cycle_runner(schedule, operating_point_callable=my_engine_op)
print(f"Total fuel : {res.total_fuel_kg*1000:.1f} g")
print(f"CO / NOx   : {res.emissions_total_kg}")

The harness midpoint-rules the (rpm, torque) state across each segment, calls the user's operating_point_callable, and accumulates fuel + work + per-species emissions. Build the steady-state engine map once (one aeolus1d run per operating point, or a regression of measured data); drive any cycle through it without re-running the 1-D solver every step.

Steady-state convergence detection¶

from aeolus1d.analysis.system_reports import is_steady_state

# Inside an on_step recorder accumulating per-cycle IMEP:
report = is_steady_state(
    cycle_imep_pa_history,
    cov_threshold=0.02,    # 2 % rule-of-thumb
    window=5,
)
if report.converged:
    print(f"converged after {len(history)} cycles "
          f"(COV = {report.cov_imep_window*100:.2f} %)")

Auto-detect when an engine sim has reached a periodic limit cycle by checking COV(IMEP) over the most-recent N cycles. The 2 % default is the rule-of-thumb convergence criterion for 1-D wave- action engine sims — saves you from over-running at steady state.

Audible-noise export¶

For listening to the engine note rather than measuring it. A probe is an AudioProbeSpec(pipe, x_frac) — build it directly or parse a CLI-style string pipe:<id>:<x_frac> (the pipe: prefix and the x_frac are both optional; x_frac defaults to 0.5, the pipe centre):

from aeolus1d.audio import (
    render_probe_audio, write_wav, waveform_diagnostics,
    parse_probe_string, AudioProbeSpec,
)

probe = parse_probe_string("pipe:runner1:0.95")   # or AudioProbeSpec("runner1", 0.95)
waveform, rate = render_probe_audio(
    "run.h5", probe, rate=44100, duration=None,
)
write_wav("engine.wav", waveform, rate)
print(waveform_diagnostics(waveform, rate))   # duration, peak, RMS dBFS, audible?

render_probe_audio reads the tracer time series from the HDF5 output (layout /time_series/t + /time_series/pipes/<id>/p), linearly interpolates a sub-cell pressure trace at x_frac, and returns (waveform, rate) — a float-32 array in [-1, 1] plus its sample rate. It does not write a file; call write_wav(path, waveform, rate) to persist the normalised 16-bit PCM. The conversion (pressure_to_audio) removes the ~1 atm DC baseline with a single-pole high-pass (default cutoff 20 Hz), resamples onto a uniform audio-rate grid — decimating with a scipy.signal.resample_poly polyphase anti-alias filter when the solver cadence is finer than the audio rate (bare point-sampling would alias super-Nyquist content into the audible band) — and normalises to a 0.95 headroom target. waveform_diagnostics reports duration, peak, RMS in dBFS, and a coarse audible heuristic (RMS below −50 dBFS usually means the trace was too short / too uniform). Pair with spl_db_a() on the source trace to know how loud the run is in psycho-acoustic terms.

Both this and the NVH metrics are also exposed as CLI subcommands over a stored run: aeolus1d audio <run.h5> --probe pipe:<id>:<x_frac> --out file.wav and aeolus1d acoustics <run.h5> --probe pipe:<id>:<x_frac> (SPL / dB(A) / dB(C) + octave bands, plus transmission loss with --probe-out). Both require an HDF5 written with aeolus1d run --trace-dt <DT> so the time-series group is populated.