Skip to content

Acoustic mode analysis

Aeolus ships two complementary tools for analysing the natural acoustic modes of a 1-D pipe / chamber geometry. Both are deliberately separate from the time-domain Euler solver: they answer the question "what natural frequencies and mode shapes does this geometry support?" rather than "how does this gas flow over time?".

Tool When to use it Cost Fidelity
aeolus1d.modes — high-order FEM eigensolver Variable-area duct, multiple modes, mode shapes needed A few ms Most physically accurate (full 1-D wave physics)
HelmholtzResonator — 0-D textbook primitive Lumped Helmholtz analysis, runner sizing, design-stage quick-checks Microseconds Textbook 2-DOF (volume + neck inertia); accurate while L_neck ≪ λ/4

The time-domain solver still runs separately. The eigensolvers are for spectral analysis; the time-domain solver is for transient analysis. They cross-validate via mode-projection initial conditions.

High-order FEM eigensolver

aeolus1d.modes.solve_modes runs a continuous-Galerkin solve on the variable-area duct eigenproblem, in the acoustic-pressure formulation

\[ \frac{\partial}{\partial x} \left( \frac{A(x)}{\rho(x)}\, \frac{\partial p}{\partial x} \right) = -\omega^2\, \frac{A(x)}{\rho(x)\, c(x)^2}\, p \]

with Lagrange basis polynomials of arbitrary order p. Returns a ModesResult carrying natural frequencies (Hz, ascending) + mode shapes p_n(x) sampled at the global DOF positions.

The unknown is the acoustic pressure, so a shared continuous- Galerkin DOF enforces pressure continuity and the natural (flux) interface condition conserves volume velocity, \(\sum A/\rho \cdot \partial p / \partial x = 0\) — the physically correct junction conditions when pipes of different temperature (hence ρ, c) meet. (An earlier velocity-potential weighting enforced potential continuity instead, which is wrong across a temperature jump; for a uniform single pipe the two formulations differ only by a constant factor on both matrices, so their eigenfrequencies are identical.) The local mass and stiffness matrices are

\[ M^e_{ij} = \int_e \frac{A}{\rho\, c^2}\, L_i\, L_j \; \mathrm{d}x, \qquad K^e_{ij} = \int_e \frac{A}{\rho}\, L_i'\, L_j' \; \mathrm{d}x, \]

integrated with Gauss-Legendre quadrature (2p points) and assembled into the generalised symmetric eigenproblem \(K\, p = \omega^2 M\, p\), solved by scipy.linalg.eigh.

Quick start

from aeolus1d.modes import solve_modes, BC_NEUMANN, BC_DIRICHLET

result = solve_modes(
    L=1.0,                        # pipe length (m)
    rho=1.225,                    # gas density (kg/m³)
    A=1.0e-3,                     # cross-section (m²)
    c=343.0,                      # sound speed (m/s)
    n_elements=8, p=4,            # mesh + polynomial order
    bc_left=BC_NEUMANN,           # closed / rigid wall (∂p/∂x = 0)
    bc_right=BC_DIRICHLET,        # open (p = 0)
    n_modes=10,
)

print(result.frequencies[:3])     # first 3 natural frequencies (Hz)
print(result.x.shape)             # (n_dof,) physical coordinates
print(result.mode_shapes.shape)   # (n_dof, n_modes)

When to bump p vs n_elements

solve_modes defaults to p=4 and n_elements=8 — that's already 1e-3 % accurate on the first 5 modes of a uniform pipe. For higher modes or sharper accuracy, bump p before n_elements: the FEM converges exponentially in p (spectral convergence) and only algebraically in n_elements (h-refinement). On a uniform pipe, bumping p from 2 to 4 drops the first-mode eigenvalue error from ~10⁻³ to ~10⁻¹⁰ — many orders per p-step — whereas doubling n_elements at fixed p buys only the fixed algebraic h-rate. So reach for p first; grow n_elements mainly to resolve sharp area transitions the elements can't span smoothly.

Variable-area geometry

Pass a callable for A:

import numpy as np

def A_helmholtz(x):
    """A neck (x < 0.1) of A_n = 1e-3, then a chamber of A_c = 2e-2."""
    return np.where(x < 0.1, 1.0e-3, 2.0e-2)

result = solve_modes(L=0.15, rho=1.225, A=A_helmholtz, c=343.0, ...)

For sharper area transitions, build the mesh by hand with distinct elements per region — aeolus1d.modes.Element exposes the per-element constructor:

from aeolus1d.modes import Element

mesh = []
for i in range(4):              # 4 neck elements
    mesh.append(Element(
        x_left=i * 0.025, x_right=(i + 1) * 0.025,
        p=6, rho_func=rho_const, A_func=A_neck, c_func=c_const,
    ))
for i in range(4):              # 4 chamber elements
    mesh.append(Element(
        x_left=0.1 + i * 0.0125, x_right=0.1 + (i + 1) * 0.0125,
        p=6, rho_func=rho_const, A_func=A_chamber, c_func=c_const,
    ))
result = solve_modes(L=0.15, rho=1.225, A=1.0, c=343.0, mesh=mesh, ...)

Multi-pipe networks

aeolus1d.modes.solve_modes_network(case, ...) extends the solve from one pipe to a whole Case. Every pipe contributes its element matrices, and pipe ends that meet at a junction share a single network-global DOF (found by union-find over the pipe-end DOFs) — enforcing acoustic-pressure continuity, with volume-velocity conservation as the natural flux condition. A junction that carries a lumped chamber volume (plenum / ehr kinds) additionally adds its acoustic compliance \(V/(\rho c^2)\) to the shared DOF's mass entry, giving Helmholtz-style storage; chamber-free kinds (corberan / passthrough / biased_plenum) stay ideal rigid connection points.

Free pipe ends — those not in any junction — take their FEM BC from the case's [[bc]] list:

Case BC kind FEM condition
wall, valve_inlet, valve_orifice Neumann (∂p/∂x = 0) — velocity held
stagnation_inlet, pressure_outlet, transient_stagnation_inlet, transient_pressure_outlet Dirichlet (p = 0) — pressure held
transmissive, nonreflecting Dirichlet (p = 0)

Valve ends idealise as closed at θ = 0 (conservative); both anechoic-outflow kinds idealise as an open end — a truly non-reflecting termination has no lossless modal equivalent. Pass bc_overrides={(pipe_id, end): kind} to force a specific BC at any end for a parameter sweep. Each pipe's reference (ρ, c) comes from its init.p / init.T (default 1 atm / 300 K); eigenfrequencies scale with c/L and are insensitive to the absolute reference value. (The reference state is read the same way as by the CLI and GUI, through aeolus1d.modes.reference.reference_state: a (p, rho) or mode-init pipe is solved at its own temperature. Before 2026-10-07 solve_modes_network read only init.p / init.T and assumed 300 K otherwise. A BC kind with no acoustic mapping, a BC on a pipe end that does not exist and a junction port on an unknown pipe now warn; the end is treated as rigid.)

The call returns a NetworkModesResult with one ascending frequencies array plus per_pipe_shapes / per_pipe_x mappings (mode shape sampled on each pipe's own DOFs) and junction_shapes (the shared-DOF value per junction), so a plot script can stitch the network back together without tracking the global indexing. This is a Python-API entry point — the CLI and GUI docks below analyse a single pipe.

from aeolus1d.io.case import load_case
from aeolus1d.modes import solve_modes_network

case = load_case("examples/38_engine_turbo_scr.toml")
res = solve_modes_network(case, p=4, n_elements=8, n_modes=8)
print(res.frequencies)              # (n_modes,) Hz, ascending
print(list(res.per_pipe_shapes))    # pipe ids with resolved mode shapes

CLI access

aeolus1d modes <case.toml> runs the eigensolver on one pipe of a TOML case and prints the lowest N frequencies; --out FILE.csv writes one row per FEM DOF for downstream plotting:

aeolus1d modes examples/01_single_pipe_sod.toml --bc-left wall --bc-right wall --n-modes 5
aeolus1d modes examples/38_engine_turbo_scr.toml --pipe intake_atm --p 6 --n-modes 8 --out /tmp/intake_modes.csv

GUI access

Launch aeolus1d gui, select a pipe, switch to the Modes dock (tabbed alongside the network view + 2D plot). The dock shows a frequency table for the selected pipe; click any row to see its mode shape on the matplotlib panel. BC selectors (wall / open), p and n_elements spinboxes, plus a "Recompute" button.

Lumped Helmholtz resonator

aeolus1d.components.helmholtz_resonator.HelmholtzResonator is the textbook 0-D resonator: a Volume connected to a far-field reservoir through a lumped inertial neck. State carried as (m_V, E_V, u_n); equations of motion are the Helmholtz 2-DOF oscillator

\[ m_n \frac{du_n}{dt} = A_n \, (p_{\mathrm{far}} - p_V), \qquad \frac{dm_V}{dt} = \rho_{\mathrm{upwind}} \, u_n \, A_n. \]

There is no 1-D pipe involved on the neck side, so the lumped mode emerges with no wave physics in the neck at all. (This page used to say a 1-D pipe + 0-D Volume coupling cannot ring at f_H. Retracted 2026-09-27: that bench had closed the neck's far end with a wall, which is not a Helmholtz resonator; opened to the surroundings, the 1-D neck + 0-D Volume measures f_H to 0.26 %, and tests/integration/test_helmholtz.py now passes.) The 0-D HelmholtzResonator's gate tests/integration/test_helmholtz_resonator.py sweeps V / A_n / L_neck and holds the simulated peak within 1 % of f_H.

Quick start

from aeolus1d.components.volume import Volume
from aeolus1d.components.helmholtz_resonator import (
    HelmholtzResonator, integrate, rayleigh_end_correction,
)

# Chamber: 1 L volume, slightly over-pressurised.
vol = Volume(V=1.0e-3, gamma=1.4, R_gas=287.05)
vol.set_state(p=1.01e5, T=300.0)

# Neck: 50 mm geometric length, 1 cm² cross-section.
A_n = 1.0e-4
L_eff = 0.05 + rayleigh_end_correction(A_n)

res = HelmholtzResonator(
    volume=vol, A_neck=A_n, L_eff=L_eff,
    p_far=1.0e5, T_far=300.0,
)

print(f"f_H = {res.f_helmholtz():.2f} Hz")

# Integrate 30 periods at 200 substeps per period (Heun RK2).
trace = integrate(res, t_end=30.0 / res.f_helmholtz(),
                  dt=1.0 / (res.f_helmholtz() * 200))

Validity envelope

The lumped formula is exact only when the neck length is much shorter than a quarter wavelength (L_eff ≪ λ_H / 4). At L_eff / (λ/4) = 0.4 the lumped formula is still good to ~5 % vs the FEM solve; beyond that, prefer the FEM. The HelmholtzResonator validation tests cover the envelope L_eff / (λ/4) ∈ [0.07, 0.41], where the simulated peak agrees with the textbook formula to 0.025 % — that's the cheapest side of the FEM/lumped tradeoff when the geometry actually is a textbook resonator.

Cross-validation: mode-projection IC

The eigensolvers and the time-domain Euler solver should agree on the spectrum of any geometry — that's the consistency check behind the mode-projection initial conditions (tests/integration/test_mode_init.py). Recipe:

  1. Pick a pipe, decide the BCs, run solve_modes to get the target eigenfrequency f_n and mode shape φ_n(x).
  2. Initialise a transient run with [pipe.init] kind = "mode", mode_index = n — the loader runs solve_modes on the pipe geometry, skips any rigid-body (f < 1 Hz) mode so mode_index = 1 is the first physical mode, normalises the shape to peak \(|\varphi_n| = 1\), resamples it at the pipe's cell centres, and lays down a quiescent (u = 0) field with an equal fractional over-pressure / over-density in that shape: \(p = p_{\mathrm{ref}}\,(1 + \varepsilon\,\varphi_n)\) and \(\rho = \rho_{\mathrm{ref}}\,(1 + \varepsilon\,\varphi_n)\). The amplitude \(\varepsilon\) is amplitude; p_ref / T_ref set the reference state and bc_left / bc_right (plus n_elements / p_order) choose the FEM problem.
  3. Run for ~20 periods, FFT the mid-pipe pressure trace.
  4. The peak frequency should land within ~5 % of f_n — that's the time-domain solver's truncation error, well below the FEM's accuracy.

A clean diagnostic for solver-versus-spectrum drift, particularly useful when adding new closure terms (friction, heat transfer) that should shift the spectrum predictably.

[case]
name = "mode_one_validation"
t_end = 0.116    # 20 / f_1
cfl = 0.4
scheme = "muscl"
limiter = "minmod"

[gas]
gamma = 1.4
R_gas = 287.05

[[pipe]]
id = "tube"
length = 1.0
n_cells = 200
area = 1.0e-3

  [pipe.init]
  kind = "mode"
  mode_index = 1
  amplitude = 1.0e-3
  p_ref = 1.0e5
  T_ref = 300.0
  bc_left = "wall"
  bc_right = "open"

[[bc]]
pipe = "tube"
end = "left"
kind = "wall"

[[bc]]
pipe = "tube"
end = "right"
kind = "pressure_outlet"
p_back = 1.0e5