Every pipe in the engine model (header, runners, muffler, tailpipe) is simulated with the full 1‑D compressible Euler equations. No linearisation is performed, ensuring accurate reproduction of shock waves, temperature‑dependent sound speed, and mean flow advection even at pressure ratios above 2.2. The pipe is a distributed storage element—its geometric volume contributes to the total gas capacity of the network.
For a duct of cross‑sectional area \(A(x)\), the conservation laws are [1,2]:
$$ \begin{aligned} \frac{\partial (\rho A)}{\partial t} + \frac{\partial (\rho u A)}{\partial x} &= 0, \\ \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^2 A + p A)}{\partial x} &= p\frac{dA}{dx} - F_{\text{fric}}, \\ \frac{\partial (\rho E A)}{\partial t} + \frac{\partial \bigl((\rho E + p) u A\bigr)}{\partial x} &= u F_{\text{fric}} + Q_{\text{wall}}, \end{aligned} $$where \(E = e + \frac12 u^2\) is total specific energy, \(e = c_v T\) internal energy, and \(p = \rho R_g T\) (ideal gas). The source terms are:
The friction factor \(f\) is obtained from the explicit Swamee–Jain approximation (Chapter 7), and the Nusselt number \(\text{Nu}\) for heat transfer from the Dittus–Boelter correlation (Chapter 7).
To faithfully reproduce sound up to the Nyquist frequency \(f_{\max} = 22\,050\text{ Hz}\) at a 44.1 kHz sample rate, the cell length \(\Delta x\) must resolve the shortest wavelength. Using a minimum of \(K = 10\) cells per wavelength and the speed of sound \(c \approx 343\text{ m/s}\) (or the local computed value), we set:
$$ \Delta x \le \frac{c}{K f_{\max}}. $$For a pipe of total length \(L\), the number of cells is \(N = \lceil L / \Delta x \rceil\) and a uniform \(\Delta x = L/N\) is adopted. For highly nonlinear waves (shocks) a finer grid may be required; the chosen \(\Delta x\) provides a baseline for linear‑to‑moderate amplitudes.
The pipe is divided into \(N\) cells of length \(\Delta x\). All conserved quantities are stored at cell centres:
Face areas \(A_{\text{face},i+1/2}\) are pre‑computed via harmonic mean of adjacent cell areas. Face velocities are not stored independently; they are obtained from the cell‑centred state when needed.
Density, velocity, internal energy and pressure are derived from the state:
\[ \rho_i = \frac{m_i}{A_i \Delta x},\quad u_i = \frac{p_{\text{mom},i}}{m_i},\quad e_i = \frac{E_{\text{tot},i} - \frac12 p_{\text{mom},i}^2 / m_i}{m_i},\quad p_i = (\gamma-1)\rho_i e_i. \]// C# state (per cell) double[] m; // mass double[] p_mom; // momentum double[] E_tot; // total energy // geometry (constant) double[] cellArea; double[] faceArea; double dx;
For cell \(i\), the ODEs are directly obtained from the integral form of the governing equations:
$$ \begin{aligned} \frac{d m_i}{d t} &= -(\Phi_{m,i+1/2} - \Phi_{m,i-1/2}), \\ \frac{d p_{\text{mom},i}}{d t} &= -(\Phi_{\text{mom},i+1/2} - \Phi_{\text{mom},i-1/2}) + S_{\text{area},i} - F_{\text{fric},i}, \\ \frac{d E_{\text{tot},i}}{d t} &= -(\Phi_{E,i+1/2} - \Phi_{E,i-1/2}) + u_i F_{\text{fric},i} + Q_{\text{wall},i}, \end{aligned} $$where \(\Phi_m, \Phi_{\text{mom}}, \Phi_E\) are the numerical fluxes (mass, momentum, energy) through the face.
Integrating the pressure force over the cell gives the exact finite‑volume expression
$$ S_{\text{area},i} = p_i \frac{A_{\text{face},i+1/2} - A_{\text{face},i-1/2}}{\Delta x}. $$This uses the pre‑computed face areas and requires only the cell‑centred pressure \(p_i\). It is second‑order accurate for smooth area changes and conserves total momentum exactly.
Boundary handling: At the first cell (\(i=0\)), \(A_{\text{face},-1/2}\) is the area of the left boundary face; for a solid wall the corresponding pressure force is zero because the wall exerts no net force in the flow direction. At the last cell (\(i=N-1\)) the right boundary face area is used analogously.
We use a Rusanov (local Lax‑Friedrichs) flux with TVD minmod limiter. For each face, left and right states \(\mathbf{q}_L, \mathbf{q}_R\) are reconstructed from the cell‑centred conserved variables per unit volume \(\mathbf{q} = (\rho, \rho u, \rho E)^T\) using slope‑limited (minmod) gradients.
$$ \mathbf{F}_{i+1/2} = \frac12 A_{\text{face},i+1/2}\bigl[ \mathbf{F}(\mathbf{q}_L) + \mathbf{F}(\mathbf{q}_R) - \alpha_{i+1/2}(\mathbf{q}_R - \mathbf{q}_L) \bigr], $$with \(\alpha = \max(|u|+c)\) and the physical flux
$$ \mathbf{F}(\mathbf{q}) = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ (\rho E + p) u \end{pmatrix}. $$This flux provides both the hyperbolic transport and the numerical dissipation required for shock stability, applied identically to all three conserved variables.
void ComputeFaceFlux(int i, double[] qL, double[] qR, double alpha, double gamma,
out double fluxMass, out double fluxMom, out double fluxEnergy)
{
double A = faceArea[i];
double fRho_L = qL[1];
double fRhoU_L = qL[1]*qL[1]/qL[0] + (gamma-1)*(qL[2] - 0.5*qL[1]*qL[1]/qL[0]);
double fRhoE_L = (qL[2] + (gamma-1)*(qL[2] - 0.5*qL[1]*qL[1]/qL[0])) * qL[1]/qL[0];
// ... same for R ...
double avgMass = 0.5 * A * (fRho_L + fRho_R);
double avgMom = 0.5 * A * (fRhoU_L + fRhoU_R);
double avgEnergy= 0.5 * A * (fRhoE_L + fRhoE_R);
double diss = 0.5 * A * alpha;
fluxMass = avgMass - diss * (qR[0] - qL[0]);
fluxMom = avgMom - diss * (qR[1] - qL[1]);
fluxEnergy= avgEnergy- diss * (qR[2] - qL[2]);
}
A pipe end connects to a volume or junction via a PHS port. The mass flow is \(\dot{m} = \rho_{\text{face}} u_{\text{face}} A_{\text{face}}\). The face density is taken from the adjacent cell (\(\rho_{\text{face}} = \rho_0\) or \(\rho_{N-1}\)), and the face velocity is obtained from the cell momentum. This zero‑order extrapolation is consistent with the first‑order accuracy of the boundary flux; a more precise characteristic‑based boundary condition can be implemented if required.
The boundary pressure effort is imposed by setting the ghost‑cell state accordingly (for a solid wall, the ghost cell mirrors the interior pressure and sets velocity to zero).
The friction force \(F_{\text{fric},i}\) and the heat‑transfer term \(Q_{\text{wall},i}\) are added explicitly to the ODEs. Both are strictly dissipative:
These terms can be formally incorporated into the discrete \(\mathbf{R}\) matrix; for implementation they remain as source terms that are easily verified to preserve the inequality \(dH/dt \le \mathbf{y}^T \mathbf{u}_{\text{ext}}\).
The numerical dissipation of the Rusanov flux contributes a symmetric positive‑semidefinite block to \(\mathbf{R}\) that couples neighbouring cells. Its precise form follows from the difference of the diffusive fluxes and is detailed in Trenchant et al. (2018). For a cell \(i\), the local contribution from the friction term is a diagonal entry \(r_{\text{fric},i} u_i\) in the momentum equation; the corresponding energy source \(u_i F_{\text{fric},i}\) is added explicitly. The overall discrete system satisfies \(dH/dt = -\nabla H^T \mathbf{R} \nabla H + \mathbf{y}^T \mathbf{u}_{\text{ext}}\) by construction of the skew‑symmetric flux split and the dissipative sources.
The entire state vector is advanced with the implicit midpoint rule. The nonlinear system is solved by Newton‑Krylov iteration (Chapter 11). Because the pipe stores mass, momentum, and internal energy in every cell, its geometric volume contributes to the total gas capacity; after transients the pipe reaches the same equilibrium pressure as connected tanks.