PaperlayerView exact revision on arXiv

Source-linked AI summary

A relaxation system for the low-Mach limit of kinetic equations: Stability and higher order asymptotic preserving scheme

Giacomo Dimarco, Axel Klar, Theresa Köfler, Lorenzo Pareschi, Sudarshan Tiwari, Yizhou Zhou

arXiv:2610.00811v1math.NAcs.NA

TL;DR

Low-Mach kinetic simulation is computationally demanding because of velocity-space dimensionality and stiffness. The paper derives a reduced hyperbolic moment relaxation system and a high-order AP scheme; experiments report close agreement with full BGK references and substantial computational speedups in stiff regimes.

  • Problem

    Low-Mach kinetic equations are costly to solve because velocity-space discretization is high-dimensional and the fluid-limit equations are stiff.

    Source: §1 · p. 2

  • Method

    The paper derives a reduced hyperbolic moment system from a micro-macro decomposition and combines IMEX Runge-Kutta time integration with high-order spatial reconstructions in an AP framework.

    Source: §1 · p. 3, §8 · p. 29

  • Results

    In continuum-dominated tests, the relaxation model agrees closely with full BGK references, with comparisons reaching relative scales of O(10^-14), while stiff-regime simulations achieve speedups of several orders of magnitude.

    Source: §8 · p. 29

  • Takeaways & Limitations

    As ε tends to zero, the scheme becomes a consistent projection method for incompressible Navier-Stokes, and the continuous moment approach applies in arbitrary spatial dimensions without velocity-space discretization.

    Source: §8 · p. 29

    For larger ε, truncation of the moment hierarchy limits the model’s ability to capture strong non-equilibrium effects, producing structural model error.

    Source: §7.3.4 · p. 29

Abstract

from arXiv · show

This work introduces a hyperbolic relaxation system designed for the simulation of kinetic equations in the low-Mach number limit. The methodology is built upon a micro-macro decomposition of the scaled BGK model, which reformulates the kinetic distribution function into a coupled system consisting of a macroscopic equilibrium part and a microscopic non-equilibrium remainder. By projecting the microscopic deviations onto a set of orthogonal polynomials, we derive a closed moment relaxation system. The resulting system is a version of Grad's 13 moment system with a linear hyperbolic part and a relaxation adapted to the incompressible limit. We prove the model's structural stability in the incompressible Navier-Stokes limit. Moreover, we develop a high order Asymptotic-Preserving (AP) numerical framework using Implicit-Explicit (IMEX) Runge Kutta schemes for temporal accuracy and finite difference WENO reconstructions as well as central difference approximations for high order spatial resolution. The proposed scheme ensures uniform stability and consistency across different physical regimes, automatically degenerating into a consistent high order discretization of the incompressible thermal Navier-Stokes limit as the scaling parameter vanishes. Numerical experiments in one and two dimensions corroborate the theoretical findings, demonstrating significant reductions in computational cost and robust performance across a wide range of Knudsen and Mach numbers.

1. Introduction

Low-Mach kinetic simulations face costly velocity-space discretization and stiffness, while moment closures must preserve physical properties. This work derives a reduced relaxation system and an AP scheme tailored to the incompressible limit.

1. Introduction

  • Velocity-space grids and the moment-closure problem make kinetic simulation computationally burdensome, especially in multidimensional settings.AP methods address stiffness across regimes, but this model is tailored to small ε and is not designed for accurate highly rarefied-gas dynamics.Source: §1 · pp. 2–3
  • Micro-macro decomposition separates equilibrium and non-equilibrium dynamics, allowing pressure to emerge as an incompressibility multiplier while retaining viscous and thermal dissipation.The introduction identifies relatively few numerical works on the incompressible Navier-Stokes limit of BGK, including and [20].Source: §1 · p. 3
  • IMEX Runge-Kutta methods treat stiff relaxation terms implicitly and non-stiff transport terms explicitly, while WENO or central differences provide high-resolution spatial discretization.The implicit treatment avoids excessively small time steps, and the explicit transport treatment avoids high-dimensional nonlinear algebraic solves.Source: §1 · p. 3
  • The paper derives a reduced moment system for arbitrary spatial dimensions and develops an arbitrarily high-order AP scheme using IMEX Runge-Kutta and WENO or central-difference discretizations.In the asymptotic limit, the methods become high-order projection methods for incompressible Navier-Stokes on a collocated grid.Source: §1 · p. 3
  • The paper presents numerical tests in both 1D and 2D, comparing its relaxation model with full kinetic simulations to assess accuracy and computational efficiency.Source: §1 · p. 4

2. Kinetic equations and the low-Mach limit

The paper scales the BGK model around a global equilibrium, identifies the incompressible Navier-Stokes limit, and derives a closed hyperbolic moment relaxation system by projecting macro- and microscopic components.

2. Kinetic equations and the low-Mach limit

  • The derivation starts from small perturbations around a global equilibrium with unit density, zero mean velocity, and unit temperature, yielding a scaled equation for the perturbation.The macroscopic fluctuations are related to density, velocity, and temperature through the stated scaling ansatz.Source: §2 · p. 4
  • As ε approaches zero, the perturbation approaches local equilibrium, and the kinetic limit yields incompressible Navier-Stokes equations with Boussinesq coupling for density and temperature.The macroscopic quantities are defined through velocity moments of the distribution.Source: §2 · pp. 4–5

2.1. Derivation of the relaxation system

  • A micro-macro decomposition splits the distribution into its projected hydrodynamic component and kinetic fluctuation, with the rescaled microscopic variable defined using (I − P)f = O(ε).The derivation applies P and I − P to the scaled kinetic equation and uses the orthogonality relations stated for the collision operator and remainder.Source: §2 · pp. 4–5, §2.1 · pp. 6–7
  • The closure neglects ε(I − P)(vg) and approximates polynomial moments of the remainder by their ε → 0 limits to obtain equations for Θ and Γ.Θ and Γ are defined as moments of the rescaled microscopic variable against polynomials A(v) and B(v).Source: §2 · p. 5, §2.1 · p. 7
  • The resulting system is a hyperbolic moment model for ρ, u, T, Θ, and Γ, with 4, 8, and 13 moments in 1D, 2D, and 3D, respectively.Its hyperbolicity follows from symmetrizability; the reported eigenvalue structure differs across dimensions.Source: §2.1 · p. 7
  • In higher dimensions, the closure differs from a full Hermite expansion and omits terms, but the authors show it remains appropriate for the incompressible Navier-Stokes limit at sufficiently small ε.The hyperbolic main part is equivalent to the linear Grad moment system in 1D, but not in higher dimensions.Source: §2.1 · pp. 7–8

3. Stability of the relaxation system

The stability analysis establishes structural stability of the relaxation system as ε tends to zero, while distinguishing that limit from a different scaling regime where the conditions may fail.

3. Stability of the relaxation system

  • The stability discussion concerns the limit ε → 0, with the relaxation variables rescaled before the system is analyzed.The variables Θ and Γ are rescaled by ε in the equations.Source: §3 · p. 8
  • The relaxation system satisfies the stated structural stability conditions, including a symmetric positive-definite symmetrizer and conditions on the source term.The analysis also verifies invertibility of the specified derivative with respect to UII.Source: §3 · pp. 10–11
  • Convergence to incompressible Navier-Stokes requires additional assumptions on initial data and asymptotic-expansion regularity that the paper does not discuss.The cited framework’s assumptions H4–H6 are identified as the missing technical conditions.Source: §3 · p. 11
  • For ε = 1 as τ tends to zero, the model’s stability conditions are generally not fulfilled, and its formal limit need not be hyperbolic or equal to nonlinear Euler equations.This is a distinct limit from the ε → 0 regime analyzed in the paper.Source: §3 · p. 11

4. Asymptotic-preserving algorithms - Time discretization

The section develops AP time discretizations for the relaxation system, treating low-Mach and small-Knudsen stiffness implicitly. It also gives a variant that treats limiting diffusion implicitly to address parabolic stiffness.

4.2. Additional implicit discretization treating the parabolic stiffness

  • The additional variant treats the diffusive terms in the limiting equations implicitly, yielding a projection scheme with implicit diffusion.This variant adds an implicit discretization of diffusive terms to account for parabolic stiffness.Source: §4 · p. 11, §4.2 · p. 15

4.1. Semi-implicit discretization for low Mach and small Knudsen numbers

  • The first-order AP discretization treats low-Mach and small-Knudsen stiffness implicitly and produces a projection scheme for incompressible thermal Navier-Stokes with explicit diffusion.The section also presents an alternative projection scheme with implicit treatment of diffusion in the limit.Source: §4.1 · pp. 12, 15, §4.2 · p. 15
  • The initial algorithm advances the relaxation system’s state variables through staged updates of the moment variables, velocity, temperature, and pressure.The algorithm takes the state at time t_n as input and computes updated variables at the next time step.Source: §4.1 · pp. 13–14

5. AP scheme - Higher order time-discretizations

The section formulates higher-order time discretizations through IMEX Runge-Kutta methods, pairing explicit and implicit tableaux. It specifies tableau classes and accuracy properties used to characterize these schemes.

5. AP scheme - Higher order time-discretizations

  • IMEX Runge-Kutta methods combine explicit and implicit tableaux to balance time-step stability against the complexity of fully implicit solves.The explicit part handles non-stiff terms, while the implicit part handles stiff terms.Source: §5 · p. 16
  • The explicit tableau is lower triangular with zero diagonal, while the implicit tableau is a diagonally implicit Runge-Kutta method.The tableaux’s weights combine the internal stages.Source: §5 · p. 16
  • Type A methods have an invertible implicit matrix, whereas CK methods use an invertible submatrix and ARS is the CK subclass with a zero initial vector.These definitions distinguish the implicit tableau structures considered.Source: §5 · p. 16

If additionally the explicit method satisfies

The section applies IMEX Runge-Kutta methods to the moment system and states a conditional asymptotic-preserving result. Under the theorem’s tableau and initial-velocity conditions, the scheme is consistent with incompressible Navier-Stokes as ε vanishes.

If additionally the explicit method satisfies

  • The IMEX stage equations separate implicit and explicit components, and the explicit method’s weights satisfy the stated stiff-accuracy condition.The final algorithm for solving the moment system is given as Algorithm 2.Source: If additionally the explicit method… · pp. 16–17
  • With a type A or CK/ARS IMEX method satisfying GSA, and divergence-free initial velocity, the scheme becomes a consistent discretization of incompressible Navier-Stokes as ε → 0.The GSA solution coincides with the last internal stage.Source: If additionally the explicit method… · pp. 16–17
  • As ε vanishes, dominant relaxation terms slave higher-order kinetic moments to macroscopic variables, yielding the macroscopic limit equations.The moment reduction mechanism is the stated link between the relaxation system and its limiting equations.Source: If additionally the explicit method… · p. 17
  • The stated uniform numerical stability claims remain bounded by the established CFL and dissipation constraints, despite the scheme’s high-order accuracy.The passage cautions against interpreting the result as unconditional uniform stability.Source: If additionally the explicit method… · p. 17

6. Numerical Implementation and Spatial Discretization Details

The implementation combines IMEX Runge-Kutta time stepping with finite-difference and WENO spatial operators. The collocated formulation uses local numerical dissipation to suppress pressure oscillations, while explicit diffusion retains a parabolic time-step restriction.

6. Numerical Implementation and Spatial Discretization Details

  • Algorithm 2 uses ARS-type IMEX Runge-Kutta tableaus with GSA properties to treat source-term stiffness efficiently.Its time-step size is independent of the microscopic scaling parameter ε.Source: §6 · p. 18, §6.1 · p. 19

6.1. Spatial Discretization and Time-Step Restrictions

  • The time step follows an advective CFL condition independent of ε, while explicit diffusion in the limit retains a parabolic restriction Δt ≤ O(Δx^2).The AP formulation removes stiffness associated with the relaxation source terms and ε, but not the diffusive stability constraint.Source: §6.1 · p. 19, §6.3 · p. 19
  • Second- or fourth-order central differences handle first derivatives and the discrete Laplacian, while LF fluxes with WENO3/WENO5 reconstruct the Θ flux.A generalized least-squares finite-difference approach for arbitrary grids was also implemented and yielded similar results.Source: §6.1 · p. 19

6.2. Collocated Grid Formulation and Oscillation Suppression

  • On a collocated grid, compact stencils and LF-flux dissipation suppress checkerboard pressure modes by coupling adjacent cells.All primitive variables share identical grid nodes; separate grid staggering is not required for the stated pressure-velocity coupling.Source: §6.2 · p. 19

6.4. Code Availability

  • The implementation scripts and numerical routines are openly accessible in the project’s GitHub repository.The passage identifies the repository as https://github.com/sudarshan825/.Source: §6.4 · p. 19

7. Numerical results

The numerical tests compare the relaxation system with BGK solutions in one and two dimensions. They examine accuracy across parameter regimes, computational cost, convergence, and higher-order behavior.

7.1. Numerical methods for the BGK equation

  • The classical BGK comparison solver [55] uses an ε-dependent time step, while the AP BGK scheme [20] avoids that dependence and tends to a projection scheme in the limit.The [20] method computes the full kinetic evolution.Source: §7.1 · p. 20

7.2. Results in 1D

  • For ε = 0.1 and 0.05, the relaxation-model approximation is much better, with only small deviations from the kinetic solution; at ε = 1, pronounced wave patterns cause its approximation quality to deteriorate.Figures 1–3 compare density, velocity, and temperature at final time tf = 0.1, for ε = 1, 0.1, and 0.05 and τ = 0.1 or 0.01.Source: §7.2 · p. 20, Figure 1 · p. 21, Figure 2 · p. 21
  • The AP relaxation scheme's computation is independent of ε and velocity-grid size Nv, whereas conventional CFL-limited BGK cost grows as ε decreases.The reported speedup applies specifically against the standard explicit BGK solver [55]; no timing comparison is made with implicit AP BGK solvers such as [20].Source: §7.2 · pp. 20, 22

7.3. 2D: Doubly periodic double shear layer

  • The 2D tests use a periodic double shear layer and compare vorticity contours from the relaxation system with BGK reference solutions.The contour field is vorticity computed from the velocity components.Source: §7.3 · p. 23, Figure 4 · p. 25

7.3.1. Accuracy

  • The convergence study confirms expected orders only for ε=10^-6 and τ=0; its refinement time step scales as Δt(h)=O(h^2).Errors are measured against a fine-grid reference, and these results do not establish universal high-order accuracy across physical regimes.Source: §7.3.1 · p. 23

7.3.2. Comparison to BGK

  • Across Figures 4–6, the relaxation-system solution generally agrees well with the BGK reference solver [55] for ε=0.1 and for ε=0.5, τ=0.01, while larger deviations occur at ε=0.5, τ=0.001.Source: §7.3.2 · pp. 23–24

7.3.3. Comparison of higher order results

  • At ε=10^-6 and τ=0, the higher-order IMEX RK3/WENO5 scheme preserves steeper gradients and higher vorticity peaks than IMEX RK2/WENO3, with less numerical diffusion.For ε=0.01, the relaxation solution also fits the kinetic solution very well; for ε=10^-6 and τ=0, the two are consistent.Source: §7.3.3 · p. 26, Figure 10 · p. 28

7.3.4. Model errors

  • Above ε=10^-3, model deviations are O(ε) as truncated moments miss strong non-equilibrium effects; below it, velocity and density errors converge monotonically while pressure is non-monotonic.The ε≈10^-3 transition is specific to these tests; pressure error reaches a minimum at ε=10^-4 before increasing at smaller ε, possibly due to round-off.Source: §7.3.4 · pp. 28–29

8. Conclusions

The paper presents a stable hyperbolic moment-relaxation model and a high-order asymptotic-preserving scheme for low-Mach kinetic equations. Experiments show close agreement with BGK references in continuum-dominated regimes, while identifying limitations in more kinetic regimes and at extremely small scaling parameters.

8. Conclusions

  • Micro-macro decomposition yields a closed moment system whose hyperbolic core simplifies Grad’s linear moment framework and bypasses high-dimensional velocity-space discretization.Source: §8 · p. 29
  • The model is structurally stable in the incompressible limit and consistent with INS; the high-order AP scheme degenerates into an INS projection method as ε → 0.The numerical framework combines IMEX Runge-Kutta schemes with WENO and central-difference spatial reconstructions.Source: §8 · p. 29
  • In continuum-dominated tests, the relaxation model agrees closely with BGK references, while deviations grow for ε > 10^-3 and pressure errors may degrade below ε < 10^-4.The authors suggest round-off effects as a possible explanation for the pressure behavior, but note that further precision tests are needed.Source: §8 · p. 29
  • As ε → 0, the continuous moment approach converges toward the high-fidelity reference while avoiding velocity-space discretization; stiff-regime speedups reach several orders of magnitude.The framework is formulated for arbitrary spatial dimensions, subject to the model-validity qualifications discussed in the paper.Source: §8 · p. 29
  • Future work includes more complex collision operators, higher kinetic moments for rarefied regimes, high-order AP boundary conditions, and uncertainty-quantification studies.Source: §8 · pp. 29–30
Loading 2610.00811v1…