# An implicit DG solver for incompressible two-phase flows with an artificial compressibility formulation

Giuseppe Orlando

► To cite this version:

Giuseppe Orlando. An implicit DG solver for incompressible two-phase flows with an artificial compressibility formulation. International Journal for Numerical Methods in Fluids, 2024, 96 (12), pp.1932-1959. 10.1002/fld.5328. hal-04504070v3

HAL Id: hal-04504070

https://hal.science/hal-04504070v3

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
An implicit DG solver for incompressible two-phase
flows with an artificial compressibility formulation

Giuseppe Orlando$^{(1),*}$

$^{(1)}$ MOX - Dipartimento di Matematica, Politecnico di Milano
Piazza Leonardo da Vinci 32, 20133 Milano, Italy
giuseppe.orlando@polimi.it

**Keywords:** Navier-Stokes equations, incompressible flows, two-phase flows, artificial compressibility, Discontinuous Galerkin methods.

<footnote>*Present address:* CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91120 Palaiseau, France, giuseppe.orlando@polytechnique.edu</footnote>
Abstract

We propose an implicit Discontinuous Galerkin (DG) discretization for incompressible two-phase flows using an artificial compressibility formulation. The conservative level set (CLS) method is employed in combination with a reinitialization procedure to capture the moving interface. A projection method based on the L-stable TR-BDF2 method is adopted for the time discretization of the Navier-Stokes equations and of the level set method. Adaptive Mesh Refinement (AMR) is employed to enhance the resolution in correspondence of the interface between the two fluids. The effectiveness of the proposed approach is shown in a number of classical benchmarks. A specific analysis on the influence of different choices of the mixture viscosity is also carried out.
# 1 Introduction

Two-phase flows are common in many engineering and industrial applications. An evolving interface delimits the bulk regions of the single phases. Many techniques have been developed over the years to capture the motion of the interface. Two classes of methods are commonly used to locate the interface: interface-tracking and interface-capturing. Interface-tracking schemes employ either Arbitrary Lagrangian–Eulerian (ALE) methods on a mesh that deforms with the interface [16, 27] or marker and cell methods [53]. Interface-capturing techniques are instead based on fixed spatial grids with an interface function which captures the interface. A full survey on interface-capturing methods goes beyond the scope of this work and we refer e.g. to [37] for a review of these techniques. Interface capturing methods include the level set (LS) method [46, 47], which represents the interface as an iso-surface of the so-called level set function. Classically, the level set function is defined as the signed distance function. However, this choice leads to non conservative methods. A number of approaches have been developed to overcome this issue; in this work, we employ the conservative level set (CLS) method, originally proposed in [38, 39], and briefly summarized in Section 2.2. CLS includes a reinitialization equation to maintain the shape of the level set, which will be also discussed in Section 2.2.

Changing fluid properties, such as density and viscosity, and surface tension at the interface lead to discontinuities that make the discretization of the Navier-Stokes equations particularly challenging. The Discontinuous Galerkin (DG) method has been widely employed in the field of Computational Fluid Dynamics (CFD), see e.g. [9, 21, 31], and is a natural candidate for the discretization of the governing equations of two-phase flows. Several approaches have been proposed in the literature combining the DG method and the level set method, see among many others [22, 23, 48, 49]. In this paper, we propose an extension of the solver for single-phase incompressible Navier-Stokes equations with an artificial compressibility formulation presented in [40, 45], so as to overcome well known issues of projection methods. The time discretization is therefore based on the TR-BDF2 scheme [5, 10, 26, 45], which is a second order two-stage method. A brief review of the TR-BDF2 method will be given in Section 3, whereas we refer to [10, 26] for a detailed analysis of the scheme. The solver is implemented in the framework of the open source numerical library deal.II [2], which supports native non-conforming h-adaptation. We will exploit these capabilities to enhance the resolution in the regions close to the interface between the two fluids.

The paper is structured as follows: the model equations and their non-dimensional formulation are reviewed in Section 2. The time discretization approach is outlined and discussed in Section 3. The spatial discretization is presented in Section 4. The application of the proposed method to a number of significant benchmarks is reported in Section 5. Here, we also analyze the impact of different possible choices for the mixture viscosity when the interface undergoes large deformations. Finally, some conclusions and perspectives for future work are presented in Section 6.

# 2 The model equations

Let $\Omega \subset \mathbb{R}^d$, $2 \le d \le 3$ be a connected open bounded set with a sufficiently
smooth boundary $\partial\Omega$ and denote by $\mathbf{x}$ the spatial coordinates and by $t$ the temporal coordinate. The two fluids in $\Omega$ are considered immiscible and they are contained in the subdomains $\Omega_1(t)$ and $\Omega_2(t)$, respectively, so that $\overline{\Omega_1(t)} \cup \overline{\Omega_2(t)} = \overline{\Omega}$. The moving interface between the two fluids is denoted by $\Gamma(t)$, defined as $\Gamma(t) = \partial\Omega_1(t) \cap \partial\Omega_2(t)$. We consider the classical unsteady, isothermal, incompressible Navier-Stokes equations with gravity, which read as follows [29]:

$$ \rho(\mathbf{x}) \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] = -\nabla p + \nabla \cdot [2\mu(\mathbf{x})\mathbf{D}(\mathbf{u})] + \rho(\mathbf{x})\mathbf{g} \\ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

for $\mathbf{x} \in \Omega$, $t \in (0, T_f]$, supplied with suitable initial and boundary conditions, which will be specified in the following. Here $T_f$ is the final time, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, $\rho$ is the fluid density and $\mu$ is the dynamic viscosity. We assume that both the density and the viscosity are discontinuous functions

$$ \rho(\mathbf{x}) = \begin{cases} \rho_1 & \text{in } \Omega_1(t) \\ \rho_2 & \text{in } \Omega_2(t) \end{cases} \quad \text{and} \quad \mu(\mathbf{x}) = \begin{cases} \mu_1 & \text{in } \Omega_1(t) \\ \mu_2 & \text{in } \Omega_2(t) \end{cases} \qquad (2) $$

with $\rho_1, \rho_2, \mu_1$, and $\mu_2$ constant values. Moreover, $\mathbf{g}$ is the gravitational acceleration and $\mathbf{D}(\mathbf{u})$ denotes the symmetric part of the gradient of the velocity, defined as

$$ \mathbf{D}(\mathbf{u}) = \frac{1}{2} \left[ \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right]. \qquad (3) $$

In the following, for the sake of simplicity in the notation, we omit the explicit dependence on space and time for the different quantities. Surface tension effects are taken into accounts through the following balance of forces at the interface $\Gamma$:

$$ [\mathbf{u}]_{\Gamma} = \mathbf{0} \qquad [-p\mathbf{I} + 2\mu\mathbf{D}(\mathbf{u})]_{\Gamma} \mathbf{n}_{\Gamma} = \sigma\kappa\mathbf{n}_{\Gamma}, \qquad (4) $$

where $\mathbf{n}_{\Gamma}$ is the outward unit normal to $\Gamma$, $[\Psi]_{\Gamma} = \Psi|_{\Gamma \cap \Omega_1} - \Psi|_{\Gamma \cap \Omega_2}$ denotes the jump of $\Psi$ across the interface $\Gamma$, $\sigma$ is the constant surface tension coefficient, and $\kappa = -\nabla \cdot \mathbf{n}_{\Gamma}$ is the curvature. The first condition implies the continuity of the velocity along $\Gamma$, whereas the second condition describes the balance of forces at the interface. A common way to handle the term with surface tension is to introduce the following volumetric force [29]:

$$ \mathbf{f}_{\sigma} = \sigma\kappa\mathbf{n}_{\Gamma}\delta(\Gamma), \qquad (5) $$

where $\delta(\Gamma)$ is the Dirac delta distribution supported on the interface. Hence, system (1) can be rewritten as follows:

$$ \rho \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] = -\nabla p + \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] + \rho \mathbf{g} + \mathbf{f}_{\sigma} \\ \nabla \cdot \mathbf{u} = 0. \qquad (6) $$

A level set approach [46, 59] is employed to capture the interface $\Gamma$. The interface between the two fluids is considered sharp and is described as the zero level set of a smooth function. Hence, the following relation holds:

$$ \frac{\partial \varphi}{\partial t} + \mathbf{u} \cdot \nabla \varphi = 0, \qquad (7) $$
where $\varphi$ is the level set function. A common choice [59] is to consider as level set the signed distance function to $\Gamma$. In order to fix the notation, we consider $\varphi < 0$ in $\Omega_2$ and $\varphi > 0$ in $\Omega_1$. Therefore, we define

$$ \varphi = \begin{cases} -\mathrm{dist}(\mathbf{x}, \Gamma) & \text{if } \mathbf{x} \in \Omega_2 \\ 0 & \text{if } \mathbf{x} \in \Gamma \\ \mathrm{dist}(\mathbf{x}, \Gamma) & \text{if } \mathbf{x} \in \Omega_1 \end{cases} \qquad (8) $$

The unit normal vector can be evaluated at each point as follows [17, 46]:

$$ \mathbf{n}_{\Gamma} = \frac{\nabla \varphi}{|\nabla \varphi|}, \quad \mathbf{x} \in \Gamma, \qquad (9) $$

so that (7) is equivalent to

$$ \frac{\partial \varphi}{\partial t} + (\mathbf{u} \cdot \mathbf{n}_{\Gamma}) |\nabla \varphi| = 0. \qquad (10) $$

Relation (10) shows that the deformation of the level set function is due only to the normal component of the velocity. Moreover, we can express the density and the dynamic viscosity through the Heaviside function $H$

$$ \rho = \rho_2 + (\rho_1 - \rho_2) H(\varphi) \qquad (11) $$

$$ \mu = \mu_2 + (\mu_1 - \mu_2) H(\varphi) \qquad (12) $$

The whole system of equations reads therefore as follows:

$$ \begin{aligned} \rho \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] &= - \nabla p + \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] + \rho g + \mathbf{f}_{\sigma} \\ \nabla \cdot \mathbf{u} &= 0 \\ \frac{\partial \varphi}{\partial t} + \mathbf{u} \cdot \nabla \varphi &= 0. \end{aligned} \qquad (13) $$

System (13) can be rewritten in conservative form. First of all, thanks to the incompressibility constraint $\nabla \cdot \mathbf{u} = 0$, we can rewrite (7) as

$$ \frac{\partial \varphi}{\partial t} + \nabla \cdot (\varphi \mathbf{u}) = 0. \qquad (14) $$

Moreover, one can verify that (7), in combination with the incompressibility constraint, implies mass conservation. Indeed, we get

$$ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = \frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho &= (\rho_1 - \rho_2) \left( \frac{\partial H(\varphi)}{\partial t} + \mathbf{u} \cdot \nabla H(\varphi) \right) \\ &= (\rho_1 - \rho_2) \delta(\varphi) \left( \frac{\partial \varphi}{\partial t} + \mathbf{u} \cdot \nabla \varphi \right) = 0, \end{aligned} \qquad (15) $$

where we exploited the relation $\frac{dH(\varphi)}{d\varphi} = \delta(\varphi)$ [56], with $\delta(\varphi)$ denoting the Dirac delta distribution with support equal to the function $\varphi$ which implicitly describes the surface. It is appropriate to stress the fact that the differential operators involving the
Heaviside function $H(\varphi)$ have to be intended in a proper distributional sense. Finally, as discussed in [33], we can rewrite

$$ \mathbf{f}_{\sigma} = \nabla \cdot [\sigma (\mathbf{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}) \delta(\Gamma)], \quad (16) $$

where, once more, the divergence operator should be intended in a distributional sense. Hence, the conservative form of (13) is

$$ \begin{aligned} \frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) &= -\nabla p + \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] + \rho \mathbf{g} + \mathbf{f}_{\sigma} \\ \nabla \cdot \mathbf{u} &= 0 \\ \frac{\partial \varphi}{\partial t} + \nabla \cdot (\varphi \mathbf{u}) &= 0. \end{aligned} \quad (17) $$

The Continuum Surface Stress (CSS) approach, introduced in [33], is employed to treat density, viscosity, and surface tension term. A regularized Heaviside $H_{\varepsilon}(\varphi)$ is introduced, so as to obtain

$$ \rho \approx \rho_2 + (\rho_1 - \rho_2) H_{\varepsilon}(\varphi) \quad (18) $$

$$ \mu \approx \mu_2 + (\mu_1 - \mu_2) H_{\varepsilon}(\varphi). \quad (19) $$

It is important at this stage to point out the relation between $\delta(\Gamma)$ and $\delta(\varphi)$. As discussed in [18], the following relation holds:

$$ \delta(\Gamma) = \delta(\varphi) |\nabla \varphi|, \quad (20) $$

so that we can rewrite

$$ \mathbf{f}_{\sigma} = \sigma \kappa \mathbf{n}_{\Gamma} \delta(\varphi) |\nabla \varphi| = \nabla \cdot [\sigma (\mathbf{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}) \delta(\varphi) |\nabla \varphi|]. \quad (21) $$

Hence, the CSS approximation of the surface tension term reads as follows:

$$ \mathbf{f}_{\sigma} \approx \nabla \cdot [\sigma (\mathbf{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}) \delta_{\varepsilon}(\varphi) |\nabla \varphi|] = \nabla \cdot \left[ \sigma (\mathbf{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}) \frac{dH_{\varepsilon}}{d\varphi}(\varphi) |\nabla \varphi| \right]. \quad (22) $$

Since the seminal proposals in [13, 60] (see also the review in [24]), projection methods have become very popular for the discretization of incompressible Navier-Stokes equations. However, difficulties arise in choosing boundary conditions for the Poisson equation which is to be solved at each time step to compute the pressure. Several strategies have been proposed to deal with this issue and, in particular, high-order pressure boundary conditions were developed in [30, 31]. This approach, however, requires to rewrite the viscous linear terms as a solenoidal part, which is approximated by an explicit scheme, and as an irrotational part, which is approximate by an implicit scheme of appropriate order [30]. An alternative that allows to avoid or reduce some of these problems is the so-called artificial compressibility formulation, originally introduced in [12] and employed in [8, 45] among many others. In this formulation, the incompressibility constraint is relaxed and a time evolution equation for the pressure is introduced. Hence, we can consider more general boundary conditions for the pressure, and, in particular, full homogenous Neumann boundary conditions can be employed.
Moreover, this formulation is consistent with the asymptotic limit of the compressible equations for vanishing Mach number [44]. Hence, the proposed method can be extended to fully-compressible two-phase flows, leading to an asymptotic-preserving scheme, as we aim to show in future work. This kind of approach has been adopted for incompressible flows with variable density, see e.g. [7, 36], and we aim here to consider an artificial compressibility formulation for immiscible, isothermal two-phase flows with gravity. The model equations can be therefore rewritten as follows:

$$ 
\begin{aligned}
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) &= -\nabla p + \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] + \rho \mathbf{g} + \mathbf{f}_{\sigma} \\
\frac{1}{\rho_0 c^2} \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} &= 0 \\
\frac{\partial \varphi}{\partial t} + \nabla \cdot (\varphi \mathbf{u}) &= 0,
\end{aligned}
\quad (23) $$

where c is the artificial speed of sound and $\rho_0$ is a reference density. Finally, since we are relaxing the incompressibility constraint, we consider (7) for the level set motion, which is valid for the transport of $\varphi$ independently of the constraints on the velocity $\mathbf{u}$. Moreover, this choice is justified by the results reported in [40] for a rising bubble test case, for which a non-conservative formulation leads to less diffusion in the treatment of the interface. Hence, the final form of the system under consideration reads as follows:

$$ 
\begin{aligned}
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) &= -\nabla p + \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] + \rho \mathbf{g} + \mathbf{f}_{\sigma} \\
\frac{1}{\rho_0 c^2} \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} &= 0 \\
\frac{\partial \varphi}{\partial t} + \mathbf{u} \cdot \nabla \varphi &= 0,
\end{aligned}
\quad (24) $$

Before proceeding to describe the time and space discretization schemes, we perform a dimensional analysis to derive a non-dimensional version of system (24).

## 2.1 Dimensional analysis

In this Section, we derive a non-dimensional formulation for system (24). We denote with the symbol $*$ non-dimensional quantities. We introduce a reference length and velocity, denoted by $L_{ref}$ and $U_{ref}$, respectively, so as to obtain

$$ \mathbf{x} = L_{ref} \mathbf{x}^* \qquad \mathbf{u} = U_{ref} \mathbf{u}^* \qquad t = \frac{L_{ref}}{U_{ref}} t^*. \quad (25) $$

Moreover, we choose as reference density and viscosity those associated to the heavier fluid, which is conventionally considered in $\Omega_1$. For the sake of simplicity, we also assume $\rho_0 = \rho_1$. The reference pressure $p_{ref}$ is taken equal to $p_{ref} = \rho_1 U_{ref}^2$. Hence, we get

$$ \rho = \rho_1 \rho^* \qquad \mu = \mu_1 \mu^* \qquad p = \rho_1 U_{ref}^2 p^* \qquad \kappa = \frac{1}{L_{ref}} \kappa^* \qquad \varphi = L_{ref} \varphi^*. \quad (26) $$
Introducing the appropriate non-dimensional quantities, we obtain

$$ 
\begin{aligned}
\frac{\rho_1 U_{ref}^2}{L_{ref}} \frac{\partial^* (\rho^* \mathbf{u}^*)}{\partial^* t^*} + \frac{\rho_1 U_{ref}^2}{L_{ref}} \nabla^* \cdot (\rho^* \mathbf{u}^* \otimes \mathbf{u}^*) &= -\frac{\rho_1 U_{ref}^2}{L_{ref}} \nabla^* p^* + \frac{\mu_1 U_{ref}}{L_{ref}} \nabla^* \cdot [2\mu^* \mathbf{D}(\mathbf{u}^*)] \\
&\quad - \rho_1 \rho^* g \mathbf{k} \\
&\quad + \frac{1}{L_{ref}^2} \nabla^* \cdot [\sigma (\mathbf{I} - \mathbf{n}_\Gamma \otimes \mathbf{n}_\Gamma) \delta_\varepsilon^*(\varphi^*) |\nabla^* \varphi^*|]
\end{aligned}
$$

$$ \frac{\rho_1 U_{ref}^3}{\rho_1 L_{ref} c^2} \frac{\partial^* p^*}{\partial^* t^*} + \frac{U_{ref}}{L_{ref}} \nabla^* \cdot \mathbf{u}^* = 0 \qquad (27) $$

$$ U_{ref} \frac{\partial^* \varphi^*}{\partial^* t^*} + U_{ref} \mathbf{u}^* \cdot \nabla^* \varphi^* = 0, $$

where $\mathbf{k}$ is the upward pointing unit vector in the standard Cartesian reference frame. System (27) reduces to

$$ 
\begin{aligned}
\frac{\partial^* (\rho^* \mathbf{u}^*)}{\partial^* t^*} + \nabla^* \cdot (\rho^* \mathbf{u}^* \otimes \mathbf{u}^*) &= - \nabla^* p^* + \frac{1}{Re} \nabla^* \cdot [2\mu^* \mathbf{D}(\mathbf{u}^*)] - \frac{1}{Fr^2} \rho^* \mathbf{k} \\
&\quad + \frac{1}{We} \nabla^* \cdot [(\mathbf{I} - \mathbf{n}_\Gamma \otimes \mathbf{n}_\Gamma) \delta_\varepsilon^*(\varphi^*) |\nabla^* \varphi^*|] \\
M^2 \frac{\partial^* p^*}{\partial^* t^*} + \nabla^* \cdot \mathbf{u}^* &= 0 \qquad (28) \\
\frac{\partial^* \varphi^*}{\partial^* t^*} + \mathbf{u}^* \cdot \nabla^* \varphi^* &= 0,
\end{aligned}
$$

where

$$ Re = \frac{\rho_1 U_{ref} L_{ref}}{\mu_1} \qquad Fr = \frac{U_{ref}}{\sqrt{g L_{ref}}} \qquad We = \frac{\rho_1 U_{ref}^2 L_{ref}}{\sigma} \qquad M = \frac{U_{ref}}{c} \qquad (29) $$

denote the Reynolds number, the Froude number, the Weber number, and the Mach number, respectively. In the following, with a slight abuse of notation, we omit the symbol $*$ to mark non-dimensional quantities and we consider therefore the following system of equations:

$$ 
\begin{aligned}
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) &= - \nabla p + \frac{1}{Re} \nabla \cdot [2\mu \mathbf{D}(\mathbf{u})] \\
&\quad - \frac{1}{Fr^2} \rho \mathbf{k} + \frac{1}{We} \nabla \cdot [(\mathbf{I} - \mathbf{n}_\Gamma \otimes \mathbf{n}_\Gamma) \delta_\varepsilon(\varphi) |\nabla \varphi|] \\
M^2 \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} &= 0 \qquad (30) \\
\frac{\partial \varphi}{\partial t} + \mathbf{u} \cdot \nabla \varphi &= 0,
\end{aligned}
$$

where

$$ \rho = \frac{\rho_2}{\rho_1} + \left(1 - \frac{\rho_2}{\rho_1}\right) H_\varepsilon(\varphi) \qquad (31) $$

$$ \mu = \frac{\mu_2}{\mu_1} + \left(1 - \frac{\mu_2}{\mu_1}\right) H_\varepsilon(\varphi). \qquad (32) $$
## 2.2 The conservative level set method

The traditional level set method lacks of volume conservation properties [19]. The conservative level set (CLS) method [38, 39, 64] is a popular alternative to add conservation properties to level set schemes. The idea is to replace the signed distance function defined in (8) with a regularized Heaviside function:

$$ \phi(\mathbf{x}, t) = \frac{1}{1 + e^{-\varphi(\mathbf{x}, t)/\varepsilon}}, \qquad (33) $$

where $\varepsilon$ helps smoothing the transition of the discontinuous physical properties between the two subdomains and it is also known as interface thickness. Since

$$ \nabla \phi = \frac{1}{\varepsilon} \frac{e^{-\varphi/\varepsilon}}{(1 + e^{-\varphi/\varepsilon})^2} \nabla \varphi \qquad (34) $$

we can compute the outward unit normal $\mathbf{n}_\Gamma$ exactly as in (9). From definition (33), it follows that

$$ \Gamma(t) = \left\{ \mathbf{x} \in \bar{\Omega} : \phi(\mathbf{x}, t) = \frac{1}{2} \right\}. \qquad (35) $$

This new level set function needs to be reinitialized in order to keep the property of being a regularized Heaviside function [39]. This goal is achieved by solving the following PDE [38, 39]:

$$ \frac{\partial \phi}{\partial \tau} + \nabla \cdot (u_c \phi (1 - \phi) \mathbf{n}_\Gamma) = \nabla \cdot (\beta \varepsilon u_c (\nabla \phi \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma), \qquad (36) $$

where $\tau$ is an artificial pseudo-time variable, $u_c$ is an artificial compression velocity, and $\beta$ is a constant. It is important to notice that $\mathbf{n}_\Gamma$ does not change during the reinitialization procedure, but is computed using the initial value of the level set function. The relation (36) has been originally introduced as an intermediate step between the level set advection and the Navier-Stokes equations to keep the shape of the profile [38] and to stabilize the advection [39]. Two fluxes are considered: a compression flux which acts where $0 < \phi < 1$ and in normal direction to the interface, represented by $u_c\phi(1 - \phi)\mathbf{n}_\Gamma$, and a diffusion flux, represented by $\beta\varepsilon u_c (\nabla\phi \cdot \mathbf{n}_\Gamma)\mathbf{n}_\Gamma$. The reinitialization is crucial for the overall stability of the algorithm, but it also introduces errors in the solution [39, 48]. Hence, it is important to avoid unnecessary reinitialization. For this purpose, unlike the formulation proposed e.g. in [39] and [48], we introduce the coefficient $\beta$ to tune the amount of diffusion so as to keep it as small as possible. The choices for the different parameters will be specified in Section 5. Finally, we stress the fact that, in this method, we are already using a smooth version of Heaviside function so that

$$ H_{\varepsilon} = \phi \qquad (37) $$

$$ \delta(\Gamma) \approx \frac{dH_{\varepsilon}}{d\phi} |\nabla\phi| = |\nabla\phi| \qquad (38) $$

## 3 The time discretization

In this Section, we outline the time discretization strategy for system (30). Our goal here is to extend the projection method based on the TR-BDF2 scheme developed in
[45]. We now briefly recall for the convenience of the reader the formulation of the TR-BDF2. This second order implicit method has been originally introduced in [5] as a combination of the Trapezoidal Rule (or Crank-Nicolson) method and of the Backward Differentiation Formula method of order 2 (BDF2). Let $\Delta t = T_f/N$ be a discrete time step and $t^n = n\Delta t, n = 0, \dots, N$, be discrete time levels for a generic time dependent problem $\mathbf{u}' = \mathcal{N}(\mathbf{u})$. Hence, the incremental form of the TR-BDF2 scheme can be described in terms of two stages, the first one from $t^n$ to $t^{n+\gamma} = t^n + \gamma\Delta t$, and the second one from $t^{n+\gamma}$ to $t^{n+1}$, as follows:

$$ \frac{\mathbf{u}^{n+\gamma} - \mathbf{u}^n}{\gamma \Delta t} = \frac{1}{2} \mathcal{N}(\mathbf{u}^{n+\gamma}) + \frac{1}{2} \mathcal{N}(\mathbf{u}^n) \qquad (39) $$

$$ \frac{\mathbf{u}^{n+1} - \mathbf{u}^{n+\gamma}}{(1-\gamma)\Delta t} = \frac{1}{2-\gamma}\mathcal{N}(\mathbf{u}^{n+1}) + \frac{1-\gamma}{2(2-\gamma)}\mathcal{N}(\mathbf{u}^{n+\gamma}) + \frac{1-\gamma}{2(2-\gamma)}\mathcal{N}(\mathbf{u}^n). \quad (40) $$

Here, $\mathbf{u}^n$ denotes the approximation at time $n = 0, \dots, N$. Notice that, in order to guarantee L-stability, one has to choose $\gamma = 2 - \sqrt{2}$ [26]. We refer to [10, 26] for a more exhaustive discussion on the TR-BDF2 method.

We start by considering the equation in system (30) associated to the level set. In order to avoid a full coupling with the Navier-Stokes equations, we perform a linearization in velocity, so that the first stage for the level set update reads as follows:

$$ \frac{\phi^{n+\gamma} - \phi^n}{\gamma \Delta t} + \frac{1}{2} \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \phi^{n+\gamma} = -\frac{1}{2} \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \phi^n, \quad (41) $$

where the approximation $\mathbf{u}^{n+\frac{\gamma}{2}}$ is defined by extrapolation as

$$ \mathbf{u}^{n+\frac{\gamma}{2}} = \left(1 + \frac{\gamma}{2(1-\gamma)}\right) \mathbf{u}^n - \frac{\gamma}{2(1-\gamma)} \mathbf{u}^{n-1}. \quad (42) $$

Following then the projection approach described in [14, 45] and applying (39), the momentum predictor equation for the first stage reads as follows:

$$ \begin{aligned} & \frac{\rho^{n+\gamma} \mathbf{u}^{n+\gamma,*} - \rho^n \mathbf{u}^n}{\gamma \Delta t} \\ & + \frac{1}{2} \nabla \cdot (\rho^{n+\gamma} \mathbf{u}^{n+\gamma,*} \otimes \mathbf{u}^{n+\frac{\gamma}{2}}) - \frac{1}{2} \frac{1}{Re} \nabla \cdot [2\mu^{n+\gamma} \mathbf{D}(\mathbf{u}^{n+\gamma,*})] = \\ & - \frac{1}{2} \nabla \cdot (\rho^n \mathbf{u}^n \otimes \mathbf{u}^{n+\frac{\gamma}{2}}) + \frac{1}{2} \frac{1}{Re} \nabla \cdot [2\mu^n \mathbf{D}(\mathbf{u}^n)] - \nabla p^n \\ & + \frac{1}{2} \frac{1}{We} \nabla \cdot \left[ (\mathbf{I} - \mathbf{n}_\Gamma^{n+\gamma} \otimes \mathbf{n}_\Gamma^{n+\gamma}) \delta_\varepsilon(\phi^{n+\gamma}) |\nabla \phi^{n+\gamma}| \right] \quad (43) \\ & + \frac{1}{2} \frac{1}{We} \nabla \cdot \left[ (\mathbf{I} - \mathbf{n}_\Gamma^n \otimes \mathbf{n}_\Gamma^n) \delta_\varepsilon(\phi^n) |\nabla \phi^n| \right] \\ & - \frac{1}{2} \frac{1}{Fr^2} \rho^{n+\gamma} \mathbf{k} - \frac{1}{2} \frac{1}{Fr^2} \rho^n \mathbf{k}. \end{aligned} $$

Notice once more that, in order to avoid solving a non-linear system at each time step, $\mathbf{u}^{n+\frac{\gamma}{2}}$ is employed in the momentum advection terms. We set then $\delta p^{n+\gamma} = p^{n+\gamma} - p^n$ and impose

$$ \begin{aligned} & \frac{\rho^{n+\gamma} \mathbf{u}^{n+\gamma} - \mathbf{u}^{n+\gamma,*}}{\gamma \Delta t} = - \nabla \delta p^{n+\gamma} \\ & M^2 \frac{\delta p^{n+\gamma}}{\gamma \Delta t} + \nabla \cdot \mathbf{u}^{n+\gamma} = 0. \end{aligned} \quad (44) $$
Substituting the first equation into the second in (44), one obtains the Helmholtz equation

$$M^2 \frac{\delta p^{n+\gamma}}{\gamma^2 \Delta t^2} - \nabla \cdot \left( \frac{\nabla \delta p^{n+\gamma}}{\rho^{n+\gamma}} \right) = - \frac{1}{\gamma \Delta t} \nabla \cdot \mathbf{u}^{n+\gamma,*}. \quad (45)$$

Once this equation is solved, the final velocity update for the first stage is given by

$$\mathbf{u}^{n+\gamma} = \mathbf{u}^{n+\gamma,*} - \gamma \Delta t \frac{\nabla \delta p^{n+\gamma}}{\rho^{n+\gamma}}. \quad (46)$$

The second TR-BDF2 stage is performed in a similar manner applying (40). We first focus on the level set update:

$$\frac{\phi^{n+1} - \phi^{n+\gamma}}{(1-\gamma)\Delta t} + a_{33} \mathbf{u}^{n+\frac{3}{2}\gamma} \cdot \nabla \phi^{n+1} = -a_{32} \mathbf{u}^{n+\gamma} \cdot \nabla \phi^{n+\gamma} - a_{31} \mathbf{u}^n \cdot \nabla \phi^n, \quad (47)$$

where

$$a_{31} = \frac{1-\gamma}{2(2-\gamma)}, \quad a_{32} = \frac{1-\gamma}{2(2-\gamma)}, \quad a_{33} = \frac{1}{2-\gamma}. \quad (48)$$

Again, in order to avoid a full coupling with the Navier-Stokes equations, an approximation is introduced in the advection term, so that $\mathbf{u}^{n+\frac{3}{2}\gamma}$ is defined by extrapolation as

$$\mathbf{u}^{n+\frac{3}{2}\gamma} = \left(1 + \frac{1+\gamma}{\gamma}\right) \mathbf{u}^{n+\gamma} - \frac{1-\gamma}{\gamma} \mathbf{u}^n. \quad (49)$$

Then, we define the second momentum predictor:

$$\begin{aligned} & \frac{\rho^{n+1} \mathbf{u}^{n+1,*} - \rho^{n+\gamma} \mathbf{u}^{n+\gamma}}{(1-\gamma)\Delta t} \\ & + a_{33} \nabla \cdot \left( \rho^{n+1} \mathbf{u}^{n+1,*} \otimes \mathbf{u}^{n+\frac{3}{2}\gamma} \right) - a_{33} \frac{1}{Re} \nabla \cdot \left[ 2\mu^{n+1} \mathbf{D}(\mathbf{u}^{n+1,*}) \right] = \\ & - a_{32} \nabla \cdot \left( \rho^{n+\gamma} \mathbf{u}^{n+\gamma} \otimes \mathbf{u}^{n+\gamma} \right) + a_{32} \frac{1}{Re} \nabla \cdot \left[ 2\mu^{n+\gamma} \mathbf{D}(\mathbf{u}^{n+\gamma}) \right] \\ & - a_{31} \nabla \cdot \left( \rho^n \mathbf{u}^n \otimes \mathbf{u}^n \right) + a_{31} \frac{1}{Re} \nabla \cdot \left[ 2\mu^n \mathbf{D}(\mathbf{u}^n) \right] - \nabla p^{n+\gamma} \\ & + a_{33} \frac{1}{We} \nabla \cdot \left[ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^{n+1} \otimes \mathbf{n}_{\Gamma}^{n+1} \right) \delta_{\varepsilon}(\phi^{n+1}) \left| \nabla \phi^{n+1} \right| \right] \\ & + a_{32} \frac{1}{We} \nabla \cdot \left[ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^{n+\gamma} \otimes \mathbf{n}_{\Gamma}^{n+\gamma} \right) \delta_{\varepsilon}(\phi^{n+\gamma}) \left| \nabla \phi^{n+\gamma} \right| \right] \\ & + a_{31} \frac{1}{We} \nabla \cdot \left[ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^n \otimes \mathbf{n}_{\Gamma}^n \right) \delta_{\varepsilon}(\phi^n) \left| \nabla \phi^n \right| \right] \\ & - a_{33} \frac{1}{Fr^2} \rho^{n+1} \mathbf{k} - a_{32} \frac{1}{Fr^2} \rho^{n+\gamma} \mathbf{k} - a_{31} \frac{1}{Fr^2} \rho^n \mathbf{k}. \end{aligned} \quad (50)$$

Notice that $\mathbf{u}^{n+\frac{3}{2}\gamma}$ is employed in the non-linear momentum advection term. We set then $\delta p^{n+1} = p^{n+1} - p^{n+\gamma}$ and impose

$$\begin{aligned} & \rho^{n+1} \frac{\mathbf{u}^{n+1} - \mathbf{u}^{n+1,*}}{(1-\gamma)\Delta t} = -\nabla \delta p^{n+1} \\ & M^2 \frac{\delta p^{n+1}}{(1-\gamma)\Delta t} + \nabla \cdot \mathbf{u}^{n+1} = 0. \end{aligned} \quad (51)$$
Substituting the first equation into the second in (51), one obtains the Helmholtz equation

$$M^2 \frac{\delta p^{n+1}}{(1-\gamma)^2 \Delta t^2} - \nabla \cdot \left( \frac{\nabla \delta p^{n+1}}{\rho^{n+1}} \right) = - \frac{1}{(1-\gamma)\Delta t} \nabla \cdot \mathbf{u}^{n+1,*}. \quad (52)$$

The final velocity update then reads as follows:

$$\mathbf{u}^{n+1} = \mathbf{u}^{n+1,*} - (1-\gamma)\Delta t \frac{\nabla \delta p^{n+1}}{\rho^{n+1}}. \quad (53)$$

Finally, we focus on the reinitialization procedure described in Equation 36, which is performed after each stage of the level set update and before computing the momentum predictor. We consider an implicit treatment of the diffusion term $\nabla \cdot (\beta \varepsilon u_c (\nabla \phi \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma)$ and a semi-implicit treatment of the compression term $u_c \nabla \cdot (\phi (1-\phi) \mathbf{n}_\Gamma)$. Hence, the semi-discrete formulation reads as follows:

$$\frac{\phi^{k+1,*} - \phi^{k,*}}{\Delta \tau} + \nabla \cdot (u_c \phi^{k+1,*} (1 - \phi^{k,*}) \mathbf{n}_\Gamma) = \nabla \cdot (\beta \varepsilon u_c (\nabla \phi^{k+1,*} \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma), \quad (54)$$

where $\Delta \tau$ is the pseudo time step. Moreover, $\phi^{0,*} = \phi^{n+\gamma}$ after the first TR-BDF2 stage and $\phi^{0,*} = \phi^{n+1}$ after the second TR-BDF2 stage. We recall once more that $\mathbf{n}_\Gamma = \frac{\nabla \phi^{0,*}}{|\nabla \phi^{0,*}|}$ and it does not change during the reinitialization. Following [48], we define the total reinitialization time $\tau_{fin}$ as a fraction of the time step $\Delta t$, namely

$$\tau_{fin} = \eta \Delta t. \quad (55)$$

$\eta = 0$ corresponds to no reinitialization, whereas $\eta = 1$ yields an amount of reinitialization which can modify the values of level set function of the same order of magnitude of which they have been modified during the previous advection step. For most applications, $\eta \approx 0.5$ seems to provide an appropriate amount of reinitialization [48]. A pseudo time step such that two to five reinitialization steps are performed typically ensures stable solutions and leads to the updated level set function [32].

# 4 The spatial discretization

For the spatial discretization, we consider discontinuous finite element approximations. We consider a decomposition of the domain $\Omega$ into a family of hexahedra $\mathcal{T}_h$ (quadrilaterals in the two-dimensional case) and denote each element by $K$. The skeleton $\mathcal{E}$ denotes the set of all element faces and $\mathcal{E} = \mathcal{E}^I \cup \mathcal{E}^B$, where $\mathcal{E}^I$ is the subset of interior faces and $\mathcal{E}^B$ is the subset of boundary faces. Suitable jump and average operators can then be defined as customary for finite element discretizations. A face $e \in \mathcal{E}^I$ shares two elements that we denote by $K^+$ with outward unit normal $\mathbf{n}^+$ and $K^-$ with outward unit normal $\mathbf{n}^-$, whereas for a face $e \in \mathcal{E}^B$ we denote by $\mathbf{n}$ the outward unit normal. For a scalar function $\Psi$ the jump is defined as

$$[[\Psi]] = \Psi^+ \mathbf{n}^+ + \Psi^- \mathbf{n}^- \quad \text{if } e \in \mathcal{E}^I \qquad [[\Psi]] = \Psi \mathbf{n} \quad \text{if } e \in \mathcal{E}^B. \quad (56)$$

The average is defined as

$$\{\{\Psi\}\} = \frac{1}{2}(\Psi^+ + \Psi^-) \quad \text{if } e \in \mathcal{E}^I \qquad \{\{\Psi\}\} = \Psi \quad \text{if } e \in \mathcal{E}^B. \quad (57)$$
Similar definitions apply for a vector function $\Psi$:

$$ [[\Psi]] = \Psi^{+} \cdot \mathbf{n}^{+} + \Psi^{-} \cdot \mathbf{n}^{-} \quad \text{if } e \in \mathcal{E}^{I} \qquad [[\Psi]] = \Psi \cdot \mathbf{n} \quad \text{if } e \in \mathcal{E}^{B} \qquad (58) $$

$$ \{\{\Psi\}\} = \frac{1}{2}(\Psi^{+} + \Psi^{-}) \quad \text{if } e \in \mathcal{E}^{I} \qquad \{\{\Psi\}\} = \Psi \quad \text{if } e \in \mathcal{E}^{B}. \qquad (59) $$

For vector functions, it is also useful to define a tensor jump as:

$$ \langle\langle\Psi\rangle\rangle = \Psi^{+} \otimes \mathbf{n}^{+} + \Psi^{-} \otimes \mathbf{n}^{-} \quad \text{if } \Gamma \in \mathcal{E}^{I} \qquad \langle\langle\Psi\rangle\rangle = \Psi \otimes \mathbf{n} \quad \text{if } \Gamma \in \mathcal{E}^{B}. \qquad (60) $$

We now introduce the following finite element spaces:

$$ Q_k = \{v \in L^2(\Omega) : v|_{K} \in Q_k \quad \forall K \in \mathcal{T}_h\} $$

and

$$ Q_k = [Q_k]^d, $$

where $Q_k$ is the space of polynomials of degree $k$ in each coordinate direction. Considering the well-posedness analyses in [57, 61], the finite element spaces that will be used for the discretization of velocity and pressure are $V_h = Q_k$ and $W_h = Q_{k-1} \cap L_0^2(\Omega)$, respectively, where $k \ge 2$. For what concerns the level set function, we consider instead $X_h = Q_r$ with $r \ge 2$, so that its gradient is at least a piecewise linear polynomial. We then denote by $\psi_i(\mathbf{x})$ the basis functions for the finite element spaces associated to the scalar variable, i.e. $W_h$ and $X_h$, and by $\psi_i(\mathbf{x})$ the basis functions for the space $V_h$, the finite element space chosen for the discretization of the velocity. Hence, we get

$$ \mathbf{u} \approx \sum_{j=1}^{\dim(V_h)} u_j(t) \psi_j(\mathbf{x}) \qquad p \approx \sum_{j=1}^{\dim(W_h)} p_j(t) \psi_j(\mathbf{x}) \qquad \phi \approx \sum_{j=1}^{\dim(X_h)} \phi_j(t) \psi_j(\mathbf{x}) \qquad (61) $$

The shape functions correspond to the products of Lagrange interpolation polynomials for the support points of $(k+1)$-order Gauss-Lobatto quadrature rule in each coordinate direction. Finally, as custom for discontinuous finite elements, we employ a weak imposition of the boundary conditions [4]. More specifically, we deduce suitable exterior values from the boundary conditions, employing the so-called mirror principle. For Dirichlet boundary conditions, we set (e.g. for the velocity) $\mathbf{u}^{-} = -\mathbf{u}^{+} + 2\mathbf{u}_{D}$ with $[\nabla \mathbf{u}^{+}] \mathbf{n} = [\nabla \mathbf{u}^{-}] \mathbf{n}$ and $\mathbf{u}_{D}$ denoting the Dirichlet value. On the other hand, for Neumann boundary conditions, we set (e.g. for the pressure) $\nabla p^{-} \cdot \mathbf{n} = -\nabla p^{+} \cdot \mathbf{n} + 2p_{N}$ with $p^{+} = p^{-}$ and $p_{N}$ denoting the Neumann value. Given these definitions, the weak formulation of the level set update for the first stage is obtained multiplying equation
(41) by a test function $w \in X_h$:

$$
\begin{align}
& \sum_{K \in \mathcal{T}_h} \int_K \frac{\phi^{n+\gamma}}{\gamma \Delta t} w d\Omega + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \phi^{n+\gamma} w d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \phi^{n+\gamma} \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[w]] d\Sigma - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[\phi^{n+\gamma} w]] d\Sigma \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[\phi^{n+\gamma}]] \cdot [[w]] d\Sigma \nonumber \\
&= \sum_{K \in \mathcal{T}_h} \int_K \frac{\phi^n}{\gamma \Delta t} w d\Omega - \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \phi^n w d\Omega \tag{62} \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \phi^n \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[w]] d\Sigma - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[\phi^n w]] d\Sigma \nonumber \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[\phi^n]] \cdot [[w]] d\Sigma, \nonumber
\end{align}
$$

where

$$
\lambda^{n+\frac{\gamma}{2}} = \max \left( \left| \left( \mathbf{u}^{n+\frac{\gamma}{2}} \right)^{+} \cdot \mathbf{n}^{+} \right|, \left| \left( \mathbf{u}^{n+\frac{\gamma}{2}} \right)^{-} \cdot \mathbf{n}^{-} \right| \right). \quad (63)
$$

Following [9], the numerical approximation of the non-conservative term is based on a
double integration by parts. The algebraic form can be obtained taking $w = \psi_i$, $i =
1, \dots, \dim(X_h)$ and exploiting the representation in (61), so as to obtain in compact
form

$$
\left( \frac{1}{\gamma \Delta t} \mathbf{M}_{\phi} + \frac{1}{2} \mathbf{A}_{\phi}^{n+\gamma} \right) \boldsymbol{\phi}^{n+\gamma} = \mathbf{F}_{\phi}^{n}, \quad (64)
$$

where $\boldsymbol{\phi}^{n+\gamma}$ denotes the vector of the degrees of freedom associated to the level set.
Moreover, we have set

$$
\mathbf{M}_{\phi_{ij}} = \sum_{K \in \mathcal{T}_h} \int_K \psi_j \psi_i d\Omega \qquad (65)
$$

$$
\begin{align}
\mathbf{A}_{\phi_{ij}}^{n+\gamma} &= \sum_{K \in \mathcal{T}_h} \int_K \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \psi_j \psi_i d\Omega \nonumber \\
&\quad + \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \mathbf{u}^{n+\frac{\gamma}{2}} \psi_j \right\} \right\} \cdot [[\psi_i]] d\Sigma - \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[\psi_j \psi_i]] d\Sigma \nonumber \\
&\quad + \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[\psi_j]] \cdot [[\psi_i]] d\Sigma \tag{66}
\end{align}
$$
and

$$
\begin{align}
\mathbf{F}_{\phi}^{n} &= \sum_{K \in \mathcal{T}_h} \int_K \frac{\phi^n}{\gamma \Delta t} \psi_i d\Omega + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla \phi^n \psi_i d\Omega \tag{67} \\
&\quad - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \phi^n \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[\psi_i]] d\Sigma + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} \cdot [[\phi^n \psi_i]] d\Sigma \nonumber \\
&\quad - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[\phi^n]] \cdot [[\psi_i]] d\Sigma. \nonumber
\end{align}
$$

Consider now the variational formulation for equation (43). Take $\mathbf{v} \in \mathbf{V}_h$ so as to obtain after integration by parts

$$
\begin{align}
& \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{\gamma \Delta t} \rho^{n+\gamma} \mathbf{u}^{n+\gamma,*} \cdot \mathbf{v} d\Omega - \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \rho^{n+\gamma} \mathbf{u}^{n+\gamma,*} \otimes \mathbf{u}^{n+\frac{\gamma}{2}} : \nabla \mathbf{v} d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \rho^{n+\gamma} \mathbf{u}^{n+\gamma,*} \otimes \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} \left\langle\left\langle \rho^{n+\gamma} \mathbf{u}^{n+\frac{\gamma}{2}} \right\rangle\right\rangle : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
& + \frac{1}{2Re} \sum_{K \in \mathcal{T}_h} \int_K 2\mu^{n+\gamma} \mathbf{D}(\mathbf{u}^{n+\gamma,*}) : \nabla \mathbf{v} d\Omega - \frac{1}{2Re} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ 2\mu^{n+\gamma} \mathbf{D}(\mathbf{u}^{n+\gamma,*}) \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
& - \frac{1}{2Re} \sum_{e \in \mathcal{E}} \int_e \left\langle\left\langle \mathbf{u}^{n+\gamma,*} \right\rangle\right\rangle : \left\{ \left\{ 2\mu^{n+\gamma} \mathbf{D}(\mathbf{v}) \right\} \right\} d\Sigma \nonumber \\
& + \frac{1}{2Re} \sum_{e \in \mathcal{E}} \int_e C_u \left\{ \left\{ \mu^{n+\gamma} \right\} \right\}_H \left\langle\left\langle \mathbf{u}^{n+\gamma,*} \right\rangle\right\rangle : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
&= \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{\gamma \Delta t} \rho^n \mathbf{u}^n \cdot \mathbf{v} d\Omega + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \rho^n \mathbf{u}^n \otimes \mathbf{u}^{n+\frac{\gamma}{2}} : \nabla \mathbf{v} d\Omega \tag{68} \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \rho^n \mathbf{u}^n \otimes \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma - \frac{1}{2} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} \left\langle\left\langle \rho^n \mathbf{u}^n \right\rangle\right\rangle : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
& - \frac{1}{2Re} \sum_{K \in \mathcal{T}_h} \int_K 2\mu^n \mathbf{D}(\mathbf{u}^n) : \nabla \mathbf{v} d\Omega + \frac{1}{2Re} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ 2\mu^n \mathbf{D}(\mathbf{u}^n) \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
& + \sum_{K \in \mathcal{T}_h} \int_K p^n \nabla \cdot \mathbf{v} d\Omega - \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ p^n \right\} \right\} [[\mathbf{v}]] d\Sigma \nonumber \\
& - \frac{1}{2Fr^2} \sum_{K \in \mathcal{T}_h} \int_K \rho^{n+\gamma} \mathbf{k} \cdot \mathbf{v} d\Omega - \frac{1}{2Fr^2} \sum_{K \in \mathcal{T}_h} \int_K \rho^n \mathbf{k} \cdot \mathbf{v} d\Omega \nonumber \\
& - \frac{1}{2We} \sum_{K \in \mathcal{T}_h} \int_K \left( \mathbf{I} - \mathbf{n}_\Gamma^{n+\gamma} \otimes \mathbf{n}_\Gamma^{n+\gamma} \right) \delta_\varepsilon(\phi^{n+\gamma}) |\nabla \phi^{n+\gamma}| : \nabla \mathbf{v} d\Omega \nonumber \\
& + \frac{1}{2We} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_\Gamma^{n+\gamma} \otimes \mathbf{n}_\Gamma^{n+\gamma} \right) \delta_\varepsilon(\phi^{n+\gamma}) |\nabla \phi^{n+\gamma}| \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma \nonumber \\
& - \frac{1}{2We} \sum_{K \in \mathcal{T}_h} \int_K \left( \mathbf{I} - \mathbf{n}_\Gamma^n \otimes \mathbf{n}_\Gamma^n \right) \delta_\varepsilon(\phi^n) |\nabla \phi^n| : \nabla \mathbf{v} d\Omega \nonumber \\
& + \frac{1}{2We} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_\Gamma^n \otimes \mathbf{n}_\Gamma^n \right) \delta_\varepsilon(\phi^n) |\nabla \phi^n| \right\} \right\} : \langle\langle \mathbf{v} \rangle\rangle d\Sigma, \nonumber
\end{align}
$$
where

$$ \{\{\mu^{n+\gamma}\}\}_H = \frac{2}{\frac{1}{\mu^{n+\gamma,+}} + \frac{1}{\mu^{n+\gamma,-}}}. \qquad (69) $$

Here, following e.g. [1], we employ the harmonic average of the viscosity coefficient for the penalization term. Notice that the approximation of the advection term employs an upwind flux, whereas the approximation of the diffusion term is based on the Symmetric Interior Penalty (SIP) [3]. Notice also that no penalization terms have been introduced for the variables computed at previous time steps in the diffusion terms. Following [20, 45], we set for each face e of a cell K

$$ \sigma_{e,K}^{\mathbf{u}} = (k+1)^2 \frac{\mathrm{diam}(e)}{\mathrm{diam}(K)} \qquad (70) $$

and we define the penalization constant for the SIP method as

$$ C_u = \frac{1}{2} (\sigma_{e,K^+}^{\mathbf{u}} + \sigma_{e,K^-}^{\mathbf{u}}) \quad \text{if } e \in \mathcal{E}^I, \quad C_u = \sigma_{e,K}^{\mathbf{u}} \quad \text{if } e \in \mathcal{E}^B. \qquad (71) $$

Finally, we stress the fact that a centered flux has been employed for the surface tension terms. The algebraic formulation is then computed considering $\mathbf{v} = \psi_i$, $i = 1, \dots, \dim(\mathbf{V}_h)$ and the representation in (61) for the velocity. Hence, we obtain

$$ \left( \frac{1}{\gamma \Delta t} \mathbf{M}_{\mathbf{u}}^{n+\gamma} + \frac{1}{2Re} \mathbf{A}_{\mathbf{u}}^{n+\gamma} + \frac{1}{2} \mathbf{C}_{\mathbf{u}}^{n+\gamma} \right) \mathbf{U}^{n+\gamma,*} = \mathbf{F}_{\mathbf{u}}^{n}, \qquad (72) $$

where $\mathbf{U}^{n+\gamma,*}$ denotes the vector of degrees of freedom for the velocity. Moreover, we have set

$$ \mathbf{M}_{\mathbf{u}_{ij}}^{n+\gamma} = \sum_{K \in \mathcal{T}_h} \int_K \rho^{n+\gamma} \psi_j \psi_i d\Omega \qquad (73) $$

$$ \begin{aligned} \mathbf{C}_{\mathbf{u}_{ij}}^{n+\gamma} = & - \sum_{K \in \mathcal{T}_h} \int_K \rho^{n+\gamma} \psi_j \otimes \mathbf{u}^{n+\frac{\gamma}{2}} : \nabla \psi_i d\Omega \\ & + \sum_{e \in \mathcal{E}} \int_e \{ \{ \rho^{n+\gamma} \psi_j \otimes \mathbf{u}^{n+\frac{\gamma}{2}} \} \} : \langle \langle \psi_j \rangle \rangle d\Sigma \quad (74) \\ & + \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} \langle \langle \rho^{n+\gamma} \psi_j \rangle \rangle : \langle \langle \psi_i \rangle \rangle d\Sigma \end{aligned} $$

$$ \begin{aligned} \mathbf{A}_{\mathbf{u}_{ij}}^{n+\gamma} = & \sum_{K \in \mathcal{T}_h} \int_K 2\mu^{n+\gamma} \mathbf{D}(\psi_j) : \nabla \psi_i d\Omega \\ & - \sum_{e \in \mathcal{E}} \int_e \{ \{ 2\mu^{n+\gamma} \mathbf{D}(\psi_j) \} \} : \langle \langle \psi_i \rangle \rangle d\Sigma \\ & - \sum_{e \in \mathcal{E}} \int_e \langle \langle \psi_j^{n+\gamma,*} \rangle \rangle : \{ \{ 2\mu^{n+\gamma} \mathbf{D}(\psi_i) \} \} d\Sigma \\ & + \sum_{e \in \mathcal{E}} \int_e C_u \{ \{ \mu^{n+\gamma} \} \}_H \langle \langle \psi_j \rangle \rangle : \langle \langle \psi_i \rangle \rangle d\Sigma \end{aligned} \qquad (75) $$
and

$$
\begin{align}
\mathbf{F}_{\mathbf{u}}^{n} &= \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{\gamma \Delta t} \rho^n \mathbf{u}^n \cdot \psi_i d\Omega + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \rho^n \mathbf{u}^n \otimes \mathbf{u}^{n+\frac{\gamma}{2}} : \nabla \psi_i d\Omega \nonumber \\
&\quad - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \rho^n \mathbf{u}^n \otimes \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} : \langle\langle \psi_i \rangle\rangle d\Sigma \nonumber \\
&\quad - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} \langle\langle \rho^n \mathbf{u}^n \rangle\rangle : \langle\langle \psi_i \rangle\rangle d\Sigma \nonumber \\
&\quad - \frac{1}{2Re} \sum_{K \in \mathcal{T}_h} \int_K 2\mu^n \mathbf{D}(\mathbf{u}^n) : \nabla \psi_i d\Omega + \frac{1}{2Re} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ 2\mu^n \mathbf{D}(\mathbf{u}^n) \right\} \right\} : \langle\langle \psi_i \rangle\rangle d\Sigma \nonumber \\
&\quad + \sum_{K \in \mathcal{T}_h} \int_K p^n \nabla \cdot \psi_i d\Omega - \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ p^n \right\} \right\} \left[ \langle\psi_i\rangle \right] d\Sigma \tag{76} \\
&\quad - \frac{1}{2Fr^2} \sum_{K \in \mathcal{T}_h} \int_K \rho^{n+\gamma} \mathbf{k} \cdot \psi_i d\Omega - \frac{1}{2Fr^2} \sum_{K \in \mathcal{T}_h} \int_K \rho^n \mathbf{k} \cdot \psi_i d\Omega \nonumber \\
&\quad - \frac{1}{2We} \sum_{K \in \mathcal{T}_h} \int_K \left( \mathbf{I} - \mathbf{n}_{\Gamma}^{n+\gamma} \otimes \mathbf{n}_{\Gamma}^{n+\gamma} \right) \delta_{\varepsilon}(\phi^{n+\gamma}) |\nabla \phi^{n+\gamma}| : \nabla \psi_i d\Omega \nonumber \\
&\quad + \frac{1}{2We} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^{n+\gamma} \otimes \mathbf{n}_{\Gamma}^{n+\gamma} \right) \delta_{\varepsilon}(\phi^{n+\gamma}) |\nabla \phi^{n+\gamma}| \right\} \right\} : \langle\langle \psi_i \rangle\rangle d\Sigma \nonumber \\
&\quad - \frac{1}{2We} \sum_{K \in \mathcal{T}_h} \int_K \left( \mathbf{I} - \mathbf{n}_{\Gamma}^n \otimes \mathbf{n}_{\Gamma}^n \right) \delta_{\varepsilon}(\phi^n) |\nabla \phi^n| : \nabla \psi_i d\Omega \nonumber \\
&\quad + \frac{1}{2We} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^n \otimes \mathbf{n}_{\Gamma}^n \right) \delta_{\varepsilon}(\phi^n) |\nabla \phi^n| \right\} \right\} : \langle\langle \psi_i \rangle\rangle d\Sigma. \nonumber
\end{align}
$$

For what concerns the projection step, we apply again the SIP method. We multiply (45) by a test function $q \in Q_h$, we apply Green's theorem and we get

$$
\begin{align}
& \sum_{K \in \mathcal{T}_h} \int_K \frac{M^2}{\gamma^2 \Delta t^2} \delta p^{n+\gamma} q d\Omega + \sum_{K \in \mathcal{T}_h} \int_K \frac{\nabla \delta p^{n+\gamma}}{\rho^{n+\gamma}} \cdot \nabla q d\Omega \nonumber \\
& - \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \frac{\nabla \delta p^{n+\gamma}}{\rho^{n+\gamma}} \right\} \right\} \cdot [[q]] d\Sigma - \sum_{e \in \mathcal{E}} \int_e \left[ [[\delta p^{n+\gamma}]] \right] \cdot \left\{ \left\{ \frac{\nabla q}{\rho^{n+\gamma}} \right\} \right\} d\Sigma \nonumber \\
& + \sum_{e \in \mathcal{E}} \int_e C_p \left\{ \left\{ \frac{1}{\rho^{n+\gamma}} \right\} \right\}_H \cdot [[\delta p^{n+\gamma}]] \cdot [[q]] d\Sigma \tag{77} \\
& = \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{\gamma \Delta t} \mathbf{u}^{n+\gamma, **} \cdot \nabla q d\Omega - \sum_{e \in \mathcal{E}} \int_e \frac{1}{\gamma \Delta t} \left\{ \left\{ \mathbf{u}^{n+\gamma, *} \right\} \right\} \cdot [[q]] d\Sigma, \nonumber
\end{align}
$$

where we set

$$
\sigma_{e,K}^{p} = k^2 \frac{\operatorname{diam}(e)}{\operatorname{diam}(K)}, \qquad (78)
$$

so that

$$
C_p = \frac{1}{2} (\sigma_{e,K^+}^p + \sigma_{e,K^-}^p) \quad \text{if } e \in \mathcal{E}^I, \qquad C_p = \sigma_{e,K}^p \quad \text{if } e \in \mathcal{E}^B. \tag{79}
$$
The algebraic formulation is once more obtained taking $q = \psi_i, i = 1, \dots, \dim(W_h)$ and considering the expansion for $p^{n+\gamma}$ reported in (61). Hence, we get

$$ \left( \frac{M^2}{\gamma^2 \Delta t^2} \mathbf{M}_p^{n+\gamma} + \mathbf{K}_p \right) \mathbf{P}^{n+\gamma} = \mathbf{F}_p^n. \quad (80) $$

Here, $\mathbf{P}^{n+\gamma}$ denotes the vector of the degrees of freedom for the pressure. Moreover, we set

$$ \mathbf{M}_{p_{ij}}^{n+\gamma} = \sum_{K \in \mathcal{T}_h} \int_K \psi_j \psi_i d\Omega \quad (81) $$

$$ \mathbf{K}_{p_{ij}} = \sum_{K \in \mathcal{T}_h} \int_K \nabla \psi_j \cdot \nabla \psi_i d\Omega - \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \frac{\nabla \psi_j}{\rho^{n+\gamma}} \right\} \right\} \cdot [[\psi_i]] d\Sigma \quad (82) \\ - \sum_{e \in \mathcal{E}} \int_e [[\psi_j]] \cdot \left\{ \left\{ \frac{\nabla \psi_i}{\rho^{n+\gamma}} \right\} \right\} d\Sigma + \sum_{e \in \mathcal{E}} \int_e C_p \left\{ \left\{ \frac{1}{\rho^{n+\gamma}} \right\} \right\}_H [[\psi_j]] \cdot [[\psi_i]] d\Sigma $$

and

$$ \mathbf{F}_p^n = \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{\gamma \Delta t} \mathbf{u}^{n+\gamma,*} \cdot \nabla q d\Omega - \sum_{e \in \mathcal{E}} \int_e \frac{1}{\gamma \Delta t} \left\{ \left\{ \mathbf{u}^{n+\gamma,**} \right\} \right\} \cdot [[q]] d\Sigma. \quad (83) $$

The second TR-BDF2 stage can be described in a similar manner according to the formulations reported in (47), (50), and (52).

Finally, we consider the weak formulation for the reinitialization equation for the level set function (54):

$$ \begin{aligned} & \sum_{K \in \mathcal{T}_h} \int_K \frac{\phi^{k+1,*}}{\Delta \tau} w d\Omega - \sum_{K \in \mathcal{T}_h} \int_K u_c \phi^{k+1,*} (1 - \phi^{k,*}) \mathbf{n}_\Gamma \cdot \nabla w d\Omega \\ & + \sum_{e \in \mathcal{E}} \int_e u_c \left\{ \left\{ \phi^{k+1,*} (1 - \phi^{k,*}) \mathbf{n}_\Gamma \right\} \right\} \cdot [[w]] d\Sigma + \sum_{e \in \mathcal{E}} \int_e \frac{\tilde{\lambda}^k}{2} \left[ \left[ \phi^{k+1,*} \right] \right] \cdot [[w]] d\Sigma \\ & + \sum_{K \in \mathcal{T}_h} \int_K u_c \beta \varepsilon \left( \nabla \phi^{k+1,*} \cdot \mathbf{n}_\Gamma \right) \mathbf{n}_\Gamma \cdot \nabla w d\Omega \quad (84) \\ & - \sum_{e \in \mathcal{E}} \int_e u_c \beta \varepsilon \left\{ \left\{ (\nabla \phi^{k+1,*} \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma \right\} \right\} \cdot [[w]] d\Sigma \\ & - \sum_{e \in \mathcal{E}} \int_e u_c \beta \varepsilon \left\{ \left\{ (\nabla v \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma \right\} \right\} \cdot \left[ \left[ \phi^{k+1,*} \right] \right] d\Sigma + \sum_{e \in \mathcal{E}} \int_e C_\phi \left[ \left[ \phi^{k+1,*} \right] \right] \cdot [[w]] d\Sigma \\ & = \sum_{K \in \mathcal{T}_h} \int_K \frac{\phi^{k,*}}{\Delta \tau} w d\Omega, \end{aligned} $$

where

$$ \tilde{\lambda}^k = \max \left( \left| \left( 1 - (\phi^{k,*})^+ \right) \mathbf{n}_\Gamma^+ \cdot \mathbf{n}^+ \right|, \left| \left( 1 - (\phi^{k,*})^- \right) \mathbf{n}_\Gamma^- \cdot \mathbf{n}^- \right| \right). \quad (85) $$
Moreover, we set

$$ \sigma_{e,K}^{\phi} = (r+1)^2 \frac{\mathrm{diam}(e)}{\mathrm{diam}(K)}, \qquad (86) $$

so that

$$ C_{\phi} = \frac{1}{2} (\sigma_{e,K^{+}}^{\phi} + \sigma_{e,K^{-}}^{\phi}) \quad \text{if } e \in \mathcal{E}^{I}, \quad C_{\phi} = \sigma_{e,K}^{\phi} \quad \text{if } e \in \mathcal{E}^{B}. \qquad (87) $$

One can notice that, following [48], an upwind flux has been employed for the compression term and the SIP has been adopted for the diffusive term. Finally, the algebraic form is obtained considering $w = \Psi_i, i = 1, \dots, \dim(X_h)$ and the representation in (61) so as to obtain

$$ \left( \frac{1}{\Delta \tau} \mathbf{M}_{\phi} + u_c \mathbf{C}_{\phi} + \mathbf{A}_{\phi} \right) = \mathbf{F}_{\phi}. \qquad (88) $$

Here

$$ \mathbf{C}_{\phi_{ij}} = - \sum_{K \in \mathcal{T}_h} \int_K (1 - \phi^{k,*}) \mathbf{n}_\Gamma \psi_j \cdot \nabla \psi_i d\Omega \\ + \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ (1 - \phi^{k,*}) \mathbf{n}_\Gamma \psi_j \right\} \right\} : [[\psi_i]] d\Sigma + \sum_{e \in \mathcal{E}} \int_e \frac{\tilde{\lambda}^k}{2} [[\psi_j]] : [[\psi_i]] d\Sigma \quad (89) $$

$$ \mathbf{A}_{\phi} = \sum_{K \in \mathcal{T}_h} \int_K u_c \beta \varepsilon (\nabla \psi_j \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma \cdot \nabla \psi_i d\Omega \\ - \sum_{e \in \mathcal{E}} \int_e u_c \beta \varepsilon \{ \{ (\nabla \psi_j \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma \} \} \cdot [[\psi_i]] d\Sigma \\ - \sum_{e \in \mathcal{E}} \int_e u_c \beta \varepsilon [[\psi_j]] \cdot \{ \{ (\nabla \psi_i \cdot \mathbf{n}_\Gamma) \mathbf{n}_\Gamma \} \} d\Sigma \\ + \sum_{e \in \mathcal{E}} \int_e C_{\phi} [[\psi_j]] \cdot [[\psi_i]] d\Sigma. \qquad (90) $$

Adaptive Mesh Refinement (AMR) will be employed for the numerical tests in Section 5. Notice that, AMR does not affect the spatial discretization strategy. More specifically, for faces between cells of different refinement level, the integration is performed from the refined side and a suitable interpolation is performed on the coarse side. Hence, no hanging nodes appear in the implementation of the discrete weak form of the equations.

# 5 Numerical experiments

The numerical method outlined in the previous Sections has been validated in a number of classical test cases for incompressible two-phase flows using the numerical library deal.II [2], whose adaptive mesh refinement capabilities will be employed to enhance resolution close to the interface. The only constraint of the library in the use of non-conforming adaptive meshes is the requirement of not having neighbouring cells with refinement levels differing by more than one. We set $h = \min \{\mathrm{diam}(K) | K \in \mathcal{T}_h\}$
and we define two Courant numbers, one based on the flow velocity, denoted by $C_u$, and one based on the Mach number, denoted by $C$:

$$C_u = k \frac{\Delta t U}{h} \qquad C = k \frac{1}{M} \frac{\Delta t}{h}, \qquad (91)$$

where $U$ is the magnitude of the flow velocity. For the sake of convenience of the reader, we recall here that $k$ and $k-1$ are the polynomial degrees of the finite element spaces chosen for the discretization of velocity and pressure, respectively, whereas $r$ is the polynomial degree of the finite element space chosen for the discretization of the level set function. We consider $k = r = 2$ in all the numerical experiments, since the TR-BDF2 is a second order time discretization method, so as to obtain a second order flow solver. Finally, we employ homogenous Neumann boundary conditions for the level set function, while the initial value of the pressure is equal to zero, so as to avoid an extra initial forcing term in the momentum balance.

## 5.1 Static bubble

We consider first the 2D stationary bubble in a zero force field described e.g. in [15, 28, 50]. A bubble with $R = 0.25$ m centered in $(x_0, y_0) = (0.5, 0.5)$ m in $\Omega = (0, 1)^2$ m$^2$ is considered. Following [28], the fluid properties are $\rho_1 = \rho_2 = 10^4$ kg m$^{-3}$, $\mu_1 = \mu_2 = 1$ kg m$^{-1}$ s$^{-1}$, and $\sigma = 1$ N m$^{-1}$. Hence, $\frac{\rho_1}{\rho_2} = \frac{\mu_1}{\mu_2} = 1$. In [33, 50], it was conjectured that the amplitude of spurious currents, which arise considering the contribution of the surface tension, is proportional to $\frac{\sigma}{\mu_1}$. Hence, we take $U_{ref} = \frac{\sigma}{\mu_1} = 1$ m s$^{-1}$, so that the Reynolds number equals the Laplace number $La = \frac{\rho_1 \sigma L_{ref}}{\mu_1^2}$. The reference length is equal to $L_{ref} = 2R = 0.5$ m, so that the computational domain is $\Omega = (0, 2)^2$, while the final time is $T_f = 250$, i.e. we consider 250 characteristic times. The artificial speed of sound is set to $c \approx 1428$ m s$^{-1}$, which is of the same order of magnitude of the speed of sound in the water. Hence, we get $Re = 5 \times 10^3$, $We = 5 \times 10^3$, and $M = 7 \times 10^{-4}$. We consider no-slip boundary conditions for the velocity and homogeneous Neumann boundary conditions for the pressure. Finally, we set $\varepsilon = \frac{1}{20}$, $\Delta\tau = 6.25 \times 10^{-3}$, $u_c = 0.05u_{max}$ and $\beta = 0.5$, where $u_{max}$ is the maximum fluid velocity. The choice to relate $u_c$ with $u_{max}$ is rather common in the literature, see e.g. [11, 54]. The expected convergence rates are obtained, showing that the numerical method is globally of second order (Table 1). Spurious currents typically appear in the form of vortices around the interface. One can notice that the treatment of the surface tension contribution using the Laplace-Beltrami operator (16) reduces the generation of spurious currents, even for a coarse mesh (Figure 1).

## 5.2 Rayleigh-Taylor instability

The Rayleigh-Taylor instability is a well known test case in which a heavier fluid penetrates a lighter fluid under the action of gravity. We consider the configuration presented e.g. in [6, 25, 50], for which $\rho_1 = 1.225$ kg m$^{-3}$ and $\rho_2 = 0.1694$ kg m$^{-3}$, corresponding to the density of air and helium, respectively, whereas $\mu_1 = \mu_2 = 0.00313$ kg m$^{-1}$ s$^{-1}$. The effect of surface tension is neglected. Following [62], we consider as reference length, the computational width of the box $W$ and as reference
<table><caption>Table 1: Static bubble test case, convergence rates at fixed $\frac{U_{ref} \Delta t}{h}$ in norm $H^1$, $L^2$, and $L^\infty$. $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th>N<sub>el</sub></th><th>Δt</th><th>||u||<sub>H<sup>1</sup>(Ω)</sub></th><th>Conv. rate H<sup>1</sup></th><th>||u||<sub>L<sup>2</sup>(Ω)</sub></th><th>Conv. rate L<sup>2</sup></th><th>||u||<sub>L<sup>∞</sup>(Ω)</sub></th><th>Conv. rate L<sup>∞</sup></th></tr></thead><tbody><tr><td>20</td><td>1</td><td>4.37 × 10<sup>-2</sup></td><td></td><td>4.81 × 10<sup>-3</sup></td><td></td><td>4.77 × 10<sup>-3</sup></td><td></td></tr><tr><td>40</td><td>0.5</td><td>7.59 × 10<sup>-3</sup></td><td>2.5</td><td>7.48 × 10<sup>-5</sup></td><td>2.7</td><td>2.81 × 10<sup>-4</sup></td><td>4.1</td></tr><tr><td>80</td><td>0.25</td><td>1.14 × 10<sup>-3</sup></td><td>2.7</td><td>7.23 × 10<sup>-6</sup></td><td>3.4</td><td>3.01 × 10<sup>-5</sup></td><td>3.2</td></tr><tr><td>160</td><td>0.125</td><td>2.43 × 10<sup>-4</sup></td><td><b>2.2</b></td><td>8.64 × 10<sup>-7</sup></td><td><b>3.1</b></td><td>3.32 × 10<sup>-6</sup></td><td><b>3.2</b></td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 1: Static bubble test case, velocity magnitude at $t = T_f = 250$ using a mesh composed by 40 × 40 elements.</figcaption></figure>

time, the time scale of wave growth, equal to $t_{ref} = \sqrt{\frac{W}{Ag}}$, where $g = 9.81 \text{ m s}^{-2}$ and $A = \frac{\rho_1 - \rho_2}{\rho_1 + \rho_2}$ is the Atwood number. Hence, we obtain the following relations:

$$U_{ref} = \sqrt{AgW} \qquad Re = \frac{\rho_1 \sqrt{AgW}W}{\mu_1} \qquad Fr = \sqrt{A}. \qquad (92)$$

We consider $W = 1 \text{ m}$, so as to obtain a computational domain $\Omega = (0, 1) \times (0, 4)$. Hence, we get $A \approx 0.757$, $t_{ref} \approx 0.367 \text{ s}$, $U_{ref} \approx 2.725 \text{ m s}^{-1}$, $Re \approx 1066.55$, and $Fr \approx 0.87$. We take $M = 0.008$, corresponding to $c \approx 343 \text{ m s}^{-1}$, namely the speed of sound in air. The final time is $T_f = 2.45$. No-slip boundary conditions are prescribed on top and bottom walls, whereas periodic boundary conditions are imposed along the horizontal direction. The pressure is prescribed to be zero on the upper wall, while homogeneous Neumann boundary conditions are employed for the bottom wall. The initial velocity field is zero, whereas the initial level set function is

$$\phi(0) = \frac{1}{1 + \exp\left(\frac{2+0.05 \cos(2\pi x) - y}{\epsilon}\right)}. \qquad (93)$$
The computational grid is composed by 160 × 640 elements, whereas the time step is Δt ≈ 1.63 × 10⁻³, yielding a maximum advective Courant number C ≈ 1.36 and an acoustic Courant number C ≈ 65.3. Finally, we set ε = h = 1/160, Δτ = 0.05h, u_c = 0.0125u_max, and β = 1. Figure 2 shows the development of the interface at t = T_f, where one can easily notice the expected main behaviour of the Rayleigh-Taylor instability: as the heavier fluid penetrates the lighter one, the interface begins to roll up along the sides of the spike giving the typical “mushroom” shape. Obtained results are similar to those in literature, see e.g. [25, 50, 51, 58]. Moreover, for the sake of completeness, we report in Figure 3 the evolution of the relative variation of the area for the lighter fluid, defined as

$$ \frac{|\Omega_2(t) - \Omega_2(0)|}{\Omega_2(0)}. \qquad (94) $$

The maximum relative variation is 0.034 %, showing that CLS method preserves the area quite well. The results compare well with those in [55], in which a loss of around 1.1% was experienced.

An interesting analysis regards the influence of the Atwood number. We fix ρ₂ = 0.408 kg m⁻³, so as to obtain A ≈ 0.5. As a consequence, we obtain t ref ≈ 0.451 s, U ref ≈ 2.215 m s⁻¹, Re ≈ 867.05, Fr ≈ 0.71, and M = 0.006. We set the final time T_f = 2, so that the same final dimensional time of the previous configuration is achieved. The chosen time step is Δt = 2.5 · 10⁻³. One can easily notice from Figure 4 that, with higher Atwood number, the roll up effect is enhanced. This points out the earlier appearance of the Kelvin-Helmholtz instability, due to the development of short wavelength perturbations along the fluid interface.

We employ now the h-adaptive version of the scheme for the latter configuration. More specifically, we define for each element K the quantity

$$ \eta_K = \max_{i \in \mathcal{N}_K} |\nabla \phi|_i, \qquad (95) $$

which acts as local refinement indicator. Here N_K denotes the set of nodes over the element K. We allow to refine when η_K exceeds 10 and to coarsen below 5. The initial grid is composed by 80 × 320 elements and we allow up to two local refinements, so as to obtain h = 1/320 and a maximum resolution which would correspond to a 320 × 1280 uniform grid. As one can notice from Figure 5, the refinement criterion is able to increase the resolution only in correspondence of the interface between the two fluids. The final grid consists of 42850 elements, corresponding to around 40 % of elements of the fixed uniform grid. A solution using the full resolution uniform grid, i.e. a 320 × 1280 grid, has been also computed for further comparison. One can easily notice that at t = T_f/2 the interfaces are indistinguishable, whereas at t = T_f a slightly different development of the instability appears between the simulation with the adaptive grid and those with the two fixed grids (Figure 6). Moreover, one can notice that the same profiles of the interface are obtained for the two fixed grids, meaning that we have achieved grid independence. Since we are analyzing a fluid mechanic instability, every small variation in the flow corresponds to large variations, and, therefore, it is difficult to say which solution is the more reliable. Similar results and considerations have been reported for a Kelvin-Helmholtz instability in [43]. Finally, for what concerns the computation efficiency, AMR leads to a computational time saving of around 90 % for a given maximum spatial resolution (Table 2).
<figure><img src="image_3.png" /><figcaption>Figure 2: Rayleigh-Taylor instability, contour plot of the level set function at $t = T_f = 2.45$.</figcaption></figure>

<table><caption>Table 2: Rayleigh-Taylor instability, wall-clock times for the different configurations at $A \approx 0.5$.</caption><thead><tr><th></th><th>Grid</th><th>WT[s]</th></tr></thead><tbody><tr><td></td><td>160 × 640 (fixed)</td><td>2110</td></tr><tr><td></td><td>320 × 1280 (fixed)</td><td>14400</td></tr><tr><td>adaptive ($N_{el}$ = 42850 at $t = T_f$, maximum resolution 320 × 1280)</td><td></td><td>1830</td></tr></tbody></table>

## 5.3 Rising bubble benchmark

The rising bubble benchmark is a well-established test case for the validation of numerical methods for incompressible two-phase flows [29]. More specifically, the evolution of the shape, position and velocity of the center of mass of a rising bubble is compared against the reference solution in [29]. Two configurations are considered with
<figure><img src="image_4.png" /><figcaption>Figure 3: Rayleigh-Taylor instability, evolution of the relative variation of the area for the lighter fluid.</figcaption></figure>

the corresponding physical parameters and non-dimensional numbers listed in Table 3 and 4, respectively. The bubble occupies the subdomain $\Omega_2$. Following [29], we set $L_{ref} = 2r_0 = 0.5$ m and $U_{ref} = \sqrt{gL_{ref}} = 0.7$ m s$^{-1}$. We consider as domain $\Omega = (0, L_x) \times (0, L_y)$, with $L_x = 2$ and $L_y = 4$, whereas the final time is $T_f = 4.2$. No-slip boundary conditions are imposed on the top and bottom boundaries, with homogeneous Neumann boundary conditions for the pressure. Periodic conditions are prescribed in the horizontal direction. The initial velocity field is zero. Finally, the initial level set function is described by the following relation:

$$ \phi(0) = \frac{1}{1 + \exp\left(\frac{R - \sqrt{(x-x_0)^2 + (y-y_0)^2}}{\varepsilon}\right)}, \quad (96) $$

with $R = 1, x_0 = y_0 = 1$. We compute as reference quantities the position $\mathbf{x}_c$, the velocity $\mathbf{u}_c$ of the center of mass, and the so-called degree of circularity $\chi$, defined respectively as

$$ \mathbf{x}_c = \frac{\int_{\Omega_2} \mathbf{x} d\Omega}{\int_{\Omega_2} d\Omega} = \frac{\int_{\Omega_2} \mathbf{x} d\Omega}{|\Omega_2|} \quad (97) $$

$$ \mathbf{u}_c = \frac{\int_{\Omega_2} \mathbf{u} d\Omega}{\int_{\Omega_2} d\Omega} = \frac{\int_{\Omega_2} \mathbf{u} d\Omega}{|\Omega_2|} \quad (98) $$

$$ \chi = \frac{2\sqrt{\pi |\Omega_2|}}{P_b}, \quad (99) $$

where $\Omega_2$ is the subdomain occupied by the bubble, $|\Omega_2|$ is the area of the bubble, and $P_b$ is its perimeter. The degree of circularity is the ratio between the perimeter of a circle with the same area of the bubble and the current perimeter of the bubble itself. For a perfectly circular bubble, the degree of circularity is equal to one and then
<figure><img src="image_5.png" /><figcaption>Figure 4: Rayleigh-Taylor instability, comparison between $A \approx 0.757$ and $A \approx 0.5$. The black line shows the interface for $A \approx 0.757$, whereas the red line refers to the interface for $A \approx 0.5$.</figcaption></figure>

decreases as the bubble deforms itself. Since $\phi$ is a regularized Heaviside function, we can compute the reference quantities as follows:

$$ \mathbf{x}_c \approx \frac{\int_{\Omega} \mathbf{x}(1-\phi) \, d\Omega}{\int_{\Omega} (1-\phi) \, d\Omega} \qquad (100) $$

$$ \mathbf{u}_c \approx \frac{\int_{\Omega} \mathbf{u}(1-\phi) \, d\Omega}{\int_{\Omega} (1-\phi) \, d\Omega} \qquad (101) $$

$$ \chi \approx \frac{2\sqrt{\pi \int_{\Omega} (1-\phi) \, d\Omega}}{\int_{\Omega} |\nabla \phi| \, d\Omega}. \qquad (102) $$

We start with the first configuration and we set $M = 0.0005$, corresponding to $c = 1400 \, \text{m s}^{-1}$, which is of the order of magnitude of the speed of sound in water. The computational grid is composed by $320 \times 640$ elements, leading to $h = \frac{1}{160}$, whereas the time step is $\Delta t = 6 \cdot 10^{-3}$, yielding a maximum advective Courant number $C_u \approx 1.4$ and an acoustic Courant number $C = 1920$. Finally, we set $\epsilon = h$, $\Delta\tau = 0.05h$, $u_c = 0.05u_{max}$ and $\beta = 0.5$. We point out here the fact that results in the Figures have been compared with the results of Group 2 in [29]. One can easily notice that we
<figure><img src="image_6.png" /><figcaption>Figure 5: Rayleigh-Taylor instability at $A \approx 0.5$. Left: contour plot of the level set function at $t = T_f = 2$. Right: computational grid at $t = T_f = 2$.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Rayleigh-Taylor instability, comparison at $A \approx 0.5$ between the fixed grid simulation and the adaptive grid one. Left: results at $t = \frac{T_f}{2} = 1$. Right: results at $t = T_f = 2$. The black lines show the interface obtained with the adaptive grid, the red lines refer to the interface obtained with the fixed grid composed by 160 × 640 elements, whereas the blue lines report the results obtained with a fixed grid composed by 320 × 1280 elements.</figcaption></figure>

are able to recover the reference shape of the bubble at $t = T_f$ (Figure 7). A good qualitative agreement is established for the evolution of the degree of circularity, with only slightly lower values for our numerical results (Figure 8). A maximum relative
<table><caption>Table 3: Physical parameters defining the configurations from rising bubble test case (data from [29]).</caption><thead><tr><th>Test case</th><th>ρ<sub>1</sub> [kg m<sup>-3</sup>]</th><th>ρ<sub>2</sub> [kg m<sup>-3</sup>]</th><th>μ<sub>1</sub> [kg m<sup>-1</sup> s<sup>-1</sup>]</th><th>μ<sub>2</sub> [kg m<sup>-1</sup> s<sup>-1</sup>]</th><th>g [m s<sup>-2</sup>]</th><th>σ [kg s<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Config. 1</td><td>1000</td><td>100</td><td>10</td><td>1</td><td>0.98</td><td>24.5</td></tr><tr><td>Config. 2</td><td>1000</td><td>1</td><td>10</td><td>0.1</td><td>0.98</td><td>1.96</td></tr></tbody></table>

<table><caption>Table 4: Non-dimensional numbers defining the configurations from rising bubble test case (data from [29]).</caption><thead><tr><th>Test case</th><th>Re</th><th>Fr</th><th>We</th><th>ρ<sub>2</sub>/ρ<sub>1</sub></th><th>μ<sub>2</sub>/μ<sub>1</sub></th></tr></thead><tbody><tr><td>Config. 1</td><td>35</td><td>1</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td></tr><tr><td>Config. 2</td><td>35</td><td>1</td><td>125</td><td>10<sup>-3</sup></td><td>10<sup>-2</sup></td></tr></tbody></table>

error of the order of $10^{-1}$ is established with respect to the results of Group 2 in [29]. The center of mass reaches $y_c = 2.156$, which is in good agreement with the value $y_c = 2.162 \pm 0.002$ reported in [29] and implies a relative error of the order of $10^{-3}$ (Figure 8). Finally, the maximum rise velocity of the center of mass is $v_c = 0.3461$, which is again in good agreement with the value $v_c = 0.3456 \pm 0.0003$ present in [29] and implies again a relative error of the order of $10^{-3}$ (Figure 8).

<figure><img src="image_10.png" /><figcaption>Figure 7: Rising bubble benchmark, configuration 1, shape of bubble at $t = T_f = 4.2$. Left: numerical simulation. Right: image from [29]. Bounds have been rescaled by $L_{ref}$ for the sake of comparison with reference results.</figcaption></figure>

We analyze now the second configuration. The time step is $\Delta t = 5 \cdot 10^{-3}$, yielding a maximum advective Courant number $C_u \approx 1.4$ and an acoustic Courant number $C = 1600$. We also set $\epsilon = h = \frac{1}{160}$, $\Delta \tau = 0.05h = 3.125 \times 10^{-4}$, $u_c = 0.0125u_{max}$ and $\beta = 2$. The bubble develops a non-convex shape with thin filaments (Figure 9). The solutions given in [29] are different and, in some cases, the thin filaments tend to break off, although it is unclear if such a phenomenon should be observed in the current two-dimensional setting. The obtained profile is however in good agreement with that
<figure><img src="image_11.png" /><figcaption>Figure 8: Rising bubble benchmark, configuration 1, a) degree of circularity, b) vertical coordinate position of the center of mass, c) vertical coordinate velocity of the center of mass. The blue lines denote the results obtained with our method, while the red dots show the reference results from [29]. Reference data for position and velocity have been rescaled by $L_{ref}$ and $U_{ref}$, respectively, for the sake of comparison.</figcaption></figure>

of Group 2 in [29]. Figure 10 shows the evolution of the degree of circularity, of the vertical coordinate of the position of the center of mass, and of the vertical coordinate of the velocity of the center of mass. A good qualitative agreement is established for the quantities of interest, even though deviations from the chosen reference solution are visible. In particular, differences appear for the degree of circularity starting from $t \approx 2.5$, when the thin filaments start developing. Moreover, the second peak for the rising velocity reaches a lower value. As mentioned above, there is no clear agreement concerning the thin filamentary regions, and, therefore, their development can strongly affect computations of the reference quantities and can lead to different numerical results.

We employ now AMR to increase the resolution in correspondence of the interface. We consider the same refinement criterion (95) and the same thresholds for $\eta_K$ adopted in Section 5.2 and we allow up two local refinements, so as to obtain $h = \frac{1}{640}$ and a maximum resolution which would correspond to a 1280 × 2560 uniform grid. One can notice that the resolution is enhanced close to the interface between the two fluids and a comparison of the shape of the bubble with that obtained with the uniform grid shows that we have reached grid independence (Figure 11). The final grid consists of 274427 elements. Figure 12 reports a comparison for the quantities of interest between the fixed grid simulation, the adaptive one and the reference results. The profiles confirm that we have reached grid independence, since only the degree of circularity slightly differs between the two simulations, whereas the profiles of the vertical coordinates of both velocity and position of the center of mass are visually indistinguishable.
<figure><img src="image_12.png" /><figcaption>Figure 9: Rising bubble benchmark, configuration 2, shape of bubble at $t = T_f = 4.2$. Left: numerical simulation. Right: image from [29]. Bounds have been rescaled by $L_{ref}$ for the sake of comparison with the reference results.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Rising bubble benchmark, configuration 2, a) degree of circularity, b) vertical coordinate position of the center of mass, c) vertical coordinate velocity of the center of mass. The blue lines denote the results obtained with our method, while the red dots show the reference results from [29]. Reference data for position and velocity have been rescaled by $L_{ref}$ and $U_{ref}$, respectively, for the sake of comparison.</figcaption></figure>

A significant difference in the development of the thin filamentary regions depends on the modelling of the viscosity coefficient $\mu$, as pointed out in [52] for diffuse interface
<figure><img src="image_14.png" /><figcaption>Figure 11: Rising bubble benchmark, configuration 2 with adaptive mesh refinement. Left: shape of bubble at $t = T_f = 4.2$. The black line reports the interface obtained with the adaptive mesh, while the red line shows the interface obtained with the fixed grid. Right: computational grid at $t = T_f = 4.2$ (close-up to thin filamentary regions).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12: Rising bubble benchmark, configuration 2 with adaptive mesh refinement, a) vertical coordinate position of the center of mass, b) vertical coordinate velocity of the center of mass, c) degree of circularity. The blue lines denote the results obtained with the uniform grid, the black dashed lines report the results obtained with the adaptive grid, while the red dots show the reference results from [29]. Reference data for position and velocity have been rescaled by $L_{ref}$ and $U_{ref}$, respectively, for the sake of comparison.</figcaption></figure>
models. A popular alternative to the linear interpolation model defined in (32) is the so-called harmonic interpolation, defined as

$$ \frac{1}{\mu} = H_{\varepsilon}(\varphi) + \frac{\mu_1}{\mu_2} (1 - H_{\varepsilon}(\varphi)). \qquad (103) $$

This choice yields results which are more similar to Group 1 in [29], where a break-up occurs (Figure 13). For what concerns the quantities of interest, we notice from Figure 14 that, since the thin elongated filaments break themselves, the degree of circularity is higher. Moreover, both the second peak of the rising velocity and the final position of the center of mass are significantly higher and closer to the reference results. The following analysis further confirms how challenging is defining a reference benchmark solution when the bubble undergoes large deformations.

<figure><img src="image_16.png" /><figcaption>Figure 13: Rising bubble benchmark, configuration 2, shape of bubble at $t = T_f = 4.2$. Left: numerical simulation using the harmonic interpolation (103) for the viscosity. Right: image from [29]. Bounds have been rescaled by $L_{ref}$ for the sake of comparison with the reference results.</figcaption></figure>

## 5.4 Falling bubble test case

Finally, we consider the test case of a falling bubble in order to further analyze the impact of the choice of the mixture viscosity. The present configuration is a modification of that presented in [23]. The domain is $\Omega = (0, 1) \times (0, \frac{3}{2})$ m$^2$. The bubble is initially an ellipse centered in $(x_0, y_0) = (0.6, 1.2)$ m with semi-axes $a = 0.2$ m and $b = \frac{\sqrt{2}}{10}$ m. The density of the heavier fluid is $\rho_1 = 1000$ kg m$^{-3}$, while its viscosity is $\mu_1 = 4$ kg m$^{-1}$ s$^{-1}$. Following [23], we take $L_{ref} = a = 0.2$ m, $U_{ref} = 0.5$ m s$^{-1}$, $\rho_1 = 10^{-2}$, and $\mu_1 = \frac{1}{8}$. Hence, the computational domain is $\Omega = (0, 5) \times (0, 7.5)$, with $Re = 25$ and $Fr \approx 0.35$. We take $M = 3.3 \times 10^{-4}$, corresponding to $c \approx 1500$ m s$^{-1}$. The computational grid is composed by 160 $\times$ 240 elements and we set $\varepsilon = h = 0.03125$, $\Delta\tau = 0.05h$, $u_c = 0.05u_{max}$, and $\beta = 0.5$. Surface tension is not included. Hence, the bubble will tend to break itself and we aim to show the different evolution
<figure><img src="image_17.png" /><figcaption>Figure 14: Rising bubble benchmark, configuration 2, comparison between linear interpolation and harmonic interpolation for the viscosity: a) vertical coordinate of the position of the center of mass, b) vertical coordinate of the velocity of the center of mass, c) degree of circularity. The blue lines denote the results with the linear interpolation, the black lines report the results with the harmonic interpolation, while the red dots show the reference results from [29]. Reference data for position and velocity have been rescaled by $L_{ref}$ and $U_{ref}$, respectively, for the sake of comparison.</figcaption></figure>

of the interface according to the choice of the mixture viscosity. No-slip boundary conditions are imposed on top and bottom boundaries, together with homogeneous Neumann boundary conditions for the pressure. Periodic boundary conditions are prescribed in the horizontal direction. The initial velocity field is zero, while the initial level set function is described by the following relation:

$$ \phi(0) = \frac{1}{1 + \exp\left(\frac{\sqrt{\frac{b^2}{a^2}\left(x - \frac{x_0}{L_{ref}}\right)^2 + \left(y - \frac{y_0}{L_{ref}}\right)^2 - \frac{b}{L_{ref}}}}{\varepsilon}\right)}. \quad (104) $$

We take as final time $T_{fin} = 0.75$, with $\Delta t = 2.5 \times 10^{-3}$. yielding an acoustic Courant number $C \approx 485$. One can easily notice that, at $t = \frac{3}{5}T_f = 0.45$, the interface obtained with the harmonic interpolation of the viscosity (103) and that obtained using the linear interpolation of the viscosity (32) are visually indistinguishable and exhibit some oscillations (Figure 15, left). Finally, at $t = T_f$, the two interfaces are significantly different (Figure 15, right). The interface obtained using (103) shows evident oscillations and numerous small scale structures, while the interface tracked employing (32) is deforming itself with reduced oscillations and less small scale structures. The presence of irregular values using the harmonic interpolation of the viscosity is confirmed by the maximum Courant number $C_u$, which is $C_u \approx 3.6$ using (103) and $C_u \approx 1$ using
(32). These results further validate the considerations depicted in Section 5.3: the use of the harmonic interpolation of the viscosity (103) tends to enhance the break-up effect and the presence of small scale structures, while the linear interpolation of the viscosity (32) is able to keep a more regular shape of the interface for longer time. Hence, the modelling of the mixture viscosity plays a key role for the description of physical phenomena in which atomization and break-up occur, like those presented e.g. in [34, 63], or when small scale structures are explicitly considered as in [35].

<figure><img src="image_18.png" /><figcaption>Figure 15: Falling bubble test case with a 160 × 240 mesh. Left: shapes of the bubble at $t = \frac{3}{5}T_f = 0.45$. Right: shapes of the bubble at $t = T_f = 0.75$. The black lines report the solutions obtained using the linear interpolation for the viscosity (32), while the red lines show the results obtained using the harmonic interpolation for the viscosity (103).</figcaption></figure>

# 6 Conclusions

Building on the experience of [40, 45], we have proposed an implicit Discontinuous Galerkin discretization for incompressible two-phase flows. While discretizations of incompressible two-phase flows equations have been proposed in many other papers, we have presented here an approach based on an artificial compressibility formulation in order to avoid some well known issues of projection methods. The time discretization is obtained by a projection method based on the L-stable TR-BDF2 method. The implementation has been carried out in the framework of the numerical library deal.II, whose mesh adaptation capabilities have been exploited to increase the resolution in correspondence of the interface between the two fluids. The effectiveness of the proposed approach has been shown in a number of classical benchmarks. In particular, the influence of some possible choices for the mixture viscosity when the interface undergoes large deformations or when break-up occurs has been established, following and extending an analysis previously carried out for diffuse interface models.
In future work, we aim to exploit the possibility of considering well resolved interfaces for an analysis on the evolution equations of interfacial quantities, following the recent contributions [41, 42], as well as to develop an extension of analogous numerical approaches to fully compressible two-phase flows.

## Acknowledgements

G.O. is part of the INdAM-GNCS National Research Group. The author would like to thank L. Bonaventura, P. Barbante, and E. Capuano for several useful discussions on related topics. The author gratefully acknowledges N. Parolini for providing the original data of the rising bubble test case discussed in Section 5. The simulations have been partly run at CINECA thanks to the computational resources made available through the NUMNETF-HP10C06Y02 ISCRA-C project. We acknowledge the CINECA award, for the availability of high-performance computing resources and support. This work has been partially supported by the ESCAPE-2 project, European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 800897).

## References

[1] P.F. Antonietti, I. Mazzieri, A. Quarteroni, and F. Rapetti. Non-conforming high order approximations of the elastodynamics equation. *Computer Methods in Applied Mechanics and Engineering*, 209:212–238, 2012.

[2] D. Arndt, W. Bangerth, M. Feder, M. Fehling, R. Gassmöller, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J-P. Pelteret, S. Sticko, B. Turcksin, and D. Wells. The deal.II Library, Version 9.4. *Journal of Numerical Mathematics*, 0, 2022.

[3] D.N. Arnold. An interior penalty finite element method with discontinuous elements. *SIAM Journal of Numerical Analysis*, 19:742–760, 1982.

[4] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. *SIAM journal on numerical analysis*, 39:1749–1779, 2002.

[5] R.E. Bank, W.M. Coughran, W. Fichtner, E.H. Grosse, D.J. Rose, and R.K. Smith. Transient simulation of silicon devices and circuits. *IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems*, 4:436–451, 1985.

[6] C. Bassi, A. Abbà, L. Bonaventura, and L. Valdettaro. Large eddy simulation of non-Boussinesq gravity currents with a DG method. *Theoretical and Computational Fluid Dynamics*, 34:231–247, 2020.

[7] F. Bassi, L.A. Botti, A. Colombo, and F.C. Massa. Assessment of an Implicit Discontinuous Galerkin Solver for Incompressible Flow Problems with Variable Density. *Applied Sciences*, 12:11229, 2022.

[8] F. Bassi, A. Crivellini, D.A. Di Pietro, and S. Rebay. An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier-Stokes equations. *Journal of Computational Physics*, 218:794–815, 2006.
[9] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. *Journal of Computational Physics*, 131:267-279, 1997.

[10] L. Bonaventura and A. Della Rocca. Unconditionally strong stability preserving extensions of the TR-BDF2 method. *Journal of Scientific Computing*, 70:859-895, 2017.

[11] P.H. Chiu and Y.T. Lin. A conservative phase field method for solving incompressible two-phase flows. *Journal of Computational Physics*, 230:185-204, 2011.

[12] A.J. Chorin. A numerical method for solving incompressible viscous flow problems. *Journal of Computational Physics*, 2:12-26, 1967.

[13] A.J. Chorin. Numerical solution of the Navier-Stokes equations. *Mathematics of computation*, 22:745-762, 1968.

[14] A. Della Rocca. *Large-Eddy Simulations of Turbulent Reacting Flows with Industrial Applications*. PhD thesis, Politecnico di Milano, 2 2018.

[15] S.S. Deshpande, L. Anumolu, and M.F. Trujillo. Evaluating the performance of the two-phase flow solver interFoam. *Computational science & discovery*, 1:014016, 2012.

[16] W. Dettmer, P.H. Saksono, and D. Perić. On a finite element formulation for incompressible Newtonian fluid flows on moving domains in the presence of surface tension. *Communications in Numerical Methods in Engineering*, 19:659-668, 2003.

[17] D.A. Di Pietro, S. Lo Forte, and N. Parolini. Mass preserving finite element implementations of the level set method. *Applied Numerical Mathematics*, 56:1179-1195, 2006.

[18] R. Estrada and R.P. Kanwal. Applications of distributional derivatives to wave propagation. *IMA Journal of Applied Mathematics*, 26:39-63, 1980.

[19] R. Fedkin and S. Osher. *The Level Set Methods and Dynamic Implicit Surfaces*. Springer Nature, 2003.

[20] N. Fehn, M. Kronbichler, C. Lehrenfeld, G. Lube, and P. Schroeder. High-order DG solvers for under-resolved turbulent incompressible flows: A comparison of $L^2$ and $H(\text{div})$ methods. *International Journal of Numerical Methods in Fluids*, 91:533-556, 2019.

[21] N. Fehn, W.A. Wall, and M. Kronbichler. On the stability of projection methods for the incompressible Navier-Stokes equations based on high-order discontinuous Galerkin discretizations. *Journal of Computational Physics*, 351:392-421, 2017.

[22] P. Gao, J. Ouyang, and W. Zhou. Development of a finite element/discontinuous Galerkin/level set approach for the simulation of incompressible two phase flow. *Advances in Engineering Software*, 118:45-59, 2018.

[23] J. Grooss and J.S. Hesthaven. A level set discontinuous Galerkin method for free surface flows. *Computer Methods in Applied Mechanics and Engineering*, 195:3406-3429, 2006.

[24] J.L. Guermond, P. Minev, and J. Shen. An overview of projection methods for incompressible flows. *Computer Methods in Applied Mechanics and Engineering*, 195:6011-6045, 2006.
[25] M. Haghshenas, J.A. Wilson, and R. Kumar. Algebraic coupled level set-volume of fluid method for surface tension dominant two-phase flows. *International Journal of Multiphase Flow*, 90:13–28, 2017.

[26] M.E. Hosea and L.F. Shampine. Analysis and implementation of TR-BDF2. *Applied Numerical Mathematics*, 20:21–37, 1996.

[27] T.J.R. Hughes, W.K. Liu, and T.K. Zimmermann. Lagrangian-Eulerian finite element formulation for incompressible viscous flows. *Computer Methods in Applied Mechanics and Engineering*, 29:329–349, 1981.

[28] S. Hysing. A new implicit surface tension implementation for interfacial flows. *International Journal for Numerical Methods in Fluids*, 51:659–672, 2006.

[29] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska. Quantitative benchmark computations of two-dimensional bubble dynamics. *International Journal for Numerical Methods in Fluids*, 60:1259–1288, 2009.

[30] G. Karniadakis, M. Israeli, and S.A. Orszag. High-order splitting methods for the incompressible Navier-Stokes equations. *Journal of computational physics*, 97(2):414–443, 1991.

[31] G. Karniadakis and S. Sherwin. *Spectral/hp element methods for computational fluid dynamics*. Oxford University Press on Demand, 2005.

[32] M. Kronbichler, A. Diagne, and H. Holmgren. A fast massively parallel two-phase flow solver for microfluidic chip simulation. *The International Journal of High Performance Computing Applications*, 32:266–287, 2018.

[33] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti. Modelling merging and fragmentation in multiphase flows with surfer. *Journal of Computational Physics*, 113:134–147, 1994.

[34] R. Lebas, T. Menard, P.-A. Beau, A. Berlemont, and F.-X. Demoulin. Numerical simulation of primary break-up and atomization: Dns and modelling study. *International Journal of Multiphase Flow*, 35:247–260, 2009.

[35] A. Loison, T. Pichard, S. Kokh, and M. Massot. Two-scale modelling of two-phase flows based on the Stationary Action Principle and a Geometric Method Of Moments, 2023.

[36] J. Manzanero, G. Rubio, D.A. Kopriva, E. Ferrer, and E. Valero. An entropy-stable discontinuous Galerkin approximation for the incompressible Navier-Stokes equations with variable density and artificial compressibility. *Journal of Computational Physics*, 408:109241, 2020.

[37] S. Mirjalili, S.S. Jain, and M. Dodd. Interface-capturing methods for two-phase flows: An overview and recent developments. *Center for Turbulence Research Annual Research Briefs*, 2017:13, 2017.

[38] E. Olsson and G. Kreiss. A conservative level set method for two phase flow. *Journal of Computational Physics*, 210:225–246, 2005.

[39] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii. *Journal of Computational Physics*, 225:785–807, 2007.
[40] G. Orlando. Modelling and simulations of two-phase flows including geometric variables. PhD thesis, Politecnico di Milano, 2023. http://hdl.handle.net/10589/198599.

[41] G. Orlando, P.F. Barbante, and L. Bonaventura. On the evolution equations of interfacial variables in two-phase flows. *International Journal of Multiphase Flow*, 169:104617, 2023.

[42] G. Orlando, P.F. Barbante, and L. Bonaventura. A comparison of different approaches to compute surface tension contribution in incompressible two-phase flows, 2024.

[43] G. Orlando, T. Benacchio, and L. Bonaventura. An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement. *Journal of Computational and Applied Mathematics*, 427:115124, 2023.

[44] G. Orlando and L. Bonaventura. An asymptotic-preserving scheme for euler equations i: non-ideal gases, 2024.

[45] G. Orlando, A. Della Rocca, P.F. Barbante, L. Bonaventura, and N. Parolini. An efficient and accurate implicit DG solver for the incompressible Navier-Stokes equations. *International Journal for Numerical Methods in Fluids*, 94:1484–1516, 2022.

[46] S. Osher and R.P. Fedkiw. Level set methods and dynamic implicit surfaces, volume 1. Springer New York, 2005.

[47] S. Osher and J.A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. *Journal of Computational Physics*, 79:12–49, 1988.

[48] M. Owkes and O. Desjardins. A discontinuous Galerkin conservative level set scheme for interface capturing in multiphase flows. *Journal of Computational Physics*, 249:275–302, 2013.

[49] F. Pochet, K. Hillewaert, P. Geuzaine, J.F. Remacle, and E. Marchandise. A 3d strongly coupled implicit discontinuous Galerkin level set-based method for modeling two-phase flows. *Computers & Fluids*, 87:144–155, 2013.

[50] S. Popinet and S. Zaleski. A front-tracking algorithm for accurate representation of surface tension. *International Journal for Numerical Methods in Fluids*, 30:775–793, 1999.

[51] E.G. Puckett, A.S. Almgren, J.B. Bell, D.L. Marcus, and W.J. Rider. A high-order projection method for tracking fluid interfaces in variable density incompressible flows. *Journal of computational physics*, 130(2):269–282, 1997.

[52] M. Řehoř. Diffuse interface models in theory of interacting continua. PhD thesis, 2018.

[53] M. Rudman. Volume-tracking methods for interfacial flow calculations. *International Journal for Numerical Methods in Fluids*, 24:671–691, 1997.

[54] H. Rusche. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. PhD thesis, Imperial College London (University of London), 2003.
[55] A. Salih and S. Ghosh Moulic. Simulation of rayleigh-taylor instability using level set method. 12 2006.

[56] S. Salsa. *Partial differential equations in action: from modelling to theory*, volume 99. Springer, 2016.

[57] D. Schötzau, C. Schwab, and A. Toselli. Stabilized DGFEM for incompressible flows. *Mathematical Models and Methods in Applied Sciences*, 13:1413–1436, 2003.

[58] T.W.H. Sheu, C.H. Yu, and P.H. Chiu. Development of a dispersively accurate conservative level set scheme for capturing interface in two-phase flows. *Journal of Computational Physics*, 228:661–686, 2009.

[59] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressible two-phase flow. *Journal of Computational Physics*, 114:146–159, 1994.

[60] R. Temam. Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (ii). *Archive for rational mechanics and analysis*, 33:377–385, 1969.

[61] A. Toselli. *hp Discontinuous Galerkin approximations for the Stokes problem*. *Mathematical Models and Methods in Applied Sciences*, 12:1565–1597, 2002.

[62] G. Tryggvason. Numerical simulations of the Rayleigh-Taylor instability. *Journal of Computational Physics*, 75:253–282, 1988.

[63] A. Vallet and R. Borghi. Modélisation Eulerienne de l'atomisation d'un jet liquide. *Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics-Physics-Astronomy*, 327:1015–1020, 1999.

[64] S. Zahedi. *Numerical methods for fluid interface problems*. PhD thesis, KTH Royal Institute of Technology, 2011.