PaperlayerView exact revision on arXiv

Source-linked AI summary

Vectorial lattice Boltzmann solver for compressible inviscid flows with generic equation of state

S. A. Hosseini, I. V. Karlin

arXiv:2609.10639v1math.NAcs.NAphysics.flu-dyn

TL;DR

Accurately and robustly treating high-Mach compressible flows remains a longstanding computational-fluid-dynamics challenge, while classical lattice Boltzmann models face difficulties at high Mach numbers. This paper develops a vectorial lattice Boltzmann solver with adaptive relaxation and time-stepping for Euler flows with generic equations of state. It reports second-order convergence to the Euler limit in the absence of shocks under acoustic scaling and demonstrates the solver on ideal and non-ideal gases, including the BZT regime.

  • Problem

    High-Mach compressible flows remain challenging to treat accurately and robustly, and classical lattice Boltzmann models still face difficulties in that regime.

    Source: §I · p. 1

  • Method

    The paper combines vectorial lattice Boltzmann dynamics with a time-step-scaled adaptive relaxation sensor and conservative adaptive time-stepping for generic-EOS Euler flows.

    Source: §I · p. 1, §V · p. 17

  • Results

    The solver is second-order convergent to the Euler limit without shocks under acoustic scaling and is validated on ideal- and non-ideal-EOS cases, including BZT wave behavior.

    Source: D. Non-ideal (van der Waals) shock… · p. 9, F. BZT (dense-gas) two-dimensional… · p. 12, §V · p. 17

  • Takeaways & Limitations

    The reported tests demonstrate the solver across one- and two-dimensional compressible-flow cases with ideal and non-ideal equations of state.

    Source: §V · p. 17

Abstract

from arXiv · show

We develop a vectorial lattice Boltzmann model with a space-time adaptive relaxation coefficient, combined with adaptive time-stepping for compressible Euler dynamics with generic equation of state. The model, due to special form of the relaxation coefficient devised here is shown to converge to the Euler limit with second-order accuracy in the absence of shocks under acoustic scaling. The solver is robust and able to properly capture compressible gas dynamics for both ideal and non-ideal equations of state, including the Bethe--Zel'dovich--Thompson regime. This is verified and demonstrated through a variety of configurations of increasing complexity.

I. INTRODUCTION

The introduction frames robust high-Mach compressible flow as a longstanding numerical challenge and motivates a kinetic alternative to conventional shock-capturing methods. The paper formulates the target as compressible Euler balance laws for conserved mass, momentum, and energy.

I. INTRODUCTION

  • Classical finite-volume shock-capturing methods remain workhorses for compressible Euler flows, while lattice Boltzmann methods offer a local kinetic route that is harder to extend to high Mach numbers.The paper motivates revisiting vectorial lattice Boltzmann methods for that regime.Source: §I · p. 1

A. Conserved variables

B. Balance equations

  • The Euler system is closed when pressure is specified as a function of thermodynamic state through an equation of state.Its flux is an algebraic function of the conserved vector, making the equations a first-order hyperbolic conservation-law system.Source: B. Balance equations · p. 2

C. Equation of state

The paper leaves the thermodynamic closure generic and constructs a vectorial lattice Boltzmann scheme whose equilibria match the target flux algebraically. This permits EOS changes through the specified flux without re-deriving the scheme’s other ingredients.

C. Equation of state

  • The EOS specifies pressure from density and specific internal energy, subject to thermodynamic consistency.The paper considers ideal-gas and van der Waals closures; the latter can exhibit non-classical BZT wave behavior.Source: C. Equation of state · pp. 2–3

A. Vectorial lattice Boltzmann formulation

  • Unlike scalar-population schemes, VLBM populations are vector-valued and each carries a share of the conserved vector directly.They are transported along a small set of discrete velocities and relax toward equilibrium.Source: A. Vectorial lattice Boltzmann… · p. 3
  • The equilibrium is fixed by exact algebraic consistency with the target flux, rather than by a truncated Maxwellian moment expansion.The equilibrium is exact for any admissible EOS specified through the flux.Source: A. Vectorial lattice Boltzmann… · p. 3
  • Changing the EOS requires no further scheme modification once its flux is specified, and mass, momentum, and total energy are conserved exactly at the discrete level.The conserved vector is invariant under relaxation for any relaxation parameter or EOS.Source: A. Vectorial lattice Boltzmann… · p. 3

B. Hyperbolicity and linear stability

The relaxation system requires a link speed that dominates reachable characteristic speeds to approximate the target hyperbolic system stably. Its relaxation parameter controls numerical dissipation, which vanishes at β = 1.

B. Hyperbolicity and linear stability

  • The link speed must dominate every reachable characteristic speed in each direction for the relaxation system to remain a valid, entropy-dissipative approximation.With one isotropic link speed, the paper chooses a safety margin above the domain-wide maximum of |u · eα| + cs.Source: B. Hyperbolicity and linear stability · pp. 3–4
  • At β = 1 the physical amplification branch is non-dissipative, while β = 1/2 recovers the classical Lax–Friedrichs scheme.The paper states a wavenumber-wise linear stability criterion for the physical branch, but its formula is not included in the supplied passage.Source: B. Hyperbolicity and linear stability · p. 4
  • For 0 < β ≤ 1 under the subcharacteristic condition, the leading direction-split dissipative correction is positive semidefinite and vanishes at β = 1.The scheme is consistent with the target Euler system to O(δt).Source: B. Hyperbolicity and linear stability · p. 4

D. Relaxation coefficient and shock sensor

A constant relaxation parameter other than β = 1 introduces a Navier–Stokes-level correction under acoustic scaling, making the scheme first-order accurate for the Euler limit. The proposed adaptive sensor scales relaxation with time-step size and gates it to compression-dominated regions.

D. Relaxation coefficient and shock sensor

  • For β ≠ 1, the acoustic-scaling deviation appears at the Navier–Stokes level, so the scheme is first-order accurate and its dissipation is not a shock-capturing feature.The paper identifies this as the motivation for adapting β.Source: D. Relaxation coefficient and shock… · p. 4

C. Equivalent macroscopic equations

  • The proposed β form approaches 1 as δt → 0 at first order, pushing the residual dissipative contribution to order three and above.The function g is constrained to be nonnegative.Source: C. Equivalent macroscopic equations · p. 4
  • The sensor converts componentwise flux-divergence estimates into a dimensionless activity measure, smooths it, then gates it with a dilatation–vorticity discriminator.The discriminator is used because flux gradients alone can be large in both shocks and smooth, strongly rotational flow.Source: C. Equivalent macroscopic equations · pp. 4–5
  • The gate tends toward one in compression-dominated regions and zero in rotation- or shear-dominated regions; in one dimension, it is identically one.Thus the gate has no effect on the one-dimensional sensor.Source: C. Equivalent macroscopic equations · p. 5
  • The two sensor variants are used without further modification, with the gated variant suppressing relaxation activity due to smooth rotation while leaving it unchanged at genuine shocks.The paper sets βmin to the stable constant-β baseline for each case and βmax = 1.Source: C. Equivalent macroscopic equations · p. 5

E. Adaptive time step

The solver retargets link speed and time step as the domain’s maximum signal speed changes, balancing stability margin against numerical dissipation. When the link speed changes, populations are reprojected while preserving their independent moments.

E. Adaptive time step

  • Retargeting link speed and time step avoids a fixed initial speed becoming under-resolving as signals accelerate or needlessly dissipative as they slow.The fixed-speed choice is based on initial data and may no longer match the fastest signal speed during the run.Source: E. Adaptive time step · p. 5
  • Each retargeting recomputes the current maximum signal speed and sets link speed using a target CFL that trades safety margin against dissipation.Smaller CFL provides a larger margin above the current maximum speed but increases dissipation; retargeting less often relaxes the effective margin between updates.Source: E. Adaptive time step · p. 5
  • Changing link speed requires rebuilding populations while preserving all independent moments, since resetting them to equilibrium would discard non-equilibrium content.For D2Q4, the moments are W, directional flux-moments J1 and J2, and the second-order moment K.Source: E. Adaptive time step · pp. 5–6
  • The reprojection formulas recover the equilibrium closure when K = 0 and Jα = Qα(W), and reduce to the D1Q2 formula in one dimension.Retargeting holds J1, J2, and K fixed while re-expressing the populations at the new link speed.Source: E. Adaptive time step · p. 6

IV. NUMERICAL APPLICATIONS

The validation proceeds through benchmarks of increasing difficulty, from exact steady solutions and one-dimensional tests to non-ideal shock tubes and two-dimensional supersonic flow. It begins with a case without genuine discontinuities to assess the gradient-adaptive sensor apart from shock-capturing demands.

IV. NUMERICAL APPLICATIONS

  • The validation spans exact steady solutions, one-dimensional Riemann problems, composite waves, non-ideal shock tubes, and two-dimensional strongly supersonic flow.Source: §IV · p. 6
  • The benchmark sequence includes van der Waals and BZT dense-gas shock tubes alongside ideal-gas cases.Source: §IV · p. 6
  • The sequence opens with no genuine discontinuity so the gradient-adaptive sensor can be assessed before shock-capturing demands are introduced.Source: §IV · p. 6

A. Static vortex convergence: the Gresho problem

The Gresho vortex tests sensor behavior on a smooth, exact steady solution, while the Shu–Osher problem tests shock capture alongside a smooth post-shock wave train. The sensor yields close-to-second-order error scaling for the vortex and reduces Shu–Osher density errors relative to the constant-relaxation baseline.

A. Static vortex convergence: the Gresho problem

  • The Gresho vortex is an exact, time-independent compressible Euler solution, so departures from its initial state measure numerical error.Source: A. Static vortex convergence: the… · p. 6
  • For the Gresho vortex, the sensor-equipped solver exhibits close-to-second-order error scaling against the best stable constant-β scheme.The comparison uses the L2 error of |u| versus grid spacing for Sensor B and fixed β = 0.995.Source: Figure 1 · p. 7

B. Shu–Osher shock/entropy-wave interaction

  • At t = 1.9, the Shu–Osher shock is captured comparably by the sensor and constant-β baseline, while the sensor largely removes the baseline’s post-shock wave-train amplitude deficit.Figure 2 shows the density profile at N = 5000 against the reference.Source: B. Shu–Osher shock/entropy-wave… · p. 7, Figure 2 · p. 7
  • The sensor reduces Shu–Osher L2(ρ) errors by at least a factor of two, although both configurations converge at a close-to-first-order rate because shock-front error dominates.Source: B. Shu–Osher shock/entropy-wave… · p. 7

C. Strong shock tube

The sensor-driven solver reproduces severe ideal-gas shock structures and tracks independent references across van der Waals shock tubes, including non-classical wave behavior. In the Shu–Osher problem, sensing improves the post-shock wave train while convergence remains close to first order because shock-front errors dominate.

C. Strong shock tube

  • The Mach ≈198 ideal-gas shock tube reproduces the expected rarefaction, contact, compressed region, and terminating shock while tracking the independent reference closely.Relative errors are 0.2%, 0.2%, and 0.1% of the reference maxima for the three fields; errors concentrate around the contact and shock.Source: C. Strong shock tube · pp. 7–8

D. Non-ideal (van der Waals) shock tubes: the Argrow benchmarks

  • At N=1000, L1(ρ)/L1(P) errors were 1.31×10−3/4.90×10−4 (Case I), 6.91×10−4/1.20×10−3 (Case II), and 1.96×10−3/2.52×10−3 (Case III), with profiles tracking the reference.The sensor dips at genuine discontinuities, including Case I’s expansion rarefaction shock, and recovers near βmax = 1 across smooth fans and plateaus.Source: D. Non-ideal (van der Waals) shock… · p. 9
  • In Case I, Γ < 0 reverses the classical wave pattern: compression remains smooth while expansion forms a rarefaction shock; Case III instead has a classical compound-wave structure.Case II remains entirely classical, with an ordinary shock and rarefaction fan.Source: D. Non-ideal (van der Waals) shock… · p. 9

E. Two-dimensional Riemann problems

The D2Q4 solver is tested on classical two-dimensional Riemann interactions and a dense-gas BZT problem with inverted wave behavior. It closely matches independent references, including for the BZT case where Γ is negative across nearly the entire domain.

E. Two-dimensional Riemann problems

  • The ideal-gas tests exercise oblique shocks, diagonal shock interactions, and contact-driven roll-up before outer waves reach the domain boundary.Configuration 12 combines contacts, a shock, and a rarefaction; Configuration 3 has shocks on all four edges.Source: E. Two-dimensional Riemann problems · pp. 9–10
  • At N=800, Configurations 12 and 3 have relative L1 density errors of 0.15% and 0.39%, reproducing their major shocks, contacts, and Kelvin–Helmholtz structures.The reference uses a two-dimensional unsplit Godunov-type finite-volume solver at N=1000.Source: E. Two-dimensional Riemann problems · p. 10, Figure 10 · p. 11, Figure 11 · p. 11

F. BZT (dense-gas) two-dimensional Riemann problem

  • The dense-gas setup reuses Configuration 12’s quadrant topology, but its classical shock and rarefaction assignments invert when Γ < 0.The two non-contact edges were independently checked as one-dimensional mini-Riemann problems before assembling the 2D initial condition.Source: F. BZT (dense-gas) two-dimensional… · p. 11
  • The BZT 2D solution reproduces a smooth compression fan and rarefaction shock with 0.12% relative L1 density error; Γ < 0 across 99.0% of the domain.The VLBM and reference also agree on the small classical-fluid pocket behind the inverted-wave intersection.Source: F. BZT (dense-gas) two-dimensional… · pp. 11–12

G. Converging BZT compression wave

The section studies whether geometric focusing can drive a BZT compression wave across Γ = 0, then compares its dynamics with an ideal gas and sweeps pulse strength. The 2D VLBM reproduces the 1D reference and reveals a sharp, thermodynamic transition that density alone does not show.

a. Governing non-dimensional number and test-case design

  • Geometric amplification can push a flow that starts deep in the BZT regime across Γ = 0, unlike the ideal-gas construction, for which such a crossing is impossible.The 4-quadrant test exhibits BZT wave-type inversions but stays close to its initial isentropes, motivating the converging-wave test.Source: G. Converging BZT compression wave · p. 12
  • The 1D cylindrical reference shows focusing raises the axis density from 0.70 to ≈1.19 and Γ from −0.26 to +3.9 near t ≈ 1.4, before reflection produces a rarefaction.The pulse is initialized as an inward-pointing Gaussian velocity annulus in a BZT van der Waals gas; M0 controls focusing strength once geometry and ambient state are fixed.Source: G. Converging BZT compression wave · p. 12, Figure 14 · p. 14

b. 2D Cartesian D2Q4 VLBM results

c. Contrast with a classical (ideal) gas

  • With matched geometry, amplitude, and ambient sound speed, the BZT pulse peaks later and higher than the ideal-gas pulse: t ≈ 1.4, ρpeak ≈ 1.10 versus t ≈ 1.1, ρpeak ≈ 0.93.The BZT sound speed decreases under compression (0.450 → 0.246), delaying characteristic catch-up; the ideal-gas sound speed increases (0.450 → 0.501).Source: G. Converging BZT compression wave · p. 13

d. Parametric threshold in pulse Mach number

  • Γpeak rises from ≈−0.3 to ≈+0.75 between M0 = 0.095 and M0 = 0.105, marking a sharp threshold for transiently crossing into Γ > 0.Below the crossing, Γpeak first falls to ≈−0.98 near M0 ≈ 0.07; above it, Γpeak peaks near 4.68 at M0 ≈ 0.17 and settles to ≈4.2 by M0 = 0.35, while focus density varies smoothly through the crossing.Source: G. Converging BZT compression wave · pp. 13–15

H. Inviscid shock-vortex interaction

The shock-vortex test combines a persistent shock and vortex in ideal-gas and BZT media, and cross-checks both against an independent finite-volume method. The conclusion summarizes the generic-EOS Euler solver and identifies dissipation and grid refinement as future extensions.

H. Inviscid shock-vortex interaction

  • The BZT shock raises Γ from −0.255 upstream to +4.47 downstream, confining Γ < 0 fluid to the unshocked upstream region and the vortex’s outer envelope.The BZT and ideal-gas cases use matched shock Mach number Ms = 1.2, with longer runtime and higher resolution for the lower-sound-speed BZT case.Source: H. Inviscid shock-vortex interaction · pp. 15–16
  • Relative L1 density errors against the independent MUSCL+Rusanov reference are 0.093% for the ideal gas and 0.114% for BZT; both runs show the vortex crossing, shock bulging, and a shed upstream disturbance.The residuals are concentrated along the shock front and vortex core.Source: H. Inviscid shock-vortex interaction · p. 16

V. CONCLUSION

  • The proposed VLBM targets compressible Euler flow with a generic single-phase EOS, using gradient-based shock-capturing relaxation and conservative adaptive time-stepping.The validation covers one- and two-dimensional cases with ideal and non-ideal equations of state.Source: §V · p. 17
  • Navier–Stokes-level dissipation and adaptive grid refinement are identified as future extensions, not components of the presented solver.Source: §V · p. 17

Appendix A: Multi-scale analysis

The multiscale analysis expands the discrete evolution under acoustic scaling, first recovering the target Euler balance equations and then examining the next order. With the specified relaxation-parameter scaling and positive g, it establishes second-order convergence to the Euler limit.

Appendix A: Multi-scale analysis

  • The expansion starts from the discrete time-evolution equations, Taylor-expands the shift, and inserts the multiscale ansatz with δt as the smallness parameter.Source: Appendix A · p. 18
  • Summing the order-δt equations recovers the targeted Euler balance equations.Source: Appendix A · p. 18
  • Under acoustic scaling, setting 1/(2β) = 1/2 + δt g(x,t) with positive-definite g yields second-order convergence to the Euler limit.The order-δt² equations are simplified using the preceding order and then summed.Source: Appendix A · p. 18
Loading 2609.10639v1…