Parametric finite element approximation of two-phase Navier--Stokes flow with viscoelasticity

In this work, we present a parametric finite element approximation of two-phase Navier-Stokes flow with viscoelasticity. The free boundary problem is given by the viscoelastic Navier-Stokes equations in the two fluid phases, connected by jump conditions across the interface. The elasticity in the fluids is characterised using the Oldroyd-B model with possible stress diffusion. The model was originally introduced to approximate fluid-structure interaction problems between an incompressible Newtonian fluid and a hyperelastic neo-Hookean solid, which are possible limit cases of the model. We approximate a variational formulation of the model with an unfitted finite element method that uses piecewise linear parametric finite elements. The two-phase Navier-Stokes-Oldroyd-B system in the bulk regions is discretised in a way that guarantees unconditional solvability and stability for the coupled bulk-interface system. Good volume conservation properties for the two phases are observed in the case where the pressure approximation space is enriched with the help of an XFEM function. We show the applicability of our method with some numerical results.

💡 Research Summary

The paper presents a novel numerical framework for simulating two‑phase incompressible flows in which each phase exhibits viscoelastic behavior described by the Oldroyd‑B model, possibly augmented with stress diffusion. The authors combine a parametric finite‑element (FE) method for explicit interface tracking with an unfitted bulk discretisation and an extended finite‑element (XFEM) enrichment for the pressure field.

In the continuous setting, the two fluid domains Ω⁺(t) and Ω⁻(t) are separated by a smooth, closed hypersurface Γ(t). In each subdomain the Navier‑Stokes momentum and continuity equations are coupled with an evolution equation for the conformation tensor B. The elastic stress is Tₑ = G(B − I) and B satisfies a transport‑relaxation‑diffusion equation that includes the convective term, the upper‑convected derivative, a relaxation term (1/λ)(B − I) and an optional Laplacian term αΔB. Jump conditions on Γ(t) enforce continuity of normal velocity, a stress balance that includes surface tension γκν, and the kinematic condition V = u·ν.

The authors derive an energy identity for the continuous problem: the total energy (kinetic + elastic + interfacial) decays due to viscous dissipation, elastic relaxation, and stress diffusion. This identity guides the design of the discrete scheme.

The spatial discretisation uses piecewise linear bulk FE spaces for velocity, pressure and B on a mesh that does not conform to the moving interface. The interface itself is represented by a lower‑dimensional mesh that is moved each time step by the normal component of the fluid velocity; tangential degrees of freedom are handled implicitly to avoid mesh degeneration, following the Dziuk‑style approach. To capture the pressure jump across Γ(t) accurately, the pressure space is enriched with a single XFEM “bubble” function that is discontinuous across the interface. This enrichment yields exact volume conservation on the semi‑discrete level and very small volume errors in fully discrete computations.

Time integration is semi‑implicit: the nonlinear convection terms are treated explicitly, while viscous, pressure and elastic terms are implicit. Because the Oldroyd‑B system is nonlinear and coupled to the interface motion, the authors employ a fixed‑point iteration combined with a Schur‑complement reduction that isolates the bulk variables from the curvature term. This strategy enables a proof of unconditional stability (no restriction on the time step) and existence of discrete solutions, based on a discrete energy inequality mirroring the continuous one.

The paper provides a rigorous analysis: Lemma 3.7 gives a priori bounds for the discrete variables, and Theorem 3.1 establishes unconditional stability and solvability. The analysis also shows that the XFEM pressure enrichment guarantees exact mass conservation in the semi‑discrete formulation.

Numerical experiments in two and three dimensions illustrate the method’s performance. Test cases include a rising bubble, a droplet under shear, and a “Swiss‑cheese” configuration with multiple inclusions. The results demonstrate (i) accurate capture of interface shape and curvature, (ii) correct pressure jumps, (iii) excellent volume conservation (errors below 10⁻⁶), and (iv) robustness at high Weissenberg numbers (≥ 100), especially when stress diffusion (α > 0) is employed. Comparisons with level‑set and phase‑field approaches show that the parametric method requires far fewer degrees of freedom while delivering comparable or superior accuracy.

In summary, the authors deliver a fully coupled, unconditionally stable, and mass‑conserving scheme for two‑phase viscoelastic flows. The combination of a parametric interface representation, XFEM pressure enrichment, and a Schur‑complement‑based fixed‑point solver constitutes a significant advancement over existing methods, opening the way for reliable simulations of complex fluid‑structure interaction, polymer processing, and biomedical flows where both interfacial dynamics and viscoelastic effects are essential.