InterPhaseChangeFoam Solver

• InterPhaseChangeFoam is an advanced CFD solver that employs conservative, sharp-interface techniques to simulate compressible multiphase flows with phase change.
• It utilizes a four-wave approximate Riemann solver to rigorously enforce jump conditions and ensure conservation of mass, momentum, and energy.
• The modular finite-volume framework supports dynamic model selection and extensive validation against canonical tests for high-fidelity simulations.

InterPhaseChangeFoam is an advanced computational fluid dynamics (CFD) solver designed for the simulation of compressible multiphase flows with phase change and sharp interface representation. This solver, originating in the OpenFOAM ecosystem, implements fully conservative, sharp-interface methodologies and incorporates phase change via robust Riemann problem solutions at fluid interfaces. Its development is directly influenced by recent research aimed at ensuring strict conservation, resolving interface physics at cell thickness without numerical smearing, and achieving thermodynamic consistency and numerical stability in phase-changing two-phase flows (Long et al., 2021, Scheufler et al., 2021).

1. Governing Equations and Interface Jump Conditions

The foundation of InterPhaseChangeFoam-style solvers is the separate solution of bulk equations for each phase, typically liquid and vapor, under the inviscid Euler equation in conservative form:

$$\mathbf{U} = (\rho, \rho u, \rho v, \rho w, \rho E)^T$$

$$\frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{U}) = 0$$

with fluxes

$$\mathbf{F}(\mathbf{U}) = \begin{bmatrix} \rho \mathbf{u} \ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \ \mathbf{u}(\rho E + p) \end{bmatrix}.$$

Each phase employs its respective equation of state (EOS), such as the stiffened-gas for liquid and ideal-gas or Helmholtz-EOS for vapor phases (Long et al., 2021).

At the interface, sharp jump conditions are imposed consistent with conservation of mass, momentum, and energy. Across an interface moving at speed $S_p$ with normal $\hat{n}$:

$[\rho(V - S_p)] = 0$

$[\rho(V - S_p) V + p] = -\sigma \kappa$

$$[\rho(V - S_p)(e + \frac{1}{2} V^2) + p V] = j Q_\text{lat} - \sigma \kappa S_p$$

where $[a] \equiv a^*_\text{vap} - a^*_\text{liq}$, $V = \mathbf{u} \cdot \hat{n}$, $j$ is the phase-change mass flux, $\sigma$ is the surface tension, $\kappa$ is the interface curvature, and $Q_\text{lat}$ is the latent heat (Long et al., 2021).

2. Four-Wave Approximate Riemann Solver with Phase Change

The interfacial fluxes are determined by solving a generalized Riemann problem incorporating phase change. Rather than tackling the exact five-wave/8D root-finding structure, a four-wave HLL-type Riemann solver is used. It consists of left and right acoustic waves (with Davies’ speed estimates $S_L = V_L - c_L$, $S_R = V_R + c_R$), a single contact/phase wave, and enforces all jump conditions strictly for conservation. The approach reduces the root-finding to a scalar iteration on the phase-change mass flux $j$ (Long et al., 2021).

All star states and wave speeds in the Riemann fan can be expressed analytically as functions of $j$:

$$R_r: \{W_L, W_R, S_L, S_R, \sigma, \kappa, Q_\text{lat}, j\} \rightarrow \{\text{star states}, S_p, S_c\}$$

Final closure is enforced by the phase flux model, typically Hertz–Knudsen: $$f_m = \frac{1}{2\pi R_v} \left[ \lambda_\text{evap} \frac{p_\text{sat}(T_\text{liq}^*)}{\sqrt{T_\text{liq}^*}} - \lambda_\text{cond} \frac{p_\text{vap}^*}{\sqrt{T_\text{vap}^*}} \right]$$ where $T^*$ is inferred from star state pressure/internal energy, ensuring EOS consistency. Steffensen’s method accelerates rootfinding for $j$ with typical convergence in 3–4 iterations. If $\Delta T_\text{contact}>50$ K, stability safeguards override the iterative solution (Long et al., 2021).

A crucial insight is that phase-change fluxes $f_m$ must always use the adjacent Riemann-fan (star) states just next to the interface—not the far-field “left/right” states—to preserve consistency under transient evolution (Long et al., 2021).

3. Finite-Volume Integration and Data Structures

InterPhaseChangeFoam-style algorithms are implemented within a finite-volume framework. The primary fields and their locations are:

Quantity	Location	Comment
$\phi$ (level-set)	Cell-center	Captures interface geometry and motion
$\alpha$ (liquid fraction)	Cell-center	Indicator, reinitialized from level set
$U_\text{liq}$, $U_\text{vap}$	Cell-center	Conservative variables per phase
$F_\text{inv}$	Face	Inviscid fluxes (WENO/HLLC reconstruction)
$X_\text{int}$	Face	Interfacial (“cut-face”) fluxes from Riemann
$\hat{n}$, $\kappa$, $\Delta\Gamma$	Face/cell	Interface geometry: normal, curvature, length
$u_\text{int}$	Cell-center	Interface velocity for level-set advection

Time-stepping is performed using a two-stage Runge–Kutta (RK2) approach. Each RK stage alternates inviscid bulk flux updates (WENO-HLLC schemes) and interfacial (Riemann) flux computation on “cut-faces”, with geometric and level-set fields updated accordingly. The advection of $\phi$ is handled by

$$\frac{\partial\phi}{\partial t} + u_\text{int} \cdot \nabla\phi = 0$$

using high-order schemes. After each time step, $\alpha$ is recomputed as a Heaviside function of $\phi$, and geometric properties are recalculated (Long et al., 2021).

4. Model Selection and Runtime Flexibility

The OpenFOAM-based TwoPhaseFlow framework supports modular model selection at runtime, governed by input dictionaries. All curvature, surface-tension, and phase-change models are implemented as derived classes with static registration, enabling instant switching between, for example, different surface tension force models (height-function, parabolic fit, reconstructed distance function, gradAlpha) and phase change models (Schrage, explicitGrad, implicitGrad) (Scheufler et al., 2021).

This modularity allows researchers to compare alternative physical models or numerical methods on identical test cases by editing configuration files, not recompiling code. Thermophysical properties are likewise selected via OpenFOAM’s mixtureThermo or purePhaseModel facilities (Scheufler et al., 2021).

5. Validation, Numerical Considerations, and Performance

Comprehensive validation is integral to the framework. The sharp-interface solver is benchmarked on canonical problems, including:

• 1D Riemann tests (n-dodecane, water): Agreement with the exact five-wave solution demonstrates star-state accuracy if $j$ uses adjacent states.
• 2D radial droplet evaporation/condensation: Comparison with 1D quasi-analytic radial solvers verifies pressure, density, temperature, and integrated mass transfer.
• Oscillating droplet under surface tension: Simulated oscillation period matches Rayleigh’s analytical formula, validating sharp interface and surface tension models.
• Shock–droplet (or bubble) interaction: Reproduction of vaporization-induced secondary shocks and interface deformation in 2D/3D.
• Conservation diagnostics: Tracking mass in each phase confirms global conservation to round-off (Long et al., 2021).

Benchmarking of the modular (TwoPhaseFlow) framework encompasses classic phase change and capillarity tests, e.g., Stefan problem, “sucking-interface” problem, static circle/sphere reconstruction, viscous sine-wave oscillation, and advected contact-angle circle. Results demonstrate that new curvature and phase change models (implicitGrad, RDF, fitParaboloid, height-function) consistently outperform legacy approaches in accuracy and reduction of spurious currents (Scheufler et al., 2021).

For numerical stability, time steps are limited by a CFL (Courant-Friedrichs-Lewy) condition, typically

$\max(|u| + c) \Delta t / \Delta x \leq 0.6$

where $c$ is the local sound speed. The Steffensen iteration for $j$ is capped at 5 iterations, with fall-back to $j=0$ or previous values if convergence fails or contact temperature jumps exceed thresholds (Long et al., 2021).

6. Implementation and Practical Structure

Implementation closely follows the pseudo-code outline proposed in (Long et al., 2021), transcribed into OpenFOAM C++ structures. Key classes include volScalarField for $\phi$, $\alpha$, $j$; volVectorField for velocity fields; and surfaceScalarField for bulk and interface fluxes. The “fourWaveRiemannSolver” member function encapsulates the core interface physics. The entire solver is structured around three segregated update loops (for $\alpha$, $\mathbf{U}$–$p$, $T$), with all surface tension and phase-change modules dynamically switchable at runtime (Scheufler et al., 2021).

Boundary and initial conditions are specified using conventional OpenFOAM field dictionaries, with phase change and surface tension module selection determined in auxiliary property files. Support for dynamic meshes and adaptive refinement is available for advanced applications (Scheufler et al., 2021).

7. Relation to Broader Two-Phase Flow Methodology

The InterPhaseChangeFoam methodology represents a significant advance over earlier ghost-fluid and smeared-interface approaches by enforcing strict conservation and using sharp interface localization. The coupling of phase change and surface tension at the interface is handled robustly at the finite-volume level, facilitating investigation of complex multiphase phenomena under strong heat and mass transfer. The modular design within the OpenFOAM-based TwoPhaseFlow library fosters rapid prototyping, comparison of models, and extension by the research community (Scheufler et al., 2021).

A plausible implication is that, by enforcing adjacent-state flux evaluation and fully conservative update schemes, the InterPhaseChangeFoam architecture mitigates several previously identified sources of numerical inconsistency, making it suitable for high-fidelity simulation of multiphase flows with sharp interfaces and complex phase change dynamics.