---
title: InterPhaseChangeFoam Solver
url: https://www.emergentmind.com/topics/interphasechangefoam-solver
type: topic
---

# InterPhaseChangeFoam Solver

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 [2110.07995][2103.00870].

## 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 [2110.07995].

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 [2110.07995].

## 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$ [2110.07995].

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 [2110.07995].

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 [2110.07995].

## 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 [2110.07995].

## 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) [2103.00870].

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 [2103.00870].

## 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 [2110.07995].

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 [2103.00870].

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 [2110.07995].

## 6. Implementation and Practical Structure

Implementation closely follows the pseudo-code outline proposed in [2110.07995], 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 [2103.00870].

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 [2103.00870].

## 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 [2103.00870].

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.

Source: https://www.emergentmind.com/topics/interphasechangefoam-solver