# A comparison of different approaches to compute surface tension contribution in incompressible two-phase flows

Giuseppe Orlando, Paolo Francesco Barbante, Luca Bonaventura

► To cite this version:

Giuseppe Orlando, Paolo Francesco Barbante, Luca Bonaventura. A comparison of different approaches to compute surface tension contribution in incompressible two-phase flows. 2025. hal-04445636v4

HAL Id: hal-04445636

https://hal.science/hal-04445636v4

Preprint submitted on 26 Aug 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
A comparison of different approaches to compute surface tension contribution in incompressible two-phase flows

Giuseppe Orlando$^{(1)}$

Paolo Francesco Barbante$^{(2)}$, Luca Bonaventura$^{(2)}$

(1) CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris
Route de Saclay, 91120 Palaiseau, France
giuseppe.orlando@polytechnique.edu

(2) MOX - Dipartimento di Matematica, Politecnico di Milano
Piazza Leonardo da Vinci 32, 20133 Milano, Italy
paolo.barbante@polimi.it, luca.bonaventura@polimi.it

**Keywords:** Navier-Stokes equations, Incompressible flows, Two-phase flows, Level set, Surface tension, Curvature.
Abstract

We perform a quantitative assessment of different strategies to compute the con-
tribution due to surface tension in incompressible two-phase flows using a conservative
level set (CLS) method. More specifically, we compare classical approaches, such as
the direct computation of the curvature from the level set or the Laplace-Beltrami
operator, with an evolution equation for the mean curvature recently proposed in lit-
erature. We consider the test case of a static bubble, for which an exact solution for
the pressure jump across the interface is available, and the test case of an oscillating
bubble, showing pros and cons of the different approaches.
# 1 Introduction

Interfacial flows with surface tension play an important role in several industrial and engineering applications [15, 17]. Many modelling approaches have been proposed to capture the motion of the interface. We consider here the conservative level set (CLS) method, originally proposed in [20], [21], to which we refer for a detailed description of the scheme. In this framework, the normal to the interface and the curvature are implicitly determined from the level set function. We compare here different approaches to compute the force due to surface tension. More specifically, we consider three possible strategies: the use of the Laplace-Beltrami operator, the estimation of the total curvature directly from the level set function, and the use of an evolution equation for the mean curvature recently proposed in [23].

The paper is structured as follows: in Section 2, we briefly recall the different formulations chosen to model the surface tension force. In Section 3, we briefly outline the numerical method employed for the analysis. Section 4 is devoted to a quantitative assessment of relations introduced in Section 2 for the test case of a static bubble and for the test case of an oscillating bubble. Finally, some conclusions and insights for future work are presented in Section 5.

# 2 Mathematical model

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 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 without gravity, which read as follows [22]:

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

for $\mathbf{x} \in \Omega$, $t \in (0, T_f]$, supplied with suitable initial and boundary conditions. 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 density and viscosity are defined as

$$ \rho(\mathbf{x}, t) = \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}, t) = \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{D}(\mathbf{u})$ denotes the symmetric part of the gradient. Finally, $\mathbf{f}_{\sigma}$ represents a volumetric force which takes into account surface tension, defined as [14]

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

where $\sigma$ is the constant surface tension coefficient, $\mathbf{n}_{\Gamma}$ is the outward unit normal to $\Gamma$, $\kappa = -\nabla \cdot \mathbf{n}_{\Gamma}$ is the total curvature, and $\delta(\Gamma(t))$ is the Dirac delta distribution supported on the interface. In the following, for the sake of simplicity in the notation,
we omit the explicit dependence on space and time for the different quantities. As discussed in [16], we can rewrite the volumetric force as

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

which corresponds to the application of the Laplace-Beltrami operator. The CLS method [20, 21] describes implicitly the interface in terms of a regularized Heaviside function $\phi$ and its evolution equation reads as follows:

$$ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0. \quad (5) $$

For the sake of completeness, we report the definition of $\phi$ in terms of the signed distance function $\varphi$, employed in the classical level set method [25]. The following relation is assumed:

$$ \phi = \frac{1}{1 + e^{-\varphi/\epsilon}}, \quad (6) $$

where $\epsilon$ helps smoothing the transition of the discontinuous physical properties between the two subdomains and it is also known as interface thickness. From definition (6), it follows that

$$ \Gamma = \{ \mathbf{x} \in \bar{\Omega} : \phi = \frac{1}{2} \}. \quad (7) $$

The Continuum Surface Force (CSF) approach, introduced in [6], is employed to treat density, viscosity, and surface tension. More specifically, we set

$$ \rho \approx \rho_2 + (\rho_1 - \rho_2) \phi \quad (8) $$

$$ \mu \approx \mu_2 + (\mu_1 - \mu_2) \phi. \quad (9) $$

For the surface tension term, we recall that [9]

$$ \delta(\Gamma) = \delta \left( \phi - \frac{1}{2} \right) |\nabla \phi|, \quad (10) $$

where $\delta(\phi - \frac{1}{2})$ is the Dirac delta distribution with support equal to the interface implicitly described by $\phi = \frac{1}{2}$. Hence, we consider the following approximations:

$$ \mathbf{f}_{\sigma} \approx \nabla \cdot \left[ \sigma \left( |\nabla \phi| - \frac{\nabla \phi \otimes \nabla \phi}{|\nabla \phi|} \right) \right] \quad (11) $$

$$ \mathbf{f}_{\sigma} \approx -\sigma \nabla \cdot \left( \frac{\nabla \phi}{|\nabla \phi|} \right) \nabla \phi = \sigma \frac{\nabla \phi \cdot \mathbf{H}_{\phi} \nabla \phi - |\nabla \phi|^2 \Delta \phi}{|\nabla \phi|^3} \quad (12) $$

where we exploit the relation $\mathbf{n}_{\Gamma} = \frac{\nabla \phi}{|\nabla \phi|}$. Notice that, since $\phi$ represents a regularized Heaviside function, we approximate

$$ \delta \left( \phi - \frac{1}{2} \right) = \frac{d\theta \left( \phi - \frac{1}{2} \right)}{d\phi} \approx \frac{d\theta \left( \theta(\varphi) - \frac{1}{2} \right)}{d\phi} = \frac{d\theta(\varphi)}{d\phi} \approx \frac{d\phi}{d\phi} = 1 \quad (13) $$

with $\theta$ denoting the Heaviside function. The second relation in (12) is based on the so-called Bonnet's formula [11] for the curvature and $\mathbf{H}_{\phi}$ denotes the Hessian matrix
of $\phi$. Finally, we propose here another strategy to evaluate $\mathbf{f}_{\sigma}$. We compute $\kappa$ from an evolution equation for the mean curvature $H = \frac{\kappa}{2}$ recently proposed in [23]:

$$ \begin{aligned} \frac{\partial H}{\partial t} + \mathbf{u} \cdot \nabla H = & H (\nabla \mathbf{u}) \mathbf{n}_{\Gamma} \cdot \mathbf{n}_{\Gamma} + \frac{1}{2} \nabla \mathbf{n}_{\Gamma} : (\nabla \mathbf{u})^T \\ & - \frac{1}{2} (\nabla \mathbf{u})^T \mathbf{n}_{\Gamma} \cdot (\nabla \mathbf{n}_{\Gamma}) \mathbf{n}_{\Gamma} \\ & - \frac{1}{2} (\mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma} - \mathbf{I}) : [\nabla [(\nabla \mathbf{u})^T \mathbf{n}_{\Gamma}]]^T, \end{aligned} \quad (14) $$

or, equivalently,

$$ \frac{\partial H}{\partial t} + \mathbf{u} \cdot \nabla H = \frac{1}{2} \nabla \mathbf{n}_{\Gamma} : (\nabla \mathbf{u})^T + \frac{1}{2} \nabla \cdot [(\mathbf{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}) (\nabla \mathbf{u})^T \mathbf{n}_{\Gamma}], \quad (15) $$

so that

$$ \mathbf{f}_{\sigma} \approx 2\sigma H \nabla \phi. \quad (16) $$

The discretization of incompressible Navier-Stokes equations poses several major computational issues. In particular, the velocity $\mathbf{u}$ and the pressure $p$ are coupled by the incompressibility constraint $\nabla \cdot \mathbf{u} = 0$. We adopt here the so-called artificial compressibility formulation, originally introduced in [7] and employed in [3], [4], [18], [22], [24] among many others. The incompressibility constraint is relaxed and a time evolution equation for the pressure is introduced. Hence, the final form of the system 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})] + \mathbf{f}_{\sigma} \\ \frac{1}{\rho_0 c^2} \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} &= 0 \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &= 0, \end{aligned} \quad (17) $$

with $c$ being the artificial speed of sound and $\rho_0$ being a reference density. A dimensional analysis can be carried out (we refer to [22] for all the details), so as to obtain 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})] + \frac{1}{We} \mathbf{f}_{\sigma} \\ M^2 \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} &= 0 \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &= 0, \end{aligned} \quad (18) $$

where, with a slight abuse of notation, we employ the same symbols to mark non-dimensional quantities. Here $Re, We$ and $M$ are the Reynolds, Weber, and Mach number, respectively, defined as

$$ Re = \frac{\rho_{ref} U_{ref} L_{ref}}{\mu_{ref}} \qquad We = \frac{\rho_{ref} U_{ref}^2 L_{ref}}{\sigma} \qquad M = \frac{U_{ref}}{c}, \quad (19) $$

with $U_{ref}$ being the reference velocity, $L_{ref}$ denoting the reference length, $\rho_{ref}$ being the reference density and $\mu_{ref}$ denoting the reference viscosity.
## 3 Numerical method

In this Section, we briefly outline the numerical method employed for the discretization of system (18). We refer to [22] for a detailed description of the numerical scheme. We consider a decomposition of the domain $\Omega$ into a family of hexahedra (quadrilaterals in the two-dimensional case) $\mathcal{T}_h$ 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. Classical jump and average operators are then 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. \qquad (20) $$

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. \qquad (21) $$

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 (22) $$

$$ \{\{\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 (23) $$

We now introduce the following finite element spaces:

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

and

$$ \mathbf{Q}_k = [Q_k]^d, \qquad (25) $$

where $\mathbb{Q}_k$ is the space of polynomials of degree $k$ in each coordinate direction. The finite element spaces that will be used for the discretization of velocity and pressure are $\mathbf{V}_h = \mathbf{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 and the curvature, we consider instead $X_h = Q_r$ with $r \ge 2$, so that its gradient is at least a piecewise linear polynomial. In the present work, we consider $k = r = 2$.

We briefly recall for the convenience of the reader the formulation of the TR-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})$. 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 (26) $$

$$ \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). \qquad (27) $$

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}$ [12].
A reinitialization procedure is adopted for the level set function. We employ the following PDE [20, 21]:

$$ \frac{\partial \phi}{\partial \tau} + \nabla \cdot (u_c \phi (1 - \phi) \mathbf{n}_{\Gamma}) = \nabla \cdot (\beta \epsilon u_c (\nabla \phi \cdot \mathbf{n}_{\Gamma}) \mathbf{n}_{\Gamma}), \quad (28) $$

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 it is computed using the initial value of the level set function. The supplementary equation for the mean curvature (14)-(15) has to be discretized using (16) to evaluate the surface tension contribution. For the sake of completeness, we report the first stage of the TR-BDF2 scheme for (14)-(15):

$$ \begin{aligned} & \frac{H^{n+\gamma} - H^n}{\gamma \Delta t} + \frac{1}{2} \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla H^{n+\gamma} + \frac{1}{2} \mathbf{u}^{n+\frac{\gamma}{2}} \cdot \nabla H^n = \\ & \frac{1}{2} \left( \frac{1}{2} \nabla \mathbf{n}_{\Gamma}^{n+\gamma} : \nabla \mathbf{u}^{n+\frac{\gamma}{2}} \right) + \frac{1}{2} \left( \frac{1}{2} \nabla \mathbf{n}_{\Gamma}^n : \nabla \mathbf{u}^n \right) \tag{29} \\ & + \frac{1}{2} \left( \frac{1}{2} \nabla \cdot \left[ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^{n+\gamma} \otimes \mathbf{n}_{\Gamma}^{n+\gamma} \right) \left( \nabla \mathbf{u}^{n+\frac{\gamma}{2}} \right)^T \mathbf{n}_{\Gamma}^{n+\gamma} \right] \right) \\ & + \frac{1}{2} \left( \frac{1}{2} \nabla \cdot \left[ \left( \mathbf{I} - \mathbf{n}_{\Gamma}^n \otimes \mathbf{n}_{\Gamma}^n \right) \left( \nabla \mathbf{u}^n \right)^T \mathbf{n}_{\Gamma}^n \right] \right), \end{aligned} $$

where $\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}$ is defined by extrapolation. Hence, the
weak formulation for (29) reads as follows:

$$
\begin{align}
& \sum_{K \in \mathcal{T}_h} \int_K \frac{H^{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 H^{n+\gamma} w d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ H^{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 [[H^{n+\gamma} w]] d\Sigma \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[H^{n+\gamma}]] \cdot [[w]] d\Sigma \tag{30} \\
&= \sum_{K \in \mathcal{T}_h} \int_K \frac{H^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 H^n w d\Omega \nonumber \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ H^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 [[H^n w]] d\Sigma \nonumber \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \frac{\lambda^{n+\frac{\gamma}{2}}}{2} [[H^n]] \cdot [[w]] d\Sigma \nonumber \\
& + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{2} \nabla \mathbf{n}_\Gamma^{n+\gamma} : \nabla \mathbf{u}^{n+\frac{\gamma}{2}} w d\Omega + \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{2} \nabla \mathbf{n}_\Gamma^n : \nabla \mathbf{u}^n w d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \nabla \mathbf{n}_\Gamma^{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\{ \nabla \mathbf{n}_\Gamma^n \mathbf{u}^n \right\} \right\} \cdot [[w]] d\Sigma \nonumber \\
& - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \nabla \mathbf{n}_\Gamma^{n+\gamma} \mathbf{u}^{n+\frac{\gamma}{2}} \right\} \right\} : \langle \langle \mathbf{u}^{n+\frac{\gamma}{2}} w \rangle \rangle d\Sigma - \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \nabla \mathbf{n}_\Gamma^n \right\} \right\} : \langle \langle w \mathbf{u}^n \rangle \rangle d\Sigma \nonumber \\
& - \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{2} \left( \mathbf{I} - \mathbf{n}_\Gamma^{n+\gamma} \otimes \mathbf{n}_\Gamma^{n+\gamma} \right) \left( \nabla \mathbf{u}^{n+\frac{\gamma}{2}} \right)^T \mathbf{n}_\Gamma^{n+\gamma} \cdot \nabla w d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_\Gamma^{n+\gamma} \otimes \mathbf{n}_\Gamma^{n+\gamma} \right) \left( \nabla \mathbf{u}^{n+\frac{\gamma}{2}} \right)^T \mathbf{n}_\Gamma^{n+\gamma} \right\} \right\} \cdot [[w]] \nonumber \\
& - \frac{1}{2} \sum_{K \in \mathcal{T}_h} \int_K \frac{1}{2} \left( \mathbf{I} - \mathbf{n}_\Gamma^n \otimes \mathbf{n}_\Gamma^n \right) \left( \nabla \mathbf{u}^n \right)^T \mathbf{n}_\Gamma^n \cdot \nabla w d\Omega \nonumber \\
& + \frac{1}{2} \sum_{e \in \mathcal{E}} \int_e \left\{ \left\{ \left( \mathbf{I} - \mathbf{n}_\Gamma^n \otimes \mathbf{n}_\Gamma^n \right) \left( \nabla \mathbf{u}^n \right)^T \mathbf{n}_\Gamma^n \right\} \right\} \cdot [[w]] \quad \forall w \in X_h, \nonumber
\end{align}
$$

with

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

The numerical approximation of non-conservative terms is based on the approach pro-
posed in [5]. We recast the non-conservative term into the sum of two contributions:
the first one takes into account the elementwise gradient, whereas the second one con-
siders its jumps across the element faces. The second TR-BDF2 stage can be described
in a similar manner. Finally, we consider the following spatial discretization for **f**<sub>σ</sub>:

$$
\mathbf{f}_{\sigma} \approx 2H\nabla\phi \approx \sum_{K \in \mathcal{T}_h} \int_K 2H\nabla\phi \cdot \varphi d\Omega + \sum_{e \in \mathcal{E}} \int_e 2\{\{H\phi\}\} [[\varphi]] d\Sigma - \sum_{e \in \mathcal{E}} \int_e 2\{\{H\}\} [[\phi\varphi]] d\Sigma, \tag{32}
$$
which is again based on the approach presented in [5]. Here, $\varphi$ is a basis function of the finite element space chosen to discretize the momentum equation, i.e. $V_h$.

# 4 Numerical results

In this Section, we present the results from simulations comparing the different approaches presented in Section 2 to evaluate $f_\sigma$. The numerical method outlined in Section 3 has been implemented in the framework of the deal.II library [1, 2].

## 4.1 Static bubble

We consider the 2D stationary bubble in a zero force field described e.g. in [8, 13, 26]. According to the Laplace-Young law [10], the pressure jump across the interface of two immiscible fluids is directly related to surface tension. Indeed, the following relation holds:

$$ \Delta p = p_{in} - p_{out} = \frac{\sigma}{R}, \qquad (33) $$

where $p_{in}$ and $p_{out}$ are the pressure inside and outside the bubble, respectively, whereas $R$ is the radius of the bubble. A bubble with $R = 0.25$ m centered in $(x_0, y_0) = (0.5, 0.5)$ in $\Omega = (0, 1)^2$ is considered. 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}$, so that $\Delta p = 4$ N m$^{-2}$. Finally, we set $T_f = 25$ s and a fixed time step $\Delta t = 0.2$ s. 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 water. The reference length is equal to $L_{ref} = 2R = 0.5$ m and the reference velocity $U_{ref}$ is chosen in such a way that the reference pressure $p_{ref} = \rho_{ref}U_{ref}^2$ is unitary. Hence, we get $U_{ref} = 0.01$ m s$^{-1}$, so as to obtain $Re = 50, We = 0.5$, and $M = 7 \times 10^{-6}$. Periodic boundary conditions are considered. Following [8], we consider three metrics to assess the pressure jump computation:

1. $\Delta p_{total} = p_{in} - p_{out}$, where the subscripts in denotes quantities inside the bubble (averaged over cells with $r = \sqrt{(x-x_0)^2 + (y-y_0)^2} \le R$) and out quantities outside the drop (averaged over cells with $r > R$),

2. $\Delta p_{partial} = p_{in} - p_{out}$, where the subscripts in denotes quantities inside the bubble (averaged over cells with $r \le \frac{R}{2}$ and out quantities outside the drop (averaged over cells with $r > \frac{3}{2}R$), so as to avoid the interface region,

3. $\Delta p_{max} = p_{max} - p_{min}$, where $p_{max}$ and $p_{min}$ are the maximum and minimum pressure on the whole domain, respectively.

We also monitor the degree of circularity, defined in [14] as

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

with $\Omega_2$ being the subdomain occupied by the bubble, $|\Omega_2|$ being the area of the bubble and $P_b$ being 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. Since, for a perfectly circular bubble, the degree of circularity is unitary, we expect values equal to one or very close to it.
### 4.1.1 Use of Laplace-Beltrami operator

In this Section, we consider relation (11) to evaluate $f_σ$ and we report the results in Table 1. Notice that we add a small number $η = 10^{-10}$ to the denominator $|∇φ|$ so as to avoid division by zero. The same workaround is adopted to evaluate the unit normal vector in (12). The approximation is very robust and no oscillations for the pressure jump arise during the time evolution (see Figure 3). However, one can easily notice that we achieve convergence towards a value which is different with respect to the analytical one. This is in agreement with what happens using the Volume of Fluid method, as reported in [8], and we will further discuss this issue in Section 4.1.3. Finally, the degree of circularity is always around 1 for all the configurations, meaning that the circular shape is preserved. The relative error for $χ$ with $N_{el}$ is around $2 × 10^{-5}$. Finally, the generation of spurious currents, which typically appear in the form of vortices around the interface, is strongly reduced, as evident from Figure 1. This further confirms the robustness of this approach.

<table><caption>Table 1: Static bubble test case, pressure jump across the interface and degree of circularity at final time $t = T_f$ with $f_σ$ computed using the Laplace-Beltrami operator (11). $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th>$N_{el}$</th><th>$\Delta p_{total}$</th><th>$\Delta p_{partial}$</th><th>$\Delta p_{max}$</th><th>$\chi$</th></tr></thead><tbody><tr><td>40</td><td>3.37</td><td>4.15</td><td>4.16</td><td>1.00</td></tr><tr><td>80</td><td>3.37</td><td>4.13</td><td>4.16</td><td>1.00</td></tr><tr><td>160</td><td>3.37</td><td>4.13</td><td>4.16</td><td>1.00</td></tr><tr><td>320</td><td>3.37</td><td>4.13</td><td>4.17</td><td>1.00</td></tr></tbody></table>

### 4.1.2 Computation of the total curvature from the level set

In this Section, we consider relations (12) for $f_σ$. Table 2 shows the results obtained with Bonnet's formula. One can easily notice that the results are much less accurate with respect to those obtained with the Laplace-Beltrami operator. Moreover, the values of pressure jump are strongly oscillating, as one can notice from Figure 3. The two relations in (12), namely the computation of $f_σ$ with Bonnet's formula and with the divergence of the unit normal are equivalent for continuous functions. However, they yield different values once projected onto discrete functional spaces. Notice also that, to apply the first relation in (12), we need to project the unit normal vector $\frac{\nabla\phi}{|\nabla\phi|}$ onto a suitable space, so as to compute the divergence. In spite of this discrepancy, analogous results are achieved as evident from Table 3 and Figure 3. The presence of spurious oscillations in the curvature field computed from the level set formulation is well known in the literature and different filtering approaches have been proposed to mitigate this issue, see e.g. [11] for a comparison different filtering strategies. For the degree of circularity, we notice a small degradation in the description of the interface, especially at coarse resolutions. However, good results are obtained overall. The relative error for $χ$ with $N_{el} = 320$ is around $4 × 10^{-4}$ for both relations (12).
<figure><img src="image_2.png" /><figcaption>Figure 1: Static bubble test case, velocity magnitude at $t = T_f$ using (11) with $N_{el} = 320$.</figcaption></figure>

<table><caption>Table 2: Static bubble test case, pressure jump across the interface and degree of circularity at final time $t = T_f$ with $f_σ$ computed using the Bonnet's formula (12). $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th>$N_{el}$</th><th>$\Delta p_{total}$</th><th>$\Delta p_{partial}$</th><th>$\Delta p_{max}$</th><th>$\chi$</th></tr></thead><tbody><tr><td>40</td><td>2.81</td><td>3.40</td><td>10.26</td><td>0.98</td></tr><tr><td>80</td><td>3.00</td><td>3.59</td><td>3.64</td><td>1.00</td></tr><tr><td>160</td><td>2.95</td><td>3.54</td><td>3.59</td><td>1.00</td></tr><tr><td>320</td><td>2.94</td><td>3.53</td><td>3.58</td><td>1.00</td></tr></tbody></table>

<table><caption>Table 3: Static bubble test case, pressure jump across the interface and degree of circularity at final time $t = T_f$ with $f_σ$ computed using the divergence of unit normal vector (12). $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th>$N_{el}$</th><th>$\Delta p_{total}$</th><th>$\Delta p_{partial}$</th><th>$\Delta p_{max}$</th><th>$\chi$</th></tr></thead><tbody><tr><td>40</td><td>2.26</td><td>2.81</td><td>9.62</td><td>0.98</td></tr><tr><td>80</td><td>2.98</td><td>3.57</td><td>3.62</td><td>1.00</td></tr><tr><td>160</td><td>2.95</td><td>3.54</td><td>3.59</td><td>1.00</td></tr><tr><td>320</td><td>2.94</td><td>3.53</td><td>3.58</td><td>1.00</td></tr></tbody></table>
### 4.1.3 Use of an evolution equation for the curvature

The use of the Laplace-Beltrami operator is a very robust technique, but, as already noticed, it yields slightly inaccurate results. The reason is that there are functions, such as the unit normal vector or the curvature, whose definitions are properly meaningful only for the points on the surface $\Gamma$ [19, 23]. The capillarity flux tensor ($\mathbf{I} - \mathbf{n}_\Gamma \otimes \mathbf{n}_\Gamma) |\nabla \phi|$ introduces instead spurious contributions far from the interface. Analogous considerations hold for the total curvature computed from the level set, for which spurious or even singular values arise if $|\nabla \phi| \approx 0$. In order to overcome this issue, we employ an evolution equation for the mean curvature, so as to evaluate $\mathbf{f}_\sigma$ with (16). This strategy is conceptually similar to the one proposed in [27], where a curvature-augmented approach to the level set method has been proposed. However, the evolution equation (14)-(15) is more general and valid for any sufficiently regular moving closed surface. A key point is the choice of the initial value for the mean curvature $H^0$, so as to avoid spurious contributions far from the interface. We set therefore

$$H^0 = \begin{cases} -\frac{1}{2} \frac{1}{\sqrt{(x-x_0)^2+(y-y_0)^2}} & \text{if } |\nabla \phi^0| > \frac{\beta}{h} \\ 0 & \text{otherwise.} \end{cases} \quad (35)$$

Here, $\beta = 5 \times 10^{-4}$, $\phi^0$ is the initial value of the level set function, and $h = \frac{L}{N_{el}r}$ is the space step size, with $L = 1$ m being the domain length, $N_{el}$ denoting the number of elements along each direction and $r$ being the polynomial degree of the finite element space chosen for the discretization of the mean curvature. Recall that, in the present work, we take $r = 2$. The resulting initial datum for the curvature with $N_{el} = 320$ is reported in Figure 2.

Relation (35) allows us to consider contributions of the curvature in $\mathbf{f}_\sigma$ only at the interface and close to it, namely when $|\nabla \phi|$ is above a certain threshold. Table 4 shows the obtained results and one can easily notice that we achieve more accurate results with respect to those obtained with the Laplace-Beltrami operator. Moreover, $\Delta p_{partial}$ and $\Delta p_{max}$ are closer as the resolution is increased, meaning that the presence of spurious contributions far from the interface is progressively reduced. No oscillations arise for the pressure jump, as evident from Figure 3, analogously to what is observed for the Laplace-Beltrami operator. Finally, for what concerns the degree of circularity, the relative error with $N_{el} = 320$ is around $1 \times 10^{-5}$, which is analogous to the value obtained with the Laplace-Beltrami operator and one order of magnitude lower than what obtained with (12).

<table><caption>Table 4: Static bubble test case, pressure jump across the interface and degree of circularity at final time $t = T_f$ with $\mathbf{f}_\sigma$ computed using (16). $N_{el}$ denotes the number of elements along each direction.</caption><thead><tr><th>$N_{el}$</th><th>$\Delta p_{total}$</th><th>$\Delta p_{partial}$</th><th>$\Delta p_{max}$</th><th>$\chi$</th></tr></thead><tbody><tr><td>40</td><td>3.53</td><td>4.36</td><td>4.40</td><td>0.99</td></tr><tr><td>80</td><td>3.38</td><td>4.13</td><td>4.15</td><td>1.00</td></tr><tr><td>160</td><td>3.36</td><td>4.10</td><td>4.11</td><td>1.00</td></tr><tr><td>320</td><td>3.35</td><td>4.05</td><td>4.06</td><td>1.00</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Static bubble test case, initialization of $H$ using (35) with $N_{el} = 320$.</figcaption></figure>

For longer times, the numerical solution of $H$ is corrupted. As discussed in [27], a reinitialization procedure is necessary, analogously to what is done for the level set and, more generally, for transport equations of interfacial quantities. This issue starts arising with the formation of spurious currents in the form of vortices around the interface, as one can notice from Figure 4. Nevertheless, the magnitude of spurious currents is significantly reduced with respect to (12), meaning that, in this test case, the present approach is more robust and accurate with respect to the one presented in Section 4.1.2. The development of suitable reinitialization techniques for (14)-(15) will be matter of future analysis.

## 4.2 Oscillating bubble

We now consider a more dynamic example. Starting from the configuration described in Section 4.1, we modify the initial shape from a circle to an ellipse by scaling the semi-axes by a factor 1.25 in the x-direction and 0.75 in the y-direction, respectively. The final time is now $T_f = 20$ s. Indeed, as already discussed in Section 4.1, for longer simulation times the approach based on the evolution equation (14)-(15) is corrupted. We consider two computational grids, composed by $N_{el} = 40$ and $N_{el} = 80$ elements, along each direction, respectively. We compare the use of the Laplace-Beltrami operator (11), the use of the Bonnet's formula (12) and the computation of surface tension contribution (16) using the evolution equation for the mean curvature (14)-(15). Figure 5 shows the interface $\phi = \frac{1}{2}$ at the final time. One can easily notice
<figure><img src="image_4.png" /><figcaption>Figure 3: Static bubble test case, evolution of $\Delta p_{partial}$ with $N_{el} = 320$. The black line reports the results obtained with the Laplace-Beltrami operator (11), the red line shows the results obtained with the Bonnet’s formula in (12), whereas the blue dots represent the results achieved with the divergence of the unit normal in (12). Finally, the green line shows the results obtained with the evolution equation (14)-(15) for the mean curvature and surface tension computed with (16).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Static bubble test case, velocity magnitude at $t = T_f$ with $N_{el} = 320$, a) curvature computed from divergence of unit normal as in (12), b) curvature computed from evolution equation (14)-(15) and surface tension evaluated with (16).</figcaption></figure>
that oscillations around the configuration captured by (11) are present using Bonnet's formula (12). These oscillations are stronger employing the evolution equation for the mean curvature (14)-(15). Analogous considerations hold using the 80 × 80 grid, even though the oscillations are reduced as long as the resolution increases, as one can notice from Figure 6.

<figure><img src="image_6.png" /><figcaption>Figure 5: Oscillating bubble test case, isoline $\phi = \frac{1}{2}$ at $t = T_f$ with $N_{el} = 40$. The black line shows the results with the use of the Laplace-Beltrami operator (11), the red line reports the results obtained with the use of the Bonnet's formula (12), whereas the blue line represents the results achieved using the evolution equation for the mean curvature (14)-(15) to compute (16). $N_{el}$ denotes the number of elements along each direction.</figcaption></figure>

# 5 Conclusions

We have performed a quantitative assessment of different strategies to compute the contribution due to surface tension in immiscible incompressible flows. The conservative level set method developed in [22] has been employed for this comparison. The results show that the use of an evolution equation for the curvature can be considered
<figure><img src="image_7.png" /><figcaption>Figure 6: Oscillating bubble test case, isoline $\phi = \frac{1}{2}$ at $t = T_f$ with $N_{el} = 80$. The black line shows the results with the use of the Laplace-Beltrami operator (11), the red line reports the results obtained with the use of the Bonnet’s formula (12), whereas the blue line represents the results achieved using the evolution equation for the mean curvature (14)-(15) to compute (16). $N_{el}$ denotes the number of elements along each direction.</figcaption></figure>

as a valid alternative to investigate in order to compute this quantity. The most accurate results presented in Section 4.1 are those obtained using the evolution equation for the mean curvature (14)-(15) to compute the surface tension contribution (16). Issues are instead present for longer times and for dynamic configurations, as discussed in Section 4.2, for which this approach, at the current stage, yields less accurate results and requires therefore further investigation. In future work, we aim to develop suitable reinitialization techniques for the evolution of the mean curvature, so as to improve the accuracy for longer simulation times and to incorporate this relation and the analogous one for the interfacial area density described in [23] into classical models for both incompressible and compressible two-phase flows, so as to improve the computation of interfacial source term in the case of a not well resolved interface.
# Acknowledgments

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] D. Arndt, W. Bangerth, M. Bergbauer, M. Feder, M. Fehling, J. Heinz, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, Pelteret J-P., B. Turksin, D. Wells, and S. Zampini. The deal.II library, version 9.5. *Journal of Numerical Mathematics*, 31(3):231–246, 2023.

[2] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II: a general-purpose object-oriented finite element library. *ACM Transactions on Mathematical Software (TOMS)*, 33:24–51, 2007.

[3] 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(21):11229, 2022.

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

[5] 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(2):267–279, 1997.

[6] J.U. Brackbill, D.B. Kothe, and C. Zemach. A continuum method for modeling surface tension. *Journal of Computational Physics*, 100(2):335–354, 1992.

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

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

[9] R. Estrada and R.P. Kanwal. Applications of distributional derivatives to wave propagation. *IMA Journal of Applied Mathematics*, 26(1):39–63, 1980.

[10] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, and M.W. Williams. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. *Journal of Computational Physics*, 213(1):141–173, 2006.

[11] D. Henneaux, W. Haegeman, P. Schrooyen, M. Massot, and T. Magin. Assessment of curvature filtering strategies for high-order sharp interface simulation of surface tension driven flows. Technical report, 2023.

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

[13] S. Hysing. A new implicit surface tension implementation for interfacial flows. *International Journal for Numerical Methods in Fluids*, 51(6):659–672, 2006.
[14] 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(11):1259-1288, 2009.

[15] J.B. Keller and M.J. Miksis. Surface tension driven flows. *SIAM Journal on Applied Mathematics*, 43(2):268-277, 1983.

[16] 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(1):134-147, 1994.

[17] V.G. Levich and V.S. Krylov. Surface-tension-driven phenomena. *Annual review of fluid mechanics*, 1(1):293-316, 1969.

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

[19] A. Nadim. A concise introduction to surface rheology with application to dilute emulsions of viscous drops. *Chemical Engineering Communications*, 148:391-407, 1996.

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

[21] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii. *Journal of Computational Physics*, 225(1):785-807, 2007.

[22] G. Orlando. An implicit DG solver for incompressible two-phase flows with an artificial compressibility formulation. *International Journal for Numerical Methods in Fluids*, 96:1932-1959, 2024.

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

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

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

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

[27] C.J. Vogl. A curvature-augmented, REA approach to the level set method. *SIAM Journal on Scientific Computing*, 38(2):A833-A855, 2016.