High‑pressure exhaust pulses steepen into shock waves that must be resolved cleanly without artificial diffusion that dulls the “crack”. This chapter details the flux‑limiting technique embedded in the staggered‑grid Euler solver. It also explains how both numerical and physical dissipation are integrated into the discrete port‑Hamiltonian structure through the R matrix.
Centred differences produce spurious oscillations near shocks (Gibbs phenomenon). First‑order upwind fluxes smear the shock over many cells, turning a sharp crack into a soft thud. We need a high‑order, oscillation‑free method.
We reconstruct left and right states at a face from the cell‑centred conserved vectors \(\mathbf{q} = (\rho, \rho u, \rho E)^T\) (obtained by dividing the cell state by \(A_i \Delta x\)). The limited slopes are applied to \(\mathbf{q}\):
$$ \begin{aligned} \mathbf{q}_L &= \mathbf{q}_i + \frac{\Delta x}{2}\,\sigma_i, \\ \mathbf{q}_R &= \mathbf{q}_{i+1} - \frac{\Delta x}{2}\,\sigma_{i+1}, \end{aligned} $$where the slope \(\sigma_i\) is the minmod of forward and backward differences:
$$ \sigma_i = \text{minmod}\!\left(\frac{\mathbf{q}_{i+1}-\mathbf{q}_i}{\Delta x},\, \frac{\mathbf{q}_i-\mathbf{q}_{i-1}}{\Delta x}\right). $$The minmod function returns zero if the signs differ, preserving monotonicity. This ensures the TVD property.
The numerical flux at face \(i+1/2\) includes the pre‑computed face area:
$$ \mathbf{F}_{i+1/2} = \frac12 A_{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_{i+1/2} = \max(|u_L|+c_L,\, |u_R|+c_R)\) and the physical flux \(\mathbf{F}(\mathbf{q})\) defined in Chapter 3. This amount of artificial viscosity is minimal but sufficient to stabilise the shock transition.
Because the solver is implicit, the flux limiter is evaluated at the midpoint state \((\mathbf{x}^{n+1}+\mathbf{x}^n)/2\), where \(\mathbf{x} = (\ldots, m_i, p_{\text{mom},i}, E_{\text{tot},i}, \ldots)\). This creates a nonlinear system whose residual involves the numeric fluxes; the Newton‑Krylov solver requires the Jacobian of these fluxes w.r.t. the state variables.
The staggered finite‑volume scheme, including the Rusanov dissipation and the physical source terms, can be cast in the port‑Hamiltonian form \(\dot{\mathbf{x}} = (\mathbf{J} - \mathbf{R})\nabla H + \mathbf{B}\mathbf{u}_{\text{ext}}\). The construction follows the method of Trenchant et al. (2018), where the skew‑symmetric matrix \(\mathbf{J}\) arises from the central (non‑dissipative) part of the flux difference and the source term \(S_{\text{area},i}\), while the positive‑semidefinite \(\mathbf{R}\) matrix collects the numerical dissipation (Rusanov) and the Darcy–Weisbach friction. The wall heat transfer also contributes a dissipative term that can be included in \(\mathbf{R}\) if the thermal port is made explicit; in our implementation it is added as an explicit source with a provably non‑positive effect on \(dH/dt\).
For the purpose of this guide we do not spell out the full algebraic matrices; the C# code and the residual evaluation logic implicitly encode the required sparsity and signs. A developer who follows the assembly procedure described in Chapter 11 will obtain a system whose Jacobian respects the PHS structure. The key guarantee is the discrete power balance \(dH/dt = -\nabla H^T \mathbf{R}\nabla H + \mathbf{y}^T\mathbf{u}_{\text{ext}}\), which holds because each flux and source term is built to satisfy summation‑by‑parts.
The Darcy–Weisbach friction factor \(f\) for smooth pipes is computed with the explicit Swamee–Jain approximation (valid for \(10^{-6} \le \epsilon/D_h \le 10^{-2}\) and \(5\,000 \le \text{Re} \le 10^8\)):
$$ f = \frac{0.25}{\left[\log_{10}\!\left(\frac{\epsilon/D_h}{3.7} + \frac{5.74}{\text{Re}^{0.9}}\right)\right]^2}. $$The Nusselt number for turbulent pipe flow is obtained from the Dittus–Boelter correlation: \(\text{Nu} = 0.023\,\text{Re}^{0.8}\,\text{Pr}^{n}\) (\(n=0.4\) for heating, \(0.3\) for cooling), where \(\text{Pr} = \mu c_p / k_f\) is the Prandtl number. For air at moderate temperatures, \(\text{Pr} \approx 0.71\). Fluid properties \(\mu\) and \(k_f\) can be treated as constants (Appendix B) or evaluated from temperature‑dependent polynomials for high‑precision work.
The Rusanov flux with the wave‑speed estimate \(\alpha = \max(|u|+c)\) is known to be entropy‑stable [3]. This means that the numerical scheme satisfies a discrete entropy inequality, preventing non‑physical expansion shocks.
Comparisons against analytical shock tube problems (Sod’s problem) confirm that the scheme captures the shock front within 2–3 cells while preserving the correct wave speeds and amplitudes. This is essential for producing the correct acoustic impulse in the FW‑H radiator (Chapter 8).