Source-linked AI summary
A Robust All-Mach Six-Equation Diffuse-Interface Method for Multiphase Flows with Surface Tension
Ghanshyam Bharate, J. C. Mandal
TL;DR
Compressible multiphase simulation must remain accurate across low and high Mach numbers while handling interfaces, capillarity, and other physical effects. The paper develops a six-equation finite-volume method with capillary-aware HLLC and pressure relaxation plus a pressure-sensor low-Mach correction; benchmark cases demonstrate accurate interface and low-Mach features, with reduced numerical currents and stable behavior near strong pressure gradients.
Problem
Diffuse-interface simulations face excessive low-Mach numerical dissipation, while existing six-equation formulations do not explicitly preserve capillary pressure jumps through both wave propagation and pressure relaxation.
Source: §1 · p. 2Method
The method combines a two-pressure six-equation model with continuum-surface-force effects, a capillary-aware HLLC solver and pressure relaxation, and a pressure-sensor-modified low-Mach scaling.
Source: §1 · p. 3, §4 · p. 16Results
Benchmark cases demonstrate accurate oscillating-droplet frequency and energy-decay metrics, capture of rising-bubble interface dynamics and low-Mach features, and reduced spurious oscillations near strong pressure gradients.
Source: §5.3 · p. 23, §6 · p. 31Takeaways & Limitations
The framework is validated for multiphase benchmarks involving surface tension, viscosity, gravity, and compressibility across a broad range of flow regimes.
Source: §1 · p. 3, §6 · p. 31The model uses a stiffened-gas equation of state for each phase.
Source: §2 · p. 4
Abstract
from arXiv · showhide
A robust finite-volume framework is presented for the simulation of compressible multiphase flows with surface tension across a wide range of Mach numbers. The method is based on a two-pressure, six-equation diffuse interface model incorporating viscous, gravitational, and capillary effects through the continuum surface force formulation. To consistently account for capillary-induced pressure jumps, an HLLC Riemann solver is developed using generalized Riemann invariant analysis, and the instantaneous pressure relaxation procedure is modified to preserve the Laplace pressure jump during phase equilibration. To overcome the excessive numerical diffusion of conventional approximate Riemann solvers in the low-Mach regime, a robust low-Mach correction is proposed by extending our previous formulation with a modified scaling strategy that remains stable in regions of strong pressure variation. The resulting method retains accuracy from nearly incompressible flows to compressible regimes while preserving the robustness of the six-equation formulation. The numerical framework is validated using a series of benchmark problems involving surface tension, viscosity, gravity, and compressibility. The results demonstrate accurate prediction of interface dynamics, capillary pressure, and low-Mach flow features, while significantly reducing numerical dissipation without compromising stability. The proposed methodology provides an efficient and reliable approach for the simulation of complex multiphase flows spanning a broad range of flow regimes.
1. Introduction
The paper develops a six-equation diffuse-interface framework for multiphase flows with capillarity across a broad range of Mach numbers. It targets low-Mach dissipation and capillary-pressure consistency while retaining robustness.
1. Introduction
- The six-equation model retains pressure non-equilibrium and, with instantaneous pressure relaxation, avoids several numerical difficulties of reduced five-equation formulations.The reduced model can struggle with volume-fraction positivity during strong expansion or compression and can produce non-monotonic mixture sound speed and spurious sonic points.Source: §1 · p. 2
- Approximate Riemann solvers can cause excessive low-Mach dissipation because their artificial dissipation scales with acoustic wave speed rather than fluid velocity.The authors’ earlier correction improves low-Mach accuracy without extra time-step constraints, but its original scaling may become unstable under strong pressure gradients.Source: §1 · pp. 2–3
- Existing six-equation formulations do not explicitly adapt both the Riemann solver and pressure relaxation to preserve capillary pressure jumps.The paper identifies this as a gap in consistently incorporating capillary effects into wave propagation and phase equilibration.Source: §1 · p. 3
- The method combines a two-pressure six-equation diffuse-interface model with viscous, gravitational, and capillary effects, modifying HLLC wave propagation and pressure relaxation to preserve Laplace pressure jumps.A modified low-Mach scaling strategy is designed to remain stable under large pressure variations; benchmarks cover surface tension, viscosity, gravity, and compressibility.Source: §1 · p. 3
2. Multiphase model
The model extends a two-pressure six-equation multiphase system with viscous, capillary, and gravitational effects. It includes phasic pressure relaxation and uses a stiffened-gas equation of state for each phase.
2. Multiphase model
- The model includes mixture momentum and total-energy equations with viscous stress, surface-tension force, and gravity contributions.The total-energy equation is included to ensure mixture total-energy conservation, which is essential in the presence of shock waves.Source: §2 · pp. 4–5
- Phasic internal-energy equations contain pressure-relaxation terms proportional to the pressure difference p1 − p2, using relaxation rate Kp and interface pressure pI.The model uses a stiffened-gas equation of state for each phase.Source: §2 · p. 4
- The viscous formulation assumes mixture and phasic stress tensors are equal, and defines mixture viscosity as the simple average of phasic viscosities.The mixture quantities are distinguished from phase-specific volume fraction, density, pressure, and internal energy.Source: §2 · pp. 4–5
- Capillary forces enter the multiphase system as a volumetric source term through the continuum surface force approach.The passage identifies surface-tension coefficient σ and interface curvature κ.Source: §2 · p. 5
3. Numerical Method
The six-equation model is solved in an evolution step followed by pressure relaxation. The numerical-method overview also identifies curvature estimation and viscous-gradient evaluation in its two figure panels.
3. Numerical Method
- The solution procedure first evolves the system without relaxation terms, then solves an ODE system retaining only pressure relaxation.The two-step procedure is described as the approach used to solve the six-equation model.Source: §3 · p. 5
- Figure 1(a) shows a quadrilateral-cell stencil for curvature estimation, while Figure 1(b) shows Green–Gauss velocity-gradient estimation on a diamond control volume.The velocity gradients are evaluated at face pq for viscous fluxes.Source: Figure 1 · p. 6
3.1. Evolution step
The evolution step discretizes convective, capillary, curvature, and viscous contributions using a finite-volume formulation. Its HLLC solver accounts for capillary jumps, while the curvature and viscous terms use mesh-based gradient procedures.
3.1. Evolution step
- The evolution step omits relaxation terms and organizes the remaining system into convective, viscous, and source contributions for finite-volume discretization.Cell area and face length and normal define the geometric quantities used in the discretization.Source: §3.1 · pp. 6–7
3.1.1. Discretization of convective terms
- The convective finite-volume discretization computes normal conservative fluxes at cell interfaces with a Riemann solver and evaluates non-conservative terms using cell-averaged variables and face-normal velocity.The system is written in a finite-volume form over each cell and its neighboring faces.Source: §3.1 · p. 7, §3.1.1 · p. 8
3.1.2. HLLC Riemann solver with capillary terms
- The HLLC solver for the two-pressure six-equation model incorporates capillary effects using generalized Riemann invariant analysis to derive the pressure jump across the contact wave.The solver constructs intermediate states between outer waves, with a shared middle-wave speed; it uses curvature at cell faces to evaluate the capillary relation.Source: §3.1.2 · pp. 8–9
3.1.3. Surface tension terms
- Surface-tension contributions use a Green–Gauss evaluation of the volume-fraction gradient, with interface quantities taken consistently from the capillary HLLC solution.The total-energy surface-tension term is expressed in divergence form for discretization.Source: §3.1.3 · p. 10
3.1.4. Curvature computation
- Curvature is computed from a smoothed volume-fraction field whose gradients use weighted least squares, then corrected iteratively to reduce errors near pure-phase regions.The correction redistributes curvature from cells near ψ = 0.5 to neighboring cells using weights that peak there.Source: §3.1.4 · pp. 10–12
3.1.5. Viscous terms
- Viscous fluxes include phasic viscous dissipation, and face velocity gradients are evaluated with Green–Gauss on a diamond-shaped control volume.Viscosity, volume fraction, and velocity at faces are obtained by averaging neighboring-cell values.Source: §3.1.5 · pp. 12–13
3.2. Relaxation step
The relaxation step projects the evolved non-equilibrium state toward instantaneous mechanical equilibrium while retaining the surface-tension pressure jump. It also corrects the relaxed pressures and energies to enforce mixture total-energy conservation.
3.2. Relaxation step
- Instantaneous pressure relaxation projects the evolved non-equilibrium state onto mechanical equilibrium, with the surface-tension condition p*_1 − p*_2 = σκ retained.The pressure and volume-fraction constraints are applied in the relaxation procedure.Source: §3.2 · pp. 13–14
- During relaxation, α_jρ_j and the velocity components u and v remain constant, reducing the phasic energy equations to a smaller ODE system.The reduced system is discretized before the relaxed state is computed.Source: §3.2 · p. 13
- The final phasic pressures and energies are corrected using the mixture total-energy equation, then the equation of state and surface-tension constraint determine the final pressures.The pressure difference is constrained by p**_1 − p**_2 = σκ.Source: §3.2 · p. 14
3.3. Second order formulation
The second-order formulation reconstructs primitive variables at cell faces, uses a dedicated interface-sharpening reconstruction for volume fraction, and advances the evolution and relaxation operators with SSPRK.
3.3. Second order formulation
- Cell-face primitive variables are reconstructed from neighboring cell averages using a Taylor expansion truncated after first-order terms.Solution-dependent weighted least squares computes gradients for primitive variables other than volume fraction.Source: §3.3 · pp. 14–15
- The volume fraction uses the Overbee limiter with distance-weighted least-squares gradients, a dedicated reconstruction that reduces interface smearing.This procedure is specific to α1.Source: §3.3 · p. 15
- Second-order temporal accuracy is obtained with SSPRK applied to the evolution and relaxation operators.The scheme is written as the combined implementation of those two operators.Source: §3.3 · p. 15
4. Low Mach correction with modified scaling factor
The modified low-Mach correction uses a pressure-sensitive scaling factor to reduce numerical diffusion while addressing instabilities near pressure jumps. Its pressure function accounts for jumps arising from shocks, surface tension, or hydrostatic gradients.
4. Low Mach correction with modified scaling factor
- The modified scaling factor is introduced because the prior low-Mach correction can become unstable near pressure jumps despite reducing excessive diffusion as M → 0.The new factor is designed to reduce the correction according to pressure-jump strength.Source: §4 · p. 16
- The multidimensional pressure function detects pressure jumps and reduces the low-Mach correction according to their strength.The face value is evaluated from local pressure functions across the faces of both adjacent cells.Source: §4 · p. 16
- The modified scaling enhances stability for low-Mach multiphase flows with pressure jumps from surface tension or hydrostatic gradients.The authors report this stability benefit in numerical experiments.Source: §4 · p. 16
5. Results
The benchmarks test low-Mach flows with surface tension, viscosity, and gravity, alongside a high-speed shock-induced bubble collapse. Across these cases, the modified low-Mach scaling improves stability and accuracy, while the method captures capillary pressure jumps and dynamic interface behavior.
5.1. Dam-break
- At t = 4.43, the modified scaling reduces dam-break front-position error from 8.22% to 2.85% against experiment.The simple scaling produces spurious pressure oscillations near the interface; the modified scaling reduces its effect where pressure variation is high.Source: §5.1 · p. 17, Figure 4 · p. 19
5.2. Static droplet
- The static-droplet pressure jump converges toward the Young–Laplace solution with almost second-order accuracy, while parasitic velocities remain below 10^-2 m/s.At Laplace number 10^11, the simple scaling produces unbounded currents and interface distortion; the modified scaling lowers parasitic currents across tested Laplace numbers.Source: §5.2 · pp. 19–21, Figure 5 · p. 20, Figure 6 · p. 21
5.3. Oscillating droplet
- The modified low-Mach correction removes spurious interface oscillations under the large pressure jump, and the oscillating-droplet frequency error is 4.38% on a 192 × 192 mesh.Numerical dissipation decreases with grid refinement and is lower with the correction on the finest grid; the reported comparison work had 66% error on a comparable 200 × 200 mesh.Source: §5.3 · pp. 22–23, Figure 7 · p. 22, Table 1 · p. 24
5.4. Rising bubble
- With the correction, rising-bubble trajectories and shape are in close agreement with the incompressible reference; results change little from Grid B to Grid C.The bubble reaches terminal velocity as viscous drag balances buoyancy, and deforms from circular to ellipsoidal during ascent.Source: §5.4 · p. 25, Figure 11 · p. 26, Figure 12 · p. 27
5.5. Shock-induced bubble collapse
- The shock-induced collapse results agree qualitatively with previous studies and reproduce the reference bubble height and width before collapse.The calculation resolves reflected rarefaction, transmitted shocks, bubble breakup, and a blast wave despite the fluids' strongly differing compressibilities.Source: §5.5 · p. 29, Figure 14 · p. 30, Figure 15 · p. 30
6. Conclusion
The framework combines capillary-aware six-equation flow modeling with a pressure-sensor low-Mach correction. Benchmark results support its accuracy and robustness across the tested multiphase regimes.
6. Conclusion
- The method couples continuum-surface-force capillarity with a two-pressure six-equation model and modifies both the HLLC solver and pressure relaxation to preserve Laplace pressure jumps.The capillary pressure jump across the contact wave is derived using generalized Riemann invariant analysis.Source: §6 · p. 31
- A multidimensional pressure-based sensor modifies the low-Mach scaling to suppress spurious oscillations near strong pressure gradients and produce smooth interface profiles.The correction extends the authors’ earlier approach while addressing its instability in regions of strong pressure variation.Source: §6 · p. 31
- Benchmarks show accurate droplet oscillation metrics, rising-bubble interface dynamics, stable static-droplet currents, improved dam-break agreement, and robust shock-induced bubble collapse without negative pressures.The modified scaling factor is identified as important for dampening parasitic currents and maintaining stability, particularly at lower physical viscosity.Source: §6 · p. 31
Appendix A. Generalised Riemann Invariants Analysis
The appendix derives generalized Riemann invariants for the six-equation system with capillary effects. These relations inform the HLLC solver’s treatment of intermediate waves and the Laplace pressure jump.
Appendix A. Generalised Riemann Invariants Analysis
- The capillary term leaves the multiphase system’s eigenvalues unchanged, with the acoustic waves surrounding six coincident intermediate eigenvalues.The appendix writes the system in quasi-linear form and identifies the mixture sound speed in the eigenvalue expressions.Source: Appendix A · pp. 32–33
- In the six-equation compressible multiphase framework, λ2 through λ6 are represented by a single contact speed S∗ in the HLLC solver; discontinuity invariants are found by GRI analysis of right eigenvectors with eigenvalues λ2,…,6 = u.Source: Appendix A · p. 34
- Across the intermediate wave, the mixture pressure follows the Laplace law while phasic densities and velocity components remain constant.The relations give dp = σκdα1; the analysis also notes that variables other than normal velocity may change across the wave.Source: Appendix A · p. 34
Appendix B. Preventing Computational Breakdown in Strong Shock Waves via Isentropic Pressure Relations
The appendix addresses negative pressures during strong-shock evolution by using an isentropic pressure relation. It also describes an implicit interface-pressure approximation as a further way to reduce their influence during relaxation.
Appendix B. Preventing Computational Breakdown in Strong Shock Waves via Isentropic Pressure Relations
- Using the isentropic relation from Eq. (20), rather than the expression in Eq. (19), prevents large negative pressures during the evolution step.The concern is that phasic pressures from the evolution step can affect volume fraction and final relaxed pressure, potentially causing negative volume fractions or pressures after relaxation.Source: Appendix B · pp. 34–35
- At the stated condition jK ρjK = 1, both expressions converge to the same value; the passage additionally refers to Figure B.16, identified in the cited caption as the pressure ratio.Source: Appendix B · p. 34, Figure B.16 · p. 35
- During relaxation, an implicit approximation of interface pressure using p∗I can further minimize the influence of large negative pressures from the evolution step.The appendix presents this as an additional measure for problems involving high-strength shock waves.Source: Appendix B · p. 35