A variational front-tracking method for multiphase flow with triple junctions

We present and analyze a variational front-tracking method for a sharp-interface model of multiphase flow. The fluid interfaces between different phases are represented by curve networks in two space dimensions (2d) or surface clusters in three space dimensions (3d) with triple junctions where three interfaces meet, and boundary points/lines where an interface meets a fixed planar boundary. The model is described by the incompressible Navier–Stokes equations in the bulk domains, with classical interface conditions on the fluid interfaces, and appropriate boundary conditions at the triple junctions and boundary points/lines. We propose a weak formulation for the model, which combines a parametric formulation for the evolving interfaces and an Eulerian formulation for the bulk equations. We employ an unfitted discretization of the coupled formulation to obtain a fully discrete finite element method, where the existence and uniqueness of solutions can be shown under weak assumptions. The constructed method admits an unconditional stability result in terms of the discrete energy. Furthermore, we adapt the introduced method so that an exact volume preservation for each phase can be achieved for the discrete solutions. Numerical examples for three-phase flow and four-phase flow are presented to show the robustness and accuracy of the introduced methods.

💡 Research Summary

This paper introduces a variational front‑tracking framework for sharp‑interface multiphase flow problems that involve triple junctions, i.e., points (in 2‑D) or lines (in 3‑D) where three fluid interfaces meet. The authors consider a bounded polygonal (or polyhedral) domain Ω that is partitioned into a finite number of bulk phases Rℓ, each occupied by an incompressible Newtonian fluid with its own density ρℓ and viscosity ηℓ. In the bulk, the standard incompressible Navier–Stokes equations govern the velocity field u and pressure p. On each interface Γi the usual continuity of velocity, the jump condition for the normal stress balanced by surface tension γi times the mean curvature κi, and the kinematic condition V·νi = u·νi are imposed. At a triple junction Tk the authors enforce the classical force‑balance condition Σj=1^3 γsj μsj = 0, where μsj are the conormals of the three meeting interfaces, together with an orientation vector ok that fixes the sign convention. Boundary points/lines where an interface meets a fixed planar wall are treated with a 90° contact‑angle condition and a no‑penetration condition.

The core contribution is a weak (variational) formulation that couples a parametric description of the moving interfaces with an Eulerian description of the bulk flow. The interface motion is represented by a velocity vector V derived from the time derivative of the parametric maps xi(·,t). By testing the Navier–Stokes equations with arbitrary velocity fields and the interface evolution equations with arbitrary test functions, the authors obtain a single coupled variational statement (equation (3.5) in the paper). This formulation is designed so that the total free energy – kinetic energy plus interfacial energy – dissipates exactly at the discrete level, mirroring the continuous energy law (2.10).

For spatial discretisation the authors adopt an unfitted finite‑element approach: a fixed background mesh (typically piecewise linear) is used for the bulk variables, while the interfaces are discretised by independent parametric meshes (e.g., piecewise quadratic curves in 2‑D or surface patches in 3‑D). The coupling between the two meshes is handled by cut‑FEM / XFEM‑type techniques that allow integration over intersected elements without remeshing. Time integration is performed with a θ‑scheme; choosing θ=1 yields a fully implicit, linear system at each time step. The resulting fully discrete scheme is linear, unconditionally stable (Theorem 4.2), and admits a unique solution under mild assumptions.

A notable challenge in multiphase flow is the preservation of the volume of each phase. The authors adapt a technique from previous work on two‑phase flow: they compute a time‑averaged interface normal (\hat n^{n+1/2}) and modify the discrete kinematic condition so that the discrete volume change of each phase vanishes exactly. This modification introduces a mild nonlinearity, but the authors show that the scheme remains solvable and retains the unconditional energy stability.

The paper provides a thorough theoretical analysis. The authors prove an energy dissipation inequality for the discrete solution, demonstrate that the discrete curvature balance at triple junctions mimics the continuous analogue, and establish existence and uniqueness of the discrete solution. They also discuss mesh quality: the tangential freedom inherent in the BGN (Barrett‑Garcke‑Nürnberg) formulation leads to asymptotically equidistributed nodes in 2‑D and well‑conditioned surface meshes in 3‑D, preventing mesh degeneration near triple junctions.

Numerical experiments illustrate the method’s robustness. In a 2‑D three‑phase test, two circular bubbles merge, forming a triple junction that evolves according to the prescribed angle condition; spurious currents are negligible and the volumes of the three phases remain constant up to machine precision. A 3‑D four‑phase example shows a cluster of three surfaces meeting along a line, intersecting a planar wall with a 90° contact angle; the energy decays monotonically and the interface geometry evolves smoothly without mesh tangling. Additional tests with contrasting densities and viscosities, as well as external body forces, confirm that the scheme handles realistic physical parameters.

In conclusion, the authors deliver the first fully discrete variational front‑tracking method that simultaneously guarantees unconditional energy stability, exact volume conservation, and accurate treatment of triple junctions for both 2‑D and 3‑D multiphase flows. The approach combines the flexibility of unfitted finite elements with the geometric robustness of the BGN parametric formulation, opening the way for reliable simulations of complex multiphase phenomena such as oil recovery, ink‑jet printing, microfluidic devices, and material science applications where triple (or higher‑order) junctions play a critical role. Future work may extend the framework to anisotropic surface tensions, quadruple junctions, higher‑order time integration, and large‑scale parallel implementations.