Port‑Hamiltonian Fluid‑Acoustic Engine Sound Synthesis
Implementation Guide – Version 2.0

Table of Contents

1. Overview & Design Philosophy
2. Port‑Hamiltonian Framework
   1. State‑space definition
   2. The input matrix $\BB(\xx)$
   3. Power balance & passivity
   4. The direct feed‑through matrix $\DD(\xx)$
3. 1D Acoustic Pipe – PHSPipe1D
   1. Continuous PHS model
   2. Staggered‑grid discretisation
   3. Boundary port variables
   4. Viscothermal losses (Zwikker–Kosten theory)
   5. CFL condition
4. 0D Volumes
   1. Acoustic compliance – PHSVolumeAcoustic
   2. Thermodynamic volume – PHSVolumeThermo
   3. Why DC coupling works automatically in PHS
5. Cylinder – PHSCylinder
   1. PHS formulation with four ports
   2. Slider‑crank kinematics
   3. Combustion heat addition
   4. Valve orifice flow (choked / subsonic)
6. Engine Network Topology
7. Global Solver – Implicit Midpoint Rule
   1. Global state vector assembly
   2. Step‑by‑step algorithm
   3. Sparse matrix structure
   4. Thomas algorithm for tridiagonal blocks
   5. CFL‑free stability
8. Audio Output – Monopole Radiation
9. Stability Guarantees
10. Implementation Roadmap
11. Component Architecture Diagram
12. References
13. Appendices
    1. Nomenclature
    2. Physical Constants
    3. Typical Engine Parameters
    4. Relationship to Digital Waveguide Networks
    5. Recommended Source File Structure

1. Overview & Design Philosophy

This document describes a physically‑based, real‑time audio synthesis engine for combustion‑engine sounds. The simulation is built on a Port‑Hamiltonian System (PHS) framework, guaranteeing passivity, energy conservation, and numerical stability for any interconnection of components.

The core architecture consists of:

• 1D lossy acoustic waveguides (pipes) discretised on a staggered grid.
• 0D thermodynamic volumes (manifolds, cylinders) that exchange mass, energy, and acoustic power with the waveguides.
• Resistive elements (valves, orifices) modelled by compressible‑flow equations.
• Open‑end radiation via a passive one‑pole filter and monopole formula.
• A global implicit‑midpoint time integrator that solves all state equations simultaneously, preserving the Dirac structure and passivity.

The system runs synchronously with audio hardware at 44.1 kHz and produces the radiated sound pressure directly from physical quantities – no arbitrary gains or sample‑based playback.

Design principles (to be repeated in every source file):

1. Hamiltonian – total energy is a quadratic form in the state variables $H = \frac12 \sum C_i p_i^2 + \frac12 \sum L_j U_j^2$ (pipes) and $H_V = \frac{V}{\gamma-1}p_V$ (volumes).
2. Dirac structure – evolution uses $\dot{\xx} = (\JJm-\RRm)\nabla\HH(\xx) + \BB\;\uu_{\text{ext}}$, with $\JJm$ skew‑symmetric, $\RRm \succeq 0$.
3. Midpoint rule – structure‑preserving time integration requires simultaneous evaluation of all state derivatives before updating.
4. Power ports – boundaries enforce effort continuity and flow continuity; flows are computed from the current state and used consistently for all connected components.
5. CFL – for any explicit sub‑step, $c\Delta t/\Delta x \le 1$ must hold.
6. Dissipation – a positive semidefinite $\RRm$ guarantees stability; physical loss (viscothermal) is realised as small velocity damping.
7. Wave scattering – boundaries can use wave variables $(p^+, p^-)$ with reflectance $|R|\le 1$; passivity is automatic.

2. Port‑Hamiltonian Framework

2.1 State‑space definition

A Port‑Hamiltonian System is defined by the equations: $$ \begin{aligned} \dot{\xx} &= \bigl(\JJm(\xx) - \RRm(\xx)\bigr)\,\nabla \HH(\xx) \;+\; \BB(\xx)\,\uu_{\text{ext}},\\[4pt] \yy &= \BB(\xx)^{\!\top}\,\nabla \HH(\xx) \;+\; \DD(\xx)\,\uu_{\text{ext}}, \end{aligned} $$ where the variables have the following meanings:

2.2 The input matrix $\BB(\xx)$

$\BB(\xx)$ (often constant) selects where external inputs enter the system. It is the Jacobian of the port variables with respect to the inputs, or equivalently the matrix that satisfies $\BB(\xx)\,\uu_{\text{ext}} = \sum_k G_k\,\uu_{\text{ext},k}$, where $G_k$ are the input vector fields.

Examples:

• Pipe boundaries: The external inputs are the volume flows $Q_{\text{left}}, Q_{\text{right}}$ entering the pipe at its ends. $\BB$ places these flows into the pressure‑update equations at nodes $0$ and $N$: $$ \BB = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ 0 & 1 \end{bmatrix} \begin{matrix} \leftarrow p_0 \\ \leftarrow U_{1/2} \\ \vdots \\ \leftarrow U_{N-1/2} \\ \leftarrow p_N \end{matrix} $$
• Thermodynamic volume: The inputs are mass flow $\dot{m}_{\text{in}}$ and heat addition $\dot{Q}_{\text{comb}}$. $\BB = I_{2\times2}$ maps them directly to $\dot{m}$ and $\dot{U}$.
• Cylinder piston work: The piston exerts a force on the gas; this appears as a $p\,\dif V/\dif t$ term in $\dot{U}$, which is part of the internal Dirac structure, not an external input.

2.3 Power balance & passivity

Taking the time derivative of the Hamiltonian along trajectories: $$ \frac{\dif\HH}{\dif t} = (\nabla\HH)^{\!\top}\,\dot{\xx} = (\nabla\HH)^{\!\top}\bigl[(\JJm-\RRm)\nabla\HH + \BB\,\uu_{\text{ext}}\bigr] = -(\nabla\HH)^{\!\top}\RRm\,\nabla\HH \;+\; \yy^{\!\top}\uu_{\text{ext}} \;-\; \uu_{\text{ext}}^{\!\top}\DD\,\uu_{\text{ext}}. $$ Since $\RRm \succeq 0$, the first term is $\le 0$ (dissipation). If $\DD + \DD^{\!\top} \succeq 0$, the third term is also non‑positive. Therefore, $$ \frac{\dif\HH}{\dif t} \le \yy^{\!\top}\uu_{\text{ext}}, $$ which is the passivity inequality: the stored energy can increase only due to external power supplied through the ports. When $\uu_{\text{ext}} = 0$, the energy is non‑increasing – the system is Lyapunov stable.

2.4 The direct feed‑through matrix $\DD(\xx)$

In physical models derived from balance equations, $\DD$ is usually zero because the output $\yy$ depends only on the state $\xx$, not directly on the input $\uu_{\text{ext}}$. For our acoustics system $\DD = 0$ everywhere. A non‑zero $\DD$ can arise when algebraic constraints are eliminated or when reduced‑order models are used.

3. 1D Acoustic Pipe – PHSPipe1D

3.1 Continuous PHS model

The linear wave equation for an acoustic duct of constant cross‑section $A$ is written in PHS form using the state variables $\alpha = \frac{p}{\rho_0 c_0^2}$ (dilatation) and $\mu = \rho_0 u$ (momentum density): $$ \begin{aligned} \frac{\pdif}{\pdif t}\begin{bmatrix} \alpha \\ \mu \end{bmatrix} &= \begin{bmatrix} 0 & -\frac{\pdif}{\pdif \xi} \\[2pt] -\frac{\pdif}{\pdif \xi} & 0 \end{bmatrix} \begin{bmatrix} \frac{\delta\HH}{\delta\alpha} \\[2pt] \frac{\delta\HH}{\delta\mu} \end{bmatrix}, \end{aligned} $$ with Hamiltonian $$ \HH = \frac12 \int_0^L \Bigl( \rho_0 c_0^2\,\alpha^2 + \frac{1}{\rho_0}\,\mu^2 \Bigr) A\,\dif\xi. $$ In terms of the more familiar pressure $p$ and volume velocity $U = A u$: $$ \frac{\pdif p}{\pdif t} = -\frac{\rho_0 c_0^2}{A}\,\frac{\pdif U}{\pdif \xi},\qquad \frac{\pdif U}{\pdif t} = -\frac{A}{\rho_0}\,\frac{\pdif p}{\pdif \xi}. $$ The boundary ports are $(p(0), U(0))$ at the left end and $(p(L), -U(L))$ at the right end (sign convention: positive flow enters the domain).

3.2 Staggered‑grid discretisation

Divide the pipe into $N$ segments of length $\Delta x = L/N$. Place pressures $p_i$ at integer nodes $i = 0,\dots,N$ and volume velocities $U_{j+1/2}$ at half‑integer midpoints $j = 0,\dots,N-1$. This is a staggered grid, which is essential for preserving the skew‑symmetry of the gradient operator.

The discrete gradient and divergence operators are: $$ (\DDm p)_{j+1/2} = \frac{p_{j+1} - p_j}{\Delta x},\qquad (\DDm\!\cdot\! U)_i = \frac{U_{i+1/2} - U_{i-1/2}}{\Delta x}\quad (i=1,\dots,N-1). $$ These satisfy $\DDm^{\!\top} = -\DDm\!\cdot$, which is the discrete analogue of integration by parts.

The semi‑discrete equations become: $$ \begin{aligned} \frac{\dif U_{j+1/2}}{\dif t} &= -\frac{A}{\rho_0}\,(\DDm p)_{j+1/2},\\[4pt] \frac{\dif p_i}{\dif t} &= -\frac{\rho_0 c_0^2}{A}\,(\DDm\!\cdot\! U)_i \qquad (i=1,\dots,N-1),\\[4pt] \frac{\dif p_0}{\dif t} &= -\frac{\rho_0 c_0^2}{A}\,\bigl(Q_{\text{left}} - U_{1/2}\bigr),\\[4pt] \frac{\dif p_N}{\dif t} &= -\frac{\rho_0 c_0^2}{A}\,\bigl(U_{N-1/2} + Q_{\text{right}}\bigr), \end{aligned} $$ where $Q_{\text{left}}$ and $Q_{\text{right}}$ are the external volume flows entering the pipe at the boundaries (positive = into the pipe). The characteristic impedance is $Z_0 = \dfrac{\rho_0 c_0}{A}$.

3.3 Boundary port variables

At the left end ($\xi=0$): effort $= p_0$, flow $= -Q_{\text{left}}$ (flow leaving the domain). At the right end ($\xi=L$): effort $= p_N$, flow $= Q_{\text{right}}$ (flow leaving the domain). The power flowing out of the pipe through the ports is $p_0(-Q_{\text{left}}) + p_N Q_{\text{right}}$.

3.4 Viscothermal losses (Zwikker–Kosten theory)

Real pipes exhibit frequency‑dependent damping due to viscous drag and thermal conduction at the walls. The classical theory of Zwikker and Kosten (1949) gives the complex wave number: $$ k(\omega) = \frac{\omega}{c_0}\sqrt{\frac{1 + (\gamma-1)\Psi_h}{1 - \Psi_v}}, $$ where $\Psi_v$ and $\Psi_h$ are functions of the shear and thermal boundary‑layer thicknesses. For audio frequencies in pipes of radius $a \gtrsim 1$ cm, the damping is weak and can be approximated by a Kelvin–Voigt model: the wall shear stress acts as a velocity‑dependent body force.

In our discrete model this is implemented by adding a diagonal dissipation matrix $\RRm$ with entries $\varepsilon \approx 0.001\text{–}0.005$ on the velocity nodes. Each velocity is multiplied by $(1-\varepsilon)$ per time step, giving a frequency‑independent loss rate that approximates the low‑frequency behaviour of the full Zwikker–Kosten model. A more accurate filter (e.g. a Padé approximant of $\sqrt{i\omega}$) can be substituted later.

3.5 CFL condition

For the explicit leapfrog used inside the global midpoint solver, the time step $\Delta t$ must satisfy $$ \frac{c_0\,\Delta t}{\Delta x} \le 1. $$ With typical choices ($L = 1$ m, $N = 20$, $\Delta t \approx 2.27\times10^{-5}$ s) the Courant number is $\approx 0.155$, well within the stable region.

4. 0D Volumes

4.1 Acoustic compliance – PHSVolumeAcoustic

A small cavity behaves like an acoustic capacitor. Linearising the ideal gas law around ambient pressure $p_0$: $$ \frac{\dif p_V}{\dif t} = \frac{\gamma p_0}{V_0}\,Q_{\text{net}}, $$ where $p_V$ is the gauge pressure, $V_0$ the cavity volume, and $Q_{\text{net}}$ the net volume flow entering.

PHS formulation: State $x = p_V$. Hamiltonian $\HH = \dfrac{V_0}{2\gamma p_0}\,p_V^{\,2}$. The port is $(p_V, Q_{\text{net}})$ – effort is pressure, flow is volume velocity. The Dirac structure is trivial ($\JJm = 0$) and the port equation is $\dot{p}_V = \frac{\gamma p_0}{V_0}\,Q_{\text{net}}$.

4.2 Thermodynamic volume – PHSVolumeThermo

For larger volumes where mass and energy balance matters (manifolds, cylinder):

• States: mass $m$ [kg], internal energy $U$ [J].
• Hamiltonian: $\HH = U$.
• Ideal gas law: $p = (\gamma-1)U / V$, $T = U / (m\,c_v)$, $c_v = R/(\gamma-1)$, $c_p = \gamma c_v$.
• Ports:
  • Acoustic: $(p, Q)$ – pressure and volume flow.
  • Thermal/enthalpy: $(h_{\text{in}}, \dot{m}_{\text{in}})$ – specific enthalpy and mass flow.
• Evolution: $$ \begin{aligned} \frac{\dif U}{\dif t} &= -p\,\frac{\dif V}{\dif t} + \sum h_{\text{in}}\,\dot{m}_{\text{in}} + \dot{Q}_{\text{comb}},\\[4pt] \frac{\dif m}{\dif t} &= \sum \dot{m}. \end{aligned} $$

The volume $V$ may be time‑varying (piston) or constant (manifold). When $V$ is constant, $\dif V/\dif t = 0$ and the $p\,\dif V/\dif t$ term vanishes.

4.3 Why DC coupling works automatically in the PHS framework

5. Cylinder – PHSCylinder

5.1 PHS formulation with four ports

The cylinder is the most complex component. It is a thermodynamic volume with a moving wall (piston) and external heat addition (combustion). It exposes four power ports:

5.2 Slider‑crank kinematics

The cylinder volume as a function of crank angle $\theta$ (radians) is: $$ V(\theta) = V_c + \frac{\pi}{4}B^2\,s(\theta), $$ where $V_c$ is the clearance volume, $B$ the bore, and $s(\theta)$ the piston displacement: $$ s(\theta) = \frac{S}{2}\Bigl(1 - \cos\theta\Bigr) + L\Bigl(1 - \sqrt{1 - \bigl(\frac{S}{2L}\sin\theta\bigr)^2}\Bigr), $$ with $S$ the stroke and $L$ the connecting‑rod length. The clearance volume is $V_c = V_d / (r_c - 1)$, where $V_d = (\pi/4)B^2 S$ is the displacement and $r_c$ the compression ratio.

5.3 Combustion heat addition

Combustion is modelled as a constant‑volume heat addition near top‑dead‑centre (TDC) of the compression stroke: $$ \dot{Q}_{\text{comb}}(t) = \frac{E_{\text{fuel}}}{\tau_{\text{burn}}}\,\, \text{rect}\Bigl(\frac{\theta - \theta_{\text{spark}}}{\Delta\theta_{\text{burn}}}\Bigr), $$ where $E_{\text{fuel}}$ is the fuel energy per cycle, $\theta_{\text{spark}} \approx \pi$ (TDC), and $\Delta\theta_{\text{burn}} \approx \pi/6$ (30° burn duration). The heat is added to the internal energy $U$ via the thermal port.

5.4 Valve orifice flow

The intake and exhaust valves are modelled as variable‑area orifices. The mass flow rate is given by the standard compressible‑flow equation: $$ \dot{m} = C_d\,A_{\text{valve}}\,\frac{p_{\text{up}}}{\sqrt{T_{\text{up}}}}\, \begin{cases} \sqrt{\dfrac{\gamma}{R}}\, \Bigl(\dfrac{2}{\gamma+1}\Bigr)^{\!\frac{\gamma+1}{2(\gamma-1)}}, & \text{if } \dfrac{p_{\text{dn}}}{p_{\text{up}}} \le p_{\text{crit}},\\[10pt] \sqrt{\dfrac{2\gamma}{R(\gamma-1)}}\, \Bigl[\Bigl(\dfrac{p_{\text{dn}}}{p_{\text{up}}}\Bigr)^{\!2/\gamma} - \Bigl(\dfrac{p_{\text{dn}}}{p_{\text{up}}}\Bigr)^{\!(\gamma+1)/\gamma}\Bigr]^{1/2}, & \text{otherwise}, \end{cases} $$ where $p_{\text{crit}} = \bigl(\frac{2}{\gamma+1}\bigr)^{\gamma/(\gamma-1)} \approx 0.528$ for air. The curtain area is $A_{\text{valve}} = \pi D_{\text{valve}}\,L_{\text{lift}}$, with $L_{\text{lift}}(\theta)$ given by a sinusoidal cam profile.

6. Engine Network Topology

(Adapted from Baldan et al. 2015, §3.3)

A multi‑cylinder engine is assembled from the following components:

The shared manifolds act as 0‑junctions (common pressure, sum of flows = 0): all cylinder ports and the connecting pipe see the same manifold pressure, and the net volume flow into the manifold determines its pressure change via the compliance equation.

7. Global Solver – Implicit Midpoint Rule

7.1 Global state vector assembly

All component states are concatenated into a single vector: $$ \xx = \bigl[\,\underbrace{p_0^{\text{pipe1}}, U_{1/2}^{\text{pipe1}}, p_1^{\text{pipe1}}, \dots, p_N^{\text{pipe1}}}_{\text{pipe 1}}, \underbrace{p_0^{\text{pipe2}}, \dots}_{\text{pipe 2}}, \underbrace{p_V^{\text{vol1}}, m^{\text{cyl}}, U^{\text{cyl}}, \dots}_{\text{volumes \& cylinder}}\,\bigr]^{\!\top}. $$ The total number of states for a typical engine with 2 pipes of $N=20$ cells each, 2 manifold volumes, and 1 cylinder is $2(2N+1) + 2 + 2 \approx 86$ – trivially small for a modern CPU.

7.2 Step‑by‑step algorithm

1. Assemble the global $\JJm$, $\RRm$ matrices and $\nabla\HH$.
   For linear acoustics $\nabla\HH = Q\,\xx$ where $Q = \operatorname{diag}(C_1^{-1}, L_1^{-1}, \dots)$. The matrices are sparse and structured.
2. Evaluate the external inputs $\uu_{\text{ext}}^{n+1/2}$.
   These come from the valve orifice equations (semi‑implicit, using the half‑step pressure) and combustion heat release.
3. Form the linear system: $$ \bigl[I - \tfrac{\Delta t}{2}(\JJm - \RRm)Q\bigr]\,\xx^{n+1} = \bigl[I + \tfrac{\Delta t}{2}(\JJm - \RRm)Q\bigr]\,\xx^{n} + \Delta t\,\BB\,\uu_{\text{ext}}^{n+1/2}. $$
4. Solve $A\,\xx^{n+1} = b$.
   $A$ is sparse and banded. For fixed valve states (between valve events) $A$ does not change and its LU factorisation can be reused.
5. Extract audio output from the new state (see §8).

7.3 Sparse matrix structure

The matrix $A$ has the following block structure:

• Pipe interior: Tridiagonal blocks coupling $p_i$ and $U_{i\pm1/2}$.
• Pipe–volume boundary: Off‑diagonal entries connecting $p_0$ (or $p_N$) of a pipe to the volume pressure $p_V$.
• Volume interior: $1\times1$ or $2\times2$ diagonal blocks.
• Cylinder–manifold: Coupling through the valve resistive elements, which appear as additional off‑diagonal terms.

The bandwidth is dominated by the pipe tridiagonal blocks, so the total bandwidth is $\approx 3$. Storage and solution cost is $O(n)$, i.e. linear in the number of states.

7.4 Thomas algorithm for tridiagonal blocks

For the pure pipe system (no volumes), $A$ is tridiagonal and can be solved by the Thomas algorithm (tridiagonal matrix algorithm, TDMA) in $O(n)$ time with $O(1)$ extra storage. When volumes are added, the system remains narrowly banded and can be solved by a banded LU factorisation without pivoting (since $A$ is diagonally dominant for small $\Delta t$).

7.5 CFL‑free stability

The implicit midpoint rule is A‑stable and symplectic for linear Hamiltonian systems. There is no CFL restriction on $\Delta t$: the scheme is unconditionally stable. In practice we choose $\Delta t = 1/44100$ s to resolve the audio bandwidth, which is far below any stability limit.

8. Audio Output – Monopole Radiation

At each open pipe end, the volume velocity $U_{\text{open}}$ is read from the state vector. The radiated pressure at a listener distance $r$ is given by the monopole formula: $$ p_{\text{rad}}(t) = \frac{\rho_0}{4\pi r}\,\frac{\dif U_{\text{open}}}{\dif t}, $$ which we approximate by the finite difference $$ p_{\text{rad}}[n] = \frac{\rho_0}{4\pi r}\,\frac{U_{\text{open}}[n] - U_{\text{open}}[n-1]}{\Delta t}. $$ This is the raw physical pressure in Pascals. It is converted to a digital audio sample by multiplying by a single global OutputScale: $$ \text{sample}[n] = \text{OutputScale} \cdot p_{\text{rad}}[n]. $$ With typical engine flows, OutputScale ≈ 5000 gives a comfortable listening level. A soft clipper ($\tanh$) can be applied as a safety limiter without affecting the physics.

9. Stability Guarantees

10. Implementation Roadmap

11. Component Architecture Diagram

12. References

1. V. Trenchant, H. Ramirez, Y. Le Gorrec, and P. Kotyczka, “Finite differences on staggered grids preserving the port‑Hamiltonian structure with application to an acoustic duct,” J. Comput. Phys., vol. 373, pp. 673–697, 2018. DOI
2. J. Alastuey, Y. Le Gorrec, and Y. Wu, “Extended Port‑Hamiltonian 1D framework for physical systems described by hyperbolic PDEs,” submitted to Systems & Control Letters, 2025. (Preprint available on HAL)
3. A. Falaize and D. Roze, “Generic passive‑guaranteed nonlinear interaction model and structure‑ preserving spatial discretization procedure with applications in musical acoustics,” Acta Acustica, 2024. Link
4. S. Baldan, H. Lachambre, S. Delle Monache, and P. Boussard, “Physically informed car engine sound synthesis for virtual and augmented environments,” Proc. DAFx‑15, Trondheim, 2015. PDF
5. J.O. Smith, “Physical Audio Signal Processing: for Virtual Musical Instruments and Audio Effects,” online book, CCRMA, Stanford University, 2010. Link
6. F. Silva, P. Guillemain, J. Kergomard, B. Mallaroni, and A.N. Norris, “Approximation formulae for the acoustic radiation impedance of a cylindrical pipe,” J. Sound Vib., vol. 322, pp. 255–263, 2009. DOI
7. G. Golo, A.J. van der Schaft, P.C. Breedveld, and B.M. Maschke, “Hamiltonian formulation of distributed‑parameter systems with boundary energy flow,” J. Geometry and Physics, vol. 42, pp. 166–194, 2004. DOI
8. G. Haine and D. Matignon, “Structure‑preserving discretisation of port‑Hamiltonian systems: Explicit and implicit schemes for the wave equation,” IFAC‑PapersOnLine, vol. 54, no. 9, pp. 478–483, 2021. DOI
9. V. Mehrmann, R. Morandin, S. Olmi, and E. Schöll, “Port‑Hamiltonian modelling of fluid‑structure interaction in pipes,” Math. Comput. Model. Dyn. Syst., 2025. [TODO: complete reference]
10. C. Zwikker and C.W. Kosten, Sound Absorbing Materials, Elsevier, Amsterdam, 1949. (Classic text on viscothermal boundary‑layer theory)
11. S. Bilbao, Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics, Wiley, Chichester, 2009. DOI
12. J.B. Heywood, Internal Combustion Engine Fundamentals, 2nd ed., McGraw‑Hill, New York, 2018. (Standard reference for engine thermodynamics and gas dynamics)
13. G.P. Blair, Design and Simulation of Four‑Stroke Engines, SAE International, Warrendale, PA, 1999. (Wave‑action modelling of engine manifolds)
14. A. F. Mills, Basic Heat and Mass Transfer, 2nd ed., Prentice Hall, Upper Saddle River, NJ, 1999. (General reference for compressible flow and heat transfer)

Appendix A. Nomenclature

Appendix B. Physical Constants

Appendix C. Typical Engine Parameters

Appendix D. Relationship to Digital Waveguide Networks

Smith (2010) showed that a digital waveguide network (DWN) is exactly equivalent to a staggered‑grid finite‑difference scheme when the Courant number is exactly 1 ($c_0\Delta t / \Delta x = 1$). In that special case the leapfrog update reduces to shifting samples along delay lines – a waveguide.

Our PHS formulation generalises this:

• For $\text{CFL} = 1$, our staggered‑grid scheme reduces to a DWN with scattering junctions at the boundaries.
• For $\text{CFL} < 1$, we retain the full PHS structure with a sparse solver, gaining the ability to use coarser grids and handle frequency‑dependent boundary conditions (like the radiation impedance) without approximating them as pure delays.
• The midpoint rule gives second‑order accuracy in time, whereas the standard waveguide uses first‑order Euler (or trapezoidal integration via scattering).

The Baldan et al. (2015) engine model uses DWNs exclusively. Our approach produces the same sound but with stronger mathematical guarantees and the flexibility to add nonlinear components without sacrificing stability.

Appendix E. Recommended Source File Structure

Scripts/
├── Physics/
│   ├── Port.cs                     # Effort/flow port descriptor
│   ├── PHSVolume.cs                # Base class (acoustic & thermo)
│   ├── PHSPipe1D.cs                # Staggered‑grid 1D duct
│   ├── PHSCylinder.cs              # Cylinder (inherits PHSVolume)
│   ├── PHSEngineNetwork.cs         # Multi‑cylinder assembly
│   ├── OpenEndRadiation.cs         # Reflection filter + monopole
│   ├── Orifice.cs                  # Compressible orifice mass flow
│   ├── FlowNoise.cs                # Turbulence noise source
│   └── GlobalSolver.cs             # Midpoint rule, sparse assembler
├── Tests/
│   ├── TestScenario.cs             # Abstract base class
│   ├── TestRunner.cs               # Sequencer node
│   ├── Test1_PipeImpulse.cs        # Waveguide validation
│   ├── Test2_Helmholtz.cs          # Volume + pipe coupling
│   ├── Test3_OrificeFlow.cs        # Choked/subsonic validation
│   ├── Test4_MotoredCylinder.cs    # Single cylinder, no combustion
│   ├── Test5_FiringCylinder.cs     # With combustion
│   ├── Test6_MultiCylinder.cs      # 4‑cylinder engine
│   └── Test7_FlowNoise.cs          # Noise injection
├── SoundEngine.cs                  # AudioStreamPlayer wrapper
└── phs_implementation_guide.html   # This document

End of document – Version 2.0