Consistent multiple-relaxation-time lattice Boltzmann method for the volume averaged Navier-Stokes equations

Recently, we notice that a pressure-based lattice Boltzmann (LB) method was established to recover the volume-averaged Navier-Stokes equations (VANSE), which serve as the cornerstone of various fluid-solid multiphase models. It decouples the pressure from density and exhibits excellent numerical performance, however, the widely adopted density-based LB scheme still suffers from significant spurious velocities and inconsistency with VANSE. To remedy this issue, a multiple-relaxation-time LB method is devised in this work, which incorporates a provisional equation of state in an adjusted density equilibrium distribution to decouple the void fraction from density. The Galilean invariance of the recovered VANSE is guaranteed by introducing a penalty source term in moment space, effectively eliminating unwanted numerical errors. Through the Chapman-Enskog analysis and detailed numerical validations, this novel method is proved to be capable of recovering VANSE with second-order accuracy consistently, and well-suited for handling void fraction fields with large gradients and spatiotemporal distributions.

💡 Research Summary

The paper addresses a long‑standing inconsistency in lattice Boltzmann (LB) models that aim to solve the volume‑averaged Navier‑Stokes equations (VANSE), which are the governing equations for coarse‑grained simulations of fluid‑solid multiphase systems. Existing density‑based LB schemes, typically using a single‑relaxation‑time (BGK) collision operator, suffer from severe spurious velocities and pressure oscillations when the void fraction (porosity) varies sharply. The root cause is the coupling of pressure to density in the equilibrium distribution and the need to discretise the void‑fraction gradient in the forcing term, which introduces high‑order numerical errors (see Eq. 12 in the manuscript).

To overcome these drawbacks, the authors propose a fully consistent multiple‑relaxation‑time (MRT) LB method for VANSE (denoted MRT‑LB‑VANSE). The key innovations are:

1. Provisional Equation of State (EOS) – An artificial EOS (p = \kappa \rho) is introduced, where (\kappa) is a constant chosen close to the minimum void fraction in the domain. This decouples pressure from the void fraction, allowing the pressure‑correction force to be written as (\mathbf{F}_p = -\nabla(\kappa\rho) + \nabla(p)). Because (\kappa) is independent of (\phi), the force term no longer contains the problematic (\phi\nabla p) discretisation, eliminating the dominant source of spurious velocities.

2. Adjusted Equilibrium Distribution – The equilibrium density distribution is constructed from a third‑order Hermite expansion of the Maxwellian (Eq. 14). Its zeroth and first moments recover the correct mass and momentum ((\phi\rho) and (\phi\rho\mathbf{u})), while the second‑order moment incorporates the provisional EOS, fixing the pressure term. The third‑order moment, however, introduces a Galilean‑invariance error in the viscous stress tensor.

3. MRT Collision Operator with Penalty Source – Following the approach of Li et al., a source term (\mathbf{C}) is added in moment space (Eq. 20). This term precisely cancels the unwanted contributions from the third‑order moment, restoring Galilean invariance and ensuring that the recovered viscous stress matches the continuous form. The relaxation matrix (\Gamma) contains distinct relaxation rates for the energy, energy‑square, and heat‑flux modes (e.g., (s_e = 1.2), (s_\epsilon = 1.4)), providing additional numerical stability compared with the BGK model.

A Chapman‑Enskog expansion (Appendix A) demonstrates that the MRT‑LB scheme exactly recovers the VANSE continuity equation (\partial_t(\phi\rho) + \nabla!\cdot(\phi\rho\mathbf{u}) = 0) and the momentum equation (\partial_t(\phi\rho\mathbf{u}) + \nabla!\cdot(\phi\rho\mathbf{u}\mathbf{u}) = -\phi\nabla p + \nabla!\cdot(\phi\boldsymbol{\tau}) + \phi\mathbf{F}). The pressure term appears as (-\phi\nabla p) without any extra spurious contribution, guaranteeing consistency even for highly non‑uniform (\phi).

The authors validate the method through three sets of numerical experiments:

• Method of Manufactured Solutions (MMS) – Spatially and temporally varying void‑fraction fields are prescribed, and the MRT‑LB solution exhibits second‑order convergence in both space and time, confirming the theoretical accuracy.
• Spurious‑Velocity Test – A porous channel with a steep void‑fraction gradient is simulated. The MRT‑LB model reduces the maximum spurious velocity by more than an order of magnitude relative to the BGK‑based density scheme and the previously proposed pressure‑based SR‑TB model.
• Viscosity and Reynolds‑Number Sweep – Simulations across a wide range of kinematic viscosities (0.01 – 10 lattice units) and Reynolds numbers demonstrate stable operation, including cases where the void fraction drops below 0.5—a regime where earlier approaches diverge.

Overall, the paper delivers a robust, second‑order accurate LB framework that faithfully reproduces VANSE for complex, time‑dependent void‑fraction distributions. By decoupling pressure from density, employing an MRT collision operator, and correcting viscous stresses via a moment‑space source term, the method overcomes the principal limitations of prior LB schemes. This advancement opens the door to reliable LB simulations of fluid‑solid multiphase flows, porous media, and coupled particle‑fluid dynamics, and provides a solid foundation for extending the approach to non‑Newtonian fluids, heat transfer, and reactive transport problems.