Multi-Phase Flow

Overview

The FluidDynamics module supports multi-phase and multi-component flow via pseudo-potential methods. These methods introduce a short-range non-ideal force between fluid nodes that drives phase separation into liquid and vapor phases without tracking an explicit interface.

Two pseudo-potential models are available:

Model	Controlled by	EOS
Benzi	Do_Benzi = true	Shan–Chen pseudo-potential
Kupershtokh	Do_Kupershtokh = true	Van der Waals (or custom EOS)

Both models are implemented inside FlowSolverLBM and are configured via the input file.

---

Benzi pseudo-potential model

Physics

The Benzi (Shan–Chen) model introduces a cohesive force between fluid nodes of the same component via a scalar pseudo-potential:

$$\begin{array}{}\text{(1)}& \psi (\rho )=1-e^{-\rho /{\rho }_0}\end{array}$$

The resulting interaction force on a node at $\mathbf{x}$ is:

$$\begin{array}{}\text{(2)}& {\mathbf{F}}_{\text{Benzi}}(\mathbf{x})=-G_b\psi (\rho (\mathbf{x}))\sum _{xyz}{w}_{xyz}\psi (\rho (\mathbf{x}+{\mathbf{e}}_{xyz})){\mathbf{e}}_{xyz}\end{array}$$

where:

• $G_b$ — Benzi interaction strength parameter (Gb matrix, per component pair)
• ${\rho }_0$ — reference density controlling the shape of $\psi$ (rho_0 per component)
• $w_{xyz}$ — D3Q27 lattice weights
• ${\mathbf{e}}_{xyz}$ — lattice velocity vectors

The force ${\mathbf{F}}_{\text{Benzi}}$ drives like-component molecules together and repels different components, producing phase separation.

BenziGas helper

BenziGas provides a lookup table of pre-computed equilibrium state values as a function of the Benzi parameter $G_b$:

struct BenziGas
{
    struct EquilibriumValues_t {
        double BenziParameter;
        double lbVaporDensity;
        double lbLiquidDensity;
        double lbSurfaceTension;
    };

    static EquilibriumValues_t EquilibriumValues(double BenziParameter);
};

The table covers $G_b\in [-4.02,-5.86]$. Values outside this range cause the simulation to abort. Linear interpolation is used between table entries.

Usage: call BenziGas::EquilibriumValues(Gb) before setting up the FlowSolverLBM to determine the correct LiquidDensity, VaporDensity, and SurfaceTension input values.

---

Kupershtokh two-phase model

Physics

The Kupershtokh model is based on the Equation of State (EOS) approach. A potential $\mathrm{\Phi }$ derived from the EOS determines a non-ideal pressure correction:

$$\begin{array}{}\text{(3)}& \mathrm{\Phi }(\rho )=\sqrt{-(G_K\cdot p_{\text{EOS}}(\rho )-\frac{\rho }{3})}\end{array}$$

The interaction force is:

$$\begin{array}{}\text{(4)}& {\mathbf{F}}_{\text{Kuper}}(\mathbf{x})=-2\mathrm{\Phi }(\rho (\mathbf{x}))\sum _{xyz}{w}_{xyz}\mathrm{\Phi }(\rho (\mathbf{x}+{\mathbf{e}}_{xyz})){\mathbf{e}}_{xyz}\end{array}$$

where $G_K$ is the Kupershtokh parameter (lbGK matrix) and $p_{\text{EOS}}$ is the reduced pressure from the Van der Waals (or a user-specified) EOS.

Optimal parameter

The static helper OptimalParaKuper computes the optimal Kupershtokh parameter for a given reduced temperature:

static double OptimalParaKuper(double ReducedTemperature,
                                double GasParameter = 0.01);

This is essential for numerical stability: choosing a sub-optimal $G_K$ causes parasitic currents or instability near phase boundaries.

VanDerWaalsGas helper

VanDerWaalsGas provides the Van der Waals EOS and a lookup table of equilibrium coexistence values:

struct VanDerWaalsGas
{
    // Reduced Van der Waals equation of state
    static double ReducedPressure(double Density, double Temperature,
                                  double CriticalDensity = 1.0,
                                  double CriticalTemperature = 1.0);

    struct EquilibriumValues_t {
        double Temperature;
        double VaporDensity;
        double LiquidDensity;
    };

    static EquilibriumValues_t EquilibriumValues(
        double Temperature,
        double CriticalTemperature = 1.0,
        double LaplacePressure = 0.0);
};

The Van der Waals EOS in reduced form:

$$\begin{array}{}\text{(5)}& p^{\ast }=\frac{8{\rho }^{\ast }{T}^{\ast }}{3-{\rho }^{\ast }}-3({\rho }^{\ast }{)}^2\end{array}$$

where starred quantities are normalized by their critical values. The equilibrium table covers $T^{\ast }\in [0.01,1.00]$ (all in reduced units).

---

Phase density profile

For initialization, a smooth hyperbolic tangent profile is provided:

double DensityProfile(double x, size_t n = 0) const;

This generates a liquid–vapor interface profile for component $n$ at position $x$.

---

Wetting at solid surfaces

Wetting between fluid components and solid phases is controlled via the Wetting parameter matrix. A non-zero wetting parameter modifies the effective density used in the force calculation at solid nodes, implementing contact-angle control.

The DensityWetting field stores this modified density:

Storage3D<double, 1> DensityWetting; ///< ρ / wetting parameter per component

---

InteractionFluidFluid

InteractionFluidFluid is a lightweight companion class for future extensions of fluid–fluid interaction models beyond what is built into FlowSolverLBM.

class InteractionFluidFluid : public OPObject
{
public:
    void Initialize(Settings& locSettings) override;
    void ReadInput(std::string InputFileName) override;
};

Currently this class provides the initialization/input interface; the core interaction forces are computed directly in FlowSolverLBM::CalculateForceTwoPhase().

---

Fluid redistribution

In two-phase flow with moving solids, small fluid parcels can become trapped inside obstacles. FluidRedistributionRange controls the distance over which fluid is redistributed when this occurs, preventing mass loss.

---

Numerical considerations

• Surface tension and interface width are coupled: sharper interfaces require larger $|G_b|$ or $G_K$ but increase numerical artifacts (parasitic currents).
• The Benzi model is stable for $G_b\in [-4.02,-5.86]$ (enforced by BenziGas).
• The Kupershtokh model requires $T^{\ast }<1$ (below critical temperature) for phase separation; at $T^{\ast }\to 1$ the phases merge.
• EnforceMassConservation() should be called every timestep in two-phase simulations to prevent drift in total fluid mass.
• Do_FixPopulations = true is recommended for high density ratio ($>10$) simulations.