  In this paper, we present a semi-implicit numerical solver for a first order
hyperbolic formulation of two-phase flow with surface tension and viscosity.
The numerical method addresses several complexities presented by the PDE system
in consideration: (i) The presence of involution constraints of curl type in
the governing equations requires explicit enforcement of the zero-curl property
of certain vector fields (an interface field and a distortion field); the
problem is eliminated entirely by adopting a set of compatible curl and
gradient discrete differential operators on a staggered grid, allowing to
preserve the Schwartz identity of cross-derivatives exactly. (ii) The evolution
equations feature highly nonlinear stiff algebraic source terms which are used
for the description of viscous interactions as emergent behaviour of an
elasto-plastic solid in the stiff strain relaxation limit; such source terms
are reliably integrated with an efficient semi-analytical technique. (iii) In
the low-Mach number regime, standard explicit Godunov-type schemes lose
efficiency and accuracy; the issue is addressed by means of a simple
semi-implicit, pressure-based, split treatment of acoustic and non-acoustic
waves, again using staggered grids that recover the implicit solution for a
single scalar field (the pressure) through a sequence of symmetric-positive
definite linear systems that can be efficiently solved via the conjugate
gradient method.
