The Hardest Equations in CFD — CMPS Strategies for Pressureless Euler (Dilute Eulerian Particles)

1. Introduction

Among all systems of equations in Computational Fluid Dynamics, the pressureless Euler equations stand out as some of the most difficult to solve. They describe the motion of a dilute dispersed phase—typically solid or liquid particles—treated as a continuum but without an internal pressure term. The absence of pressure removes the natural wave propagation mechanism that exists in gas dynamics, leaving the system weakly hyperbolic or even non-hyperbolic.

This lack of hyperbolicity causes serious numerical challenges. The equations admit no real sound speed, so traditional Riemann solvers cannot define proper wave speeds. Even small perturbations may lead to non-physical solutions such as δ-shocks, vacuum regions, or undefined contact discontinuities. These behaviors make the system extremely sensitive to discretization errors and often unstable under standard CFD schemes.

In the dilute limit, where the particle volume fraction and density approach zero, the situation becomes even more challenging. Conservation must still be maintained, yet the convective fluxes dominate completely and may collapse into singular structures. This is why the pressureless Euler equations are often regarded as the hardest equations in CFD.

CMPS addresses these issues with a fully coupled, fully implicit framework that combines a controlled artificial compressibility, regularized fluxes, and exact Jacobians, transforming an almost ill-posed system into a numerically and physically consistent model for dilute particle transport.

2. The Mathematical Difficulty of the Pressureless Euler Equations

The core pressureless Euler system in one spatial dimension (without the volume-fraction term) is

$$ \frac{\partial \rho_p}{\partial t} + \frac{\partial (\rho_p u_p)}{\partial x} = 0, \qquad \frac{\partial (\rho_p u_p)}{\partial t} + \frac{\partial (\rho_p u_p^2)}{\partial x} = 0. $$

The flux Jacobian of this system is

$$ A = \frac{\partial \mathbf{F}}{\partial \mathbf{U}} = \begin{bmatrix} 0 & 1 \\ -u_p^2 & 2u_p \end{bmatrix}, \qquad \text{with eigenvalues } \; \lambda_1 = \lambda_2 = u_p. $$

Both eigenvalues are identical and the matrix has only one independent eigenvector, so it cannot be diagonalized. The system is therefore not strictly hyperbolic. All characteristics collapse to the particle velocity \(u_p\), eliminating wave separation and breaking the foundation of classical upwind CFD methods. The result is well-known pathologies: δ-shocks, spontaneous vacuum formation, and extreme sensitivity to discretization.

3. Restoring Hyperbolicity in CMPS

CMPS restores a workable wave structure by introducing a small, controlled artificial pressure for the particle phase,

$$ p_p = \varepsilon_p \; \rho_p^{\gamma_p}, $$

which yields a finite particle-wave speed

$$ c_d = \sqrt{\varepsilon_p \, \gamma_p \, \rho_p^{\gamma_p - 1}}. $$

This converts the degenerate spectrum into two distinct real eigenvalues,

$$ \lambda_{1,2} = u_p \pm c_d, $$

re-establishing numerical hyperbolicity and enabling stable upwind fluxes. The compressibility is localized and bounded: a frictional threshold in local concentration ensures that in the truly dilute regime the artificial pressure vanishes, while a Mach-based limiter caps \(\varepsilon_p\) to avoid unphysical stiffness in slow flows.

3.1 Computing \(\varepsilon_p\): Activation and Limiting

The artificial pressure must be inactive in very dilute regions and bounded elsewhere. CMPS accomplishes this with two controls:

• Activation by concentration: Let \(\alpha = \rho_p/\rho_{p,mat}\) be the local particle volume fraction (\(\rho_{p,mat}\) is material density). Define a frictional threshold \(\alpha_{fr}\) and a maximum packing \(\alpha_{max}\). The face-average \(\alpha_{face} = (\alpha_L + \alpha_R)/2\) is clipped to \([0,\alpha_{max}]\). If \(\alpha_{face} < \alpha_{fr}\), then set \(\varepsilon_p^{eff} = 0\) and the particle phase remains pressureless.
• Mach-based limiter: To keep the pseudo-sound speed bounded by the local particle speed, enforce $$ c_d \le \frac{|\mathbf{u}_p|}{M_{min}}. $$ Using \(c_d^2 = \varepsilon_p \, \gamma_p \, \rho_p^{\gamma_p-1}\), the effective coefficient becomes $$ \varepsilon_p^{eff} = \min\!\left( \varepsilon_p,\; \frac{|\mathbf{u}_p|^2}{\gamma_p\, \rho_p^{\gamma_p-1}\, M_{min}^2} \right). $$ This guarantees that the artificial wave speed never dominates the convective scale.

These rules are applied separately to left/right states at each face, producing \(p_{p,L/R} = \varepsilon_{p,L/R}^{eff} \rho_{p,L/R}^{\gamma_p}\) and \(c_{d,L/R} = \sqrt{\varepsilon_{p,L/R}^{eff} \gamma_p \rho_{p,L/R}^{\gamma_p-1}}\). In truly dilute zones the method reverts to the pressureless limit by design.

3.2 Regularized Flux Evaluation

Using the artificial wave speeds, the face flux is computed with a robust Rusanov/HLL-type formula,

$$ \mathbf{F}_{\text{face}} = \tfrac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \tfrac{1}{2}\, \lambda_{\max} (\mathbf{U}_R - \mathbf{U}_L), \quad \lambda_{\max} = \max\big( |u_L| + c_{d,L},\, |u_R| + c_{d,R} \big). $$

This adapts dissipation to local gradients and smoothly reverts to the pressureless limit where the artificial pressure is inactive. In smooth regions the scheme remains low-diffusive; near sharp gradients it provides the dissipation needed to prevent δ-shocks and negative densities.

3.3 Coupling and Conservation

The particle equations remain consistently coupled to the gas phase via interphase momentum and heat exchange with equal-and-opposite source terms, ensuring strict global conservation. The artificial pressure affects only the particle-phase flux evaluation and does not disturb the conservation balance between phases.

3.4 Implicit Stabilization and Exact Jacobians

All fluxes and sources are linearized and solved in a fully coupled implicit system, using exact Jacobians obtained via automatic differentiation. This suppresses short-wavelength instabilities and preserves Newton convergence even under strong interphase drag and extreme density ratios.