Why the original PVU scheme smears stationary contacts
Kinetic or Boltzmann schemes for inviscid gas dynamics are attractive because the Euler equations appear as moments of the Boltzmann equation when the velocity distribution is Maxwellian, and because the advective terms in the molecular-velocity Boltzmann equation are linear, so upwinding is simpler and the scheme does not need the eigenstructure of the macroscopic system. Traditional kinetic schemes are dynamically consistent, positive (L1 stable), entropy stable, and robust at high Mach numbers, but they produce overly diffusive contact discontinuities in 1D and slip surfaces in 2D, artificially broaden the laminar boundary layer, and suffer spurious pressure and velocity oscillations near isolated contact discontinuities.
The paper starts from the peculiar velocity based upwind (PVU) Boltzmann scheme, which separates fluid and peculiar velocities and applies upwinding to each, yielding a Zha-Bilgen-type convection-pressure splitting at the macroscopic level. The original PVU scheme is robust but overly diffusive and cannot recognize any stationary discontinuity. The authors' modified partial differential equation (MPDE) analysis of the semi-discrete convection step shows why: the numerical diffusion coefficient associated with the peculiar velocity part does not vanish as the Mach number tends to zero, which is the primary reason the scheme diffuses stationary contact discontinuities.
Redesigning peculiar velocity diffusion to switch off at contacts
The Boltzmann equation rewritten in terms of peculiar velocity lets the fluid-velocity flux be treated independently of the peculiar-velocity flux. Upwinding the convective flux and treating the pressure flux centrally at the Boltzmann level would resolve an isolated stationary contact discontinuity exactly, but that discretization is not stable across strong nonlinear discontinuities. The paper therefore keeps the stable fully upwind treatment in general and instead introduces diffusion control parameters into the numerical flux, requiring them to satisfy a dynamical consistency condition so that equal conserved states on both sides give equal parameters.
The goal is a peculiar-velocity numerical diffusion that vanishes as the Mach number goes to zero and stays low but nonzero elsewhere, sharpening linear discontinuities without losing robustness across shocks and expansion fans. Variant 1 sets the advective control to 1 and the peculiar-velocity control to the local sound speed over the interface-averaged sound speed. Variant 2 uses a velocity-ratio advective control and sets the peculiar-velocity control to exp(-c^2/u^2), which is zero when the bulk velocity is zero, so the numerical flux reproduces the exact Riemann solution connecting left and right states through an isolated stationary contact discontinuity.
Variant 3 removes the residual energy flux at contacts
Across a stationary contact discontinuity, variant 1 still produces a nonzero numerical energy flux term proportional to the difference of p a^2 divided by twice the averaged sound speed, which is what diffuses the 1D stationary contact. Variant 3 approximates the density in that nonzero term as p over the squared averaged sound speed, turning the term into one proportional to (p_j - p_{j+1}) times the averaged sound speed. That term shuts off automatically whenever the pressures are equal, so the scheme recognizes the one-dimensional stationary contact discontinuity.
The same idea extends to two dimensions, where the density in the numerical energy flux is approximated as p_K over the squared averaged sound speed. The paper notes that this two-dimensional diffusion does not completely switch off across stationary slip surfaces, so variant 3 resolves them less well than variant 2. The stability condition for all variants follows from a positivity analysis of the compactly supported PVU scheme and takes a CFL-type form involving the normal velocity plus three times the square root of RT, summed over cell faces.
Exact 1D contacts and 2D slip surfaces in benchmarks
The paper reports several 1D and 2D benchmark compressible flow test cases. In the horizontal slip flow test, a Mach 3.0 stream slides over a Mach 2.0 stream in a unit square discretized with 20 by 20 Cartesian cells, with the same density and pressure throughout and supersonic inflow and outflow at left and right. Variant 2 reproduces exactly a grid-aligned, stationary, isolated slip surface free from numerical diffusion. The authors use this case to make a point that matters for readers comparing schemes: a scheme that accurately captures the one-dimensional stationary contact discontinuity need not capture the two-dimensional slip surface with the same precision.
The odd-even decoupling test moves a Mach 6 normal shock through a rectangular channel divided into 800 by 20 square cells with a perturbed centerline, evolving to 55 time units. Both variants remain free from this shock instability, which the paper says causes complete breakdown in some schemes such as HLLC. The positivity test uses large initial Mach numbers with Mx = 10 and My = -3 or 2, pressure 0.1 bar, temperature 35 K, a 0.02 by 0.02 cell size, CFL 0.25, and 100 iterations; both variants preserve pressure and density positivity without breaking down.
Shock tests show robustness is retained
In the shock diffraction test, a Mach 5.09 normal shock diffracts around a 90 degree corner in a unit square divided into 200 by 200 Cartesian cells, with the corner at x = 0.05 and y = 0.45 and the solution evolving to 0.1561 time units. The numerical results are free from entropy-violating solutions and oscillations, whereas the first-order HLLC solution shows an instability ahead of the shock near the upper boundary that becomes more pronounced at higher order.
The double Mach reflection test propagates an oblique Mach 10 shock at 60 degrees to the bottom wall through a four-by-one domain divided into 480 by 160 cells, evolving to 0.2 time units. Neither variant produces an unphysical kink in the Mach stem, an instability known in the HLLC scheme. Repeating the test on a 1920 by 640 mesh with the Barth-Jespersen limiter and primitive-variable reconstruction, variant 2 cleanly captures the complex vortical structure along the primary slip line. In the hypersonic half-cylinder test at Mach 20 on a 160 by 20 mesh, the scheme also avoids the carbuncle instability.
What the scheme does not settle
The evidence is numerical and benchmark-based. The paper demonstrates exact capture of stationary contact discontinuities in 1D and slip surfaces in 2D through test cases, and proves positivity under CFL-type conditions for the compactly supported PVU scheme, but the Maxwellian-based PVU scheme does not inherit that positivity proof because the Maxwellian has infinite support in velocity space. The authors adopt a 3 sqrt(RT) velocity cutoff for the stability condition, noting that integrating the Maxwellian first moment over that interval captures about 99.73 percent of its total mass, while the second through fifth moments cover approximately 99.73, 97.92, 97.43, and 93.2 percent of the total area in Mandal and Deshpande's numerical experiments.
The paper also states that the study did not involve the use of any data and reports no competing interests. The variants are presented as first- and second-order schemes, with the second-order extension using three-stage strong stability preserving Runge-Kutta in time. The reported comparisons are against the original PVU scheme and, in several shock tests, against HLLC behavior described in the paper.
Which variant to prefer, and when
The paper's own priority is clear from its comparisons. Variant 2 is the one that reproduces an exact grid-aligned stationary slip surface in the horizontal slip flow test, because its peculiar-velocity diffusion vanishes when the bulk velocity is zero. Variant 3 is the one that removes the residual numerical energy flux responsible for diffusing the 1D stationary contact, but its two-dimensional diffusion does not fully switch off across stationary slip surfaces, so it resolves them less well than variant 2. The authors conclude that schemes capturing steady two-dimensional slip surfaces exactly are more effective at capturing the negative skin friction coefficient in the shock wave-boundary layer interaction test case.
For a reader deciding what to take from this paper, the practical split is therefore between exactness at stationary linear discontinuities and behavior across strong nonlinear ones. Variant 2 buys exact slip-surface capture at zero bulk velocity; variant 3 buys a cleaner 1D contact by shutting off the energy-flux term when pressures match. Both retain the original scheme's robustness in the odd-even decoupling, positivity, shock diffraction, double Mach reflection, and hypersonic half-cylinder tests, so the choice turns on which discontinuity the application must resolve most sharply rather than on a general accuracy ranking.
Sources used in this answer
Accurate contact discontinuity resolving Boltzmann schemes based on peculiar velocity
MPDE analysis shows the original PVU Boltzmann scheme's peculiar-velocity numerical diffusion does not vanish as Mach number goes to zero, so it cannot recognize stationary contact discontinuities; redesigning that diffusion to depend only on bulk fluid velocity yields exact 1D contact capture and exact 2D slip-surface capture in variant 2, while variant 3 removes the residual energy flux at 1D contacts.