Theory

Overview

MultiFlow is a finite-volume numerical framework for the simulation of multiphase flows in arbitrary domains. It solves various formulations of the incompressible and compressible Navier-Stokes equations on unstructured meshes, using a fully coupled pressure-based algorithm. The framework has been under continuous development since 2000 and is now in its third generation.

Governing Equations

The framework solves the conservation laws for mass, momentum and energy:

\[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0\]

\[\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \mathbf{S}\]

\[\frac{\partial \rho h}{\partial t} + \nabla \cdot (\rho \mathbf{u} h) = \frac{\partial p}{\partial t} + \nabla \cdot (\boldsymbol{\tau} \cdot \mathbf{u}) - \nabla \cdot \mathbf{q}\]

where \(\rho\) is the density, \(\mathbf{u}\) the velocity, \(p\) the pressure, \(\boldsymbol{\tau}\) the viscous stress tensor, \(h\) the specific total enthalpy, and \(\mathbf{q} = -k \nabla T\) the heat flux via Fourier's law.

Finite-Volume Discretisation

The governing equations are discretised using a collocated finite-volume method on unstructured meshes, applicable to structured and unstructured grids with cells of arbitrary topology. The discretisation is consistently second-order accurate in space, employing:

Gradient evaluation via the divergence theorem
Advection terms discretised with TVD schemes (Minmod, van Leer) for robustness in the presence of shocks and sharp gradients
Diffusion terms with deferred correction for non-orthogonal meshes
Transient terms discretised with first- or second-order backward Euler schemes

Fully Coupled Pressure-Based Algorithm

A key feature of MultiFlow is that the discretised governing equations — continuity, momentum and energy — are solved simultaneously as a single linear system for pressure, velocity and temperature at each time step:

\[\begin{pmatrix} A_{\rho,p} & A_{\rho,\mathbf{u}} \\ A_{\rho\mathbf{u},p} & A_{\rho\mathbf{u},\mathbf{u}} \\ A_{\rho h,p} & A_{\rho h,\mathbf{u}} & A_{\rho h,T} \end{pmatrix} \begin{pmatrix} \psi_p \\ \psi_\mathbf{u} \\ \psi_T \end{pmatrix} = \begin{pmatrix} \sigma_\rho \\ \sigma_{\rho\mathbf{u}} \\ \sigma_{\rho h} \end{pmatrix}\]

This fully coupled approach — in contrast to segregated methods such as SIMPLE or PISO — provides superior robustness and convergence, particularly for flows with large density gradients, strong source terms, and multiphase interfaces. The linear system is solved using the BiCGSTAB solver and Block-Jacobi preconditioner from the PETSc library.

Momentum-Weighted Interpolation

To avoid pressure-velocity decoupling on collocated meshes, MultiFlow employs a momentum-weighted interpolation (MWI) to determine the advecting velocity at cell faces. The face velocity is given by:

\[\vartheta_f = u_{i,f} n_{i,f} - \hat{d}_f \left( \frac{p_Q - p_P}{s_f} - \rho^*_f \left[ \frac{1-l_{Pf}}{\rho_P} \frac{\partial p}{\partial x_i}\bigg|_P + \frac{l_{Pf}}{\rho_Q} \frac{\partial p}{\partial x_i}\bigg|_Q \right] s_{i,f} \right)\]

This formulation provides robust pressure-velocity coupling across all flow regimes, from the incompressible limit to supersonic flows.

All-Speed Capability

MultiFlow is capable of simulating flows across the full Mach number range — from creeping incompressible flows to supersonic and hypersonic regimes — within a single, unified framework. This is achieved through:

A unified thermodynamic closure combining incompressible fluid models with the Noble-Abel stiffened-gas (NASG) equation of state for ideal and real gases
Newton linearisation of the continuity equation, which acts as a transport equation for density in compressible flows and as a divergence-free constraint in the incompressible limit
A conservative formulation of the energy equation ensuring accurate prediction across strong shocks

The framework has been validated for Mach numbers ranging from \(M \to 0\) to \(M = 239\), including acoustic wave propagation, shock tubes, and supersonic flows with complex shock structures.

Multiphase Flow Modelling

Beyond single-phase flows, MultiFlow supports a range of multiphase flow configurations:

Interfacial flows are modelled using a Volume-of-Fluid (VOF) approach, with a balanced-force formulation ensuring accurate treatment of surface tension. The curvature of the interface is evaluated using a least-squares method from volume fractions, and the momentum-weighted interpolation is extended to handle the large density gradients and source terms present at fluid interfaces.

Dispersed flows are modelled using an Eulerian-Lagrangian approach, in which the continuous phase is solved on the Eulerian mesh while individual particles or droplets are tracked as Lagrangian entities. The coupling between the phases accounts for drag, volume fraction effects, and interphase momentum transfer.

The third generation of MultiFlow employs a combination of DMDA and DMPlex with the DMForest/p4est framework from PETSc, enabling adaptive octree mesh refinement. This allows computational resources to be concentrated in regions of interest — such as near interfaces, shocks, or particles — while maintaining efficiency across the rest of the domain.

Parallelisation

MultiFlow is designed for high-performance computing. Domain decomposition and inter-processor communication are handled through the PETSc library, enabling efficient scaling across large numbers of processors. The framework has been applied to simulations with tens of millions of cells on distributed-memory HPC clusters.

Key References

[1] B. van Wachem, V.R. Gopala, A coupled solver approach for multiphase flow calculations on collocated grids, in: European Conference on Computational Fluid Dynamics, ECCOMAS CFD, TU Delft, 2006, pp. 1–16.

[2] B. van Wachem, A. Benavides, V.R. Gopala, A coupled solver approach for multiphase flow problems, in: 6th International Conference on Multiphase Flows 2007, Leipzig, Germany, 2007, Paper No 183.

[3] F. Denner, B. van Wachem, Fully-coupled balanced-force VOF framework for arbitrary meshes with least-squares curvature evaluation from volume fractions, Numerical Heat Transfer Part B: Fundamentals 65 (2014) 218–255. https://doi.org/10.1080/10407790.2013.849996

[4] C.-N. Xiao, F. Denner, B. van Wachem, Fully-coupled pressure-based finite-volume framework for the simulation of fluid flows at all speeds in complex geometries, Journal of Computational Physics 346 (2017) 91–130. https://doi.org/10.1016/j.jcp.2017.06.009

[5] P. Bartholomew, F. Denner, M.H. Abdol-Azis, A. Marquis, B. van Wachem, Unified formulation of the momentum-weighted interpolation for collocated variable arrangements, Journal of Computational Physics 375 (2018) 177–208. https://doi.org/10.1016/j.jcp.2018.08.030

[6] M.H. Abdol Azis, F. Evrard, B. van Wachem, An immersed boundary method for incompressible flows in complex domains, Journal of Computational Physics 378 (2019) 770–795. https://doi.org/10.1016/j.jcp.2018.10.048

[7] F. Denner, F. Evrard, B. van Wachem, Breaching the capillary time-step constraint using a coupled VOF method with implicit surface tension, Journal of Computational Physics 459 (2022) 111128. https://doi.org/10.1016/j.jcp.2022.111128

[8] R. Janodet, B. van Wachem, F. Denner, A fully-coupled algorithm with implicit surface tension treatment for interfacial flows with large density ratios, Journal of Computational Physics 520 (2025) 113520. https://doi.org/10.1016/j.jcp.2024.113520