An unstructured finite element model for incompressible two-phase flow based on a monolithic conservative level set method

TL;DR

This paper introduces a robust, accurate, and efficient unstructured finite element model for incompressible two-phase flows using an improved monolithic level set method with new redistancing and stabilization techniques.

Contribution

It develops a novel PDE-based redistancing method and enhances the level set scheme for better accuracy and robustness in two-phase flow simulations.

Findings

The model conserves mass and volume on unstructured meshes.

It achieves high accuracy with minimal numerical tuning.

The method performs well across various benchmark problems.

Abstract

We present a robust numerical method for solving incompressible, immiscible two-phase flows. The method extends the monolithic phase conservative level set method with embedded redistancing by Quezada de Luna et al. [38] and a semi-implicit high-order projection scheme for variable-density flows by Guermond and Salgado [17]. The level set method can be initialized conveniently via a simple phase indicator field, which is pre-processed to obtain an approximate signed distance function. To do this, we propose a new PDE-based redistancing method. We also improve the scheme in [38] to provide more accuracy and robustness in full two-phase flow simulations. Specifically, we perform an extra step to ensure convergence to the signed distance level set function and simplify other aspects of the original scheme. Lastly, we introduce consistent artificial viscosity to stabilize the momentum equations in the context of the projection scheme. This stabilization is algebraic, has no tunable parameters and is suitable for unstructured meshes and arbitrary refinement levels. The overall methodology includes few numerical tuning parameters; however, for the wide range of problems that we solve, we identify only one parameter that strongly affects performance of the computational model and provide a value that provides accurate results across all the benchmarks presented. The result is a robust, accurate, and efficient two-phase flow model, which is mass- and volume-conserving on unstructured meshes and has low user input requirements for real applications.

Figures

Tables

Table 1. Table 1: Convergence of redistancing from discontinuous configuration.

h	N-DOFs	${\Vert {\varphi }_h^0-{\varphi }^{\text{exact}}\Vert }_{L^2\left(\Omega \right)}$	Rate
2.50E-2	1,681	6.75E-3	–
1.25E-2	6,561	2.68E-3	1.33
6.25E-3	25,921	1.15E-3	1.22
3.12E-3	103,041	5.88E-4	0.96
1.56E-3	410,881	3.05E-4	0.94
7.81E-4	1,640,961	1.28E-4	1.24

Table 2. Table 2: Convergence of two dimensional periodic vortex problem.

h	N-DOFs	${\Vert {\varphi }_h-{\varphi }^{\text{exact}}\Vert }_{L^1\left(\Omega \right)}$	Rate	${\Vert {\varphi }_h-{\varphi }^{\text{exact}}\Vert }_{L^2\left(\Omega \right)}$	Rate
2.50E-2	1,681	5.77E-2	–	7.52E-2	–
1.25E-2	6,561	1.87E-2	1.62	2.10E-2	1.84
6.25E-3	25,921	4.51E-3	2.05	4.78E-3	2.13
3.12E-3	103,041	1.33E-3	1.75	1.40E-3	1.77
1.56E-3	410,881	3.82E-4	1.80	3.53E-4	1.98

Table 3. Table 3: Error metrics ( 12 ) for the two dimensional periodic vortex problem. The initial length of the interface is L ≈ 9.42 × 10 − 1 𝐿 9.42 superscript 10 1 L\approx 9.42\times 10^{-1} .

$h$	N-DOFs	${\text{I}}_{\text{err}}$	${\text{V}}_{\text{err}}$	${\text{D}}_{\text{err}}$
2.50E-2	1,681	6.08E-2	2.26E-2	2.48E-3
1.25E-2	6,561	1.82E-2	7.68E-3	4.10E-4
6.25E-3	25,921	4.51E-3	3.76E-3	7.12E-5
3.12E-3	103,041	1.30E-3	1.39E-3	1.91E-5
1.56E-3	410,887	3.60E-4	1.74E-4	8.42E-6

Table 4. Table 4: Commonly used boundary conditions for the velocity.

Slip	Non-slip	Open top	Inflow
$𝐮\cdot 𝐧=0,{\left[\epsilon \left(𝐮\right)\cdot 𝐧\right]}_{\tau }=0$	$𝐮=0$	$\epsilon \left(𝐮\right)\cdot 𝐧=0,p=0$	$𝐮={𝐮}^{\text{inlet}}$, s.t.

Table 5. Table 5: Material parameters for the rising bubble problems. We consider two common test cases.

Test case	${\rho }_W$	${\rho }_A$	${\mu }_W$	${\mu }_A$	$g$	$\sigma$
1	1000	100	10	1	0.98	24.5
2	1000	1	10	0.1	0.98	1.96

Table 6. Table 6: Material parameters for the buckling flow problem.

${\rho }_W$	${\rho }_A$	${\mu }_W$	${\mu }_A$	$g$	$\sigma$
1800	1	500	$2\times {10}^{-5}$	9.8	0

Equations

$${\rho }^{n+1}$$

$${\mu }^{n+1}$$

$${\int }_{{\Omega }^{\ast }}(H(\varphi )-{\chi }_{\alpha })wd\mathbf{x}=0\forall w\in W({\Omega }^{\ast })$$

$${\int }_{{\Omega }^{\ast }}\left[(H(\varphi )-{\chi }_{\alpha })w+G\nabla w\right]d\mathbf{x}=0\forall w\in W({\Omega }^{\ast })$$

$${\partial }_t{S}_{ϵ}(\varphi )+\nabla \cdot \left[\mathbf{v}{S}_{ϵ}(\varphi )-\lambda (\nabla \varphi -\mathbf{q})\right]$$

$$\sqrt{|\nabla \varphi {|}^2+{\delta }^2}\mathbf{q}$$

$$S_{ϵ}(\varphi )$$

$$(\nabla \varphi -\mathbf{q})\cdot \mathbf{n}$$

$$S_{ϵ}(\varphi )$$

$$H_{ϵ}(\varphi )$$

$${\partial }_t{\int }_{\Omega }{S}_{ϵ}\left(\varphi (\mathbf{x},t)\right)d\mathbf{x}=0⟹{\partial }_t{\int }_{\Omega }{H}_{ϵ}\left(\varphi (\mathbf{x},t)\right)d\mathbf{x}=0,$$

$$\lambda =\tilde{\lambda }\left(\frac{v^{\max }}{C}\right)\frac{h(\mathbf{x})}{|||\varphi |-\overline{|\varphi |}|{|}_{L^{\infty }(\Omega )}},$$

$$\hat{\varphi }=\{\begin{array}{ll}1,& \text{ if }\text{leavevmode}\text{dist}(\mathbf{x},{\mathbf{x}}_0)<R,\\ 0,& \text{ if }\text{leavevmode}\text{dist}(\mathbf{x},{\mathbf{x}}_0)=0,\\ -1,& \text{ if }\text{leavevmode}\text{dist}(\mathbf{x},{\mathbf{x}}_0)>R,\end{array}$$

$${\int }_{\Omega }{S}_{ϵ}({\hat{\varphi }}_h)\left[|\nabla {\varphi }_h^{\ast }|-1\right]w(\mathbf{x})d\mathbf{x}+c{\int }_{\Omega }h(\mathbf{x})\nabla {\varphi }_h^{\ast }\cdot \nabla wd\mathbf{x}+{\int }_{\Omega }\tau [\nabla {\varphi }_h^{\ast }-{\mathbf{q}}_h({\varphi }_h^{\ast })]\cdot \nabla w(\mathbf{x})d\mathbf{x}=0,$$

$${\int }_{\Omega }{S}_{ϵ}({\varphi }_h^{\ast })\left[|\nabla {\varphi }_h^{\ast \ast }|-1\right]w(\mathbf{x})d\mathbf{x}+c{\int }_{\Omega }h(\mathbf{x})\nabla {\varphi }_h^{\ast \ast }\cdot \nabla wd\mathbf{x}+{\int }_{\Omega }\tau [\nabla {\varphi }_h^{\ast \ast }-{\mathbf{q}}_h({\varphi }_h^{\ast \ast })]\cdot \nabla w(\mathbf{x})d\mathbf{x}=0,$$

$$\begin{array}{rl}\alpha m_i^{{\Gamma }^{\ast \ast }}({\Phi }_i^0-{\Phi }_i^{\ast \ast })+{\int }_{\Omega }{S}_{ϵ}({\varphi }_h^{\ast \ast })\left[|\nabla {\varphi }_h^0|-1\right]w(\mathbf{x})d\mathbf{x}+c{\int }_{\Omega }h(\mathbf{x})\nabla {\varphi }_h^0\cdot \nabla wd\mathbf{x}& \\ +{\int }_{\Omega }\tau [\nabla {\varphi }_h^0-{\mathbf{q}}_h({\varphi }_h^0)]\cdot \nabla w(\mathbf{x})d\mathbf{x}=0,& \forall w(\mathbf{x})\in X_h^1.\end{array}$$

$$J_{ij}^k:=\alpha m_i^{{\Gamma }^{\ast \ast }}{\delta }_{ij}+{\int }_{\Omega }{S}_{ϵ}({\varphi }_h^{\ast \ast })\frac{(\nabla {\varphi }_h^k\cdot \nabla w_j)}{|\nabla {\varphi }^k|+10^{-10}}{w}_{i}d\mathbf{x}+{\int }_{\Omega }\left[ch(\mathbf{x})+\tau \right]\nabla w_i\cdot \nabla w_{j}d\mathbf{x},$$

$${\mathbf{q}}_h({\varphi }_h)=\sum _j{\mathbf{Q}}_j({\varphi }_h)w_j(\mathbf{x}),$$

$$Q_i^{(k)}({\varphi }_h)=\frac{{\int }_{\Omega }{\partial }_k{\varphi }_h{w}_{i}d\mathbf{x}}{{\int }_{\Omega }\sqrt{|\nabla {\varphi }_h{|}^2+{\delta }^2}{w}_{i}d\mathbf{x}}.$$

$$\alpha m_i^{{\Gamma }^n}({\Phi }_i^{\ast }-{\Phi }_i^n)+{\int }_{\Omega }{S}_{ϵ}({\varphi }_h^n)\left[|\nabla {\varphi }_h^{\ast }|-1\right]w(\mathbf{x})d\mathbf{x}+{\int }_{\Omega }\tau [\nabla {\varphi }_h^{\ast }-{\mathbf{q}}_h({\varphi }_h^n)]\cdot \nabla w(\mathbf{x})d\mathbf{x}$$

$$\begin{array}{r}R({\varphi }_h^{n+1},w):={\int }_{\Omega }\frac{S_{ϵ}({\varphi }_h^{n+1})-S_{ϵ}({\varphi }_h^n)}{\Delta t}wd\mathbf{x}-\frac{1}{2}{\int }_{\Omega }\left[S_{ϵ}({\varphi }_h^n){\mathbf{v}}^n+S_{ϵ}({\varphi }_h^{n+1}){\mathbf{v}}^{n+1}\right]\cdot \nabla wd\mathbf{x}\\ +{\int }_{\Omega }\lambda [\nabla {\varphi }_h^{n+1}-{\mathbf{q}}_h({\varphi }_h^{\ast })]\cdot \nabla wd\mathbf{x}+\frac{1}{2}{\int }_{\partial {\Omega }^{-}}\left[S_{ϵ}({\varphi }_h^n){\mathbf{v}}^n+S_{ϵ}({\varphi }_h^{n+1}){\mathbf{v}}^{n+1}\right]\cdot \mathbf{n}wd\mathbf{s}=0,\forall w\in X_h^1,\end{array}$$

$$\frac{\partial R({\varphi }_h^k,w_i)}{\partial {\Phi }_j}=:J_{ij}^k={\int }_{\Omega }{S}_{ϵ}^{\prime }({\varphi }_h^k)\left[\frac{1}{\Delta t}{w}_i{w}_j-\frac{1}{2}\nabla w_i\cdot ({\mathbf{v}}^{n+1}{w}_j)\right]d\mathbf{x}+{\int }_{\Omega }\lambda \nabla w_i\cdot \nabla w_{j}d\mathbf{x}$$

$$+\frac{1}{2}{\int }_{\partial \Omega }{S}_{ϵ}^{\prime }({\varphi }_h^k)w_i{w}_j{\mathbf{v}}^{n+1}\cdot \mathbf{n}d\mathbf{s}.$$

$$\varphi (\mathbf{x},0)$$

$$\mathbf{v}(x,y,t)$$

$${\text{I}}_{\text{err}}({\varphi }_h)$$

$${\text{V}}_{\text{err}}({\varphi }_h)$$

$${\text{D}}_{\text{err}}({\varphi }_h)$$

$$\rho ({\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u})-\nabla \cdot \left(2\mu \epsilon (\mathbf{u})\right)+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\text{BDF2}({\mathbf{u}}^{n+1})$$

$${\mathbf{u}}^{\ast }$$

$$p^{\ast }$$

$${\rho }^{n+1}\left[\text{BDF2}({\mathbf{u}}^{n+1})+({\mathbf{u}}^{\ast }\cdot \nabla ){\mathbf{u}}^{n+1}\right]-\nabla \cdot (2{\mu }^{n+1}\epsilon ({\mathbf{u}}^{n+1}))+\nabla p^{\ast }$$

$$\mathbf{u}$$

$$\epsilon ({\mathbf{u}}^{n+1})\mathbf{n}$$

$$\mathbf{u}\cdot \mathbf{n}=0,[\epsilon ({\mathbf{u}}^{n+1})\mathbf{n}{]}_{\tau }$$

$$-\Delta \delta {\psi }^{n+1}$$

$$\delta {\psi }^{n+1}$$

$$\nabla \delta {\psi }^{n+1}\cdot \mathbf{n}{|}_{\partial \Omega }$$

Full text

1 Introduction

Understanding the behavior of immiscible fluids with distinct material characteristics (e.g. water and air) is important for many engineering and industrial applications. Water-oil-gas interaction within subterranean reservoirs, combustion engines, and open channel hydraulics are all examples of flow phenomena involving multiple fluids, i.e. multiphase flow. In order to accurately describe the evolving interface between individual fluid subdomains, high fidelity numerical models are required. The methods involved should aim for computational efficiency along with preservation of qualitative properties of the material continuum. In the case of multiphase simulations concerning divergence-free velocity fields, maintaining the conservation of mass and volume is imperative.

In general, multiphase flow is modeled by solving the Navier-Stokes equations with spatially varying material parameters. We consider an Eulerian description for the fluid motion along with continuous Galerkin finite elements for discretization of all equations in space. For Eulerian methods, the fluid motion is calculated from a stationary frame of reference, a fixed computational grid. These techniques can resolve large fluid deformations, however, a representation of the phase interface is required.

Two routinely implemented methods for representing and evolving interfaces between discrete phases are the Volume of Fluid (VOF) method by [20] and level set techniques by [37, 42]. The VOF method assigns phase identities, e.g. fluids A and B, via a characteristic function that equals one in fluid A and zero in fluid B. Grid cell averaging of this phase indicator function results in the designation of a cell volume fraction in the range [0, 1]. An interface is then reconstructed in cells with intermediate volume fraction and the VOF re-initialized. For many multiphase flow applications, the characteristic function is transported using velocities obtained from the (incompressible) Navier-Stokes equations. This solution produces density and viscosity fields requiring subsequent interfacial reconstruction. Once the phase boundary location is accurately reconstructed, material properties can then be designated respectively for each fluid subdomain.

In the level set method, the interface between the phases is represented implicitly by a level set of a scalar function defined on the entire domain; e.g., as the zero level set of the signed distance function (SDF) to the interface. The SDF returns a zero value at the interface and either a positive or a negative distance for subdomains fluid A and fluid B respectively. The level set function is advected by the fluid velocity field produced by the (incompressible) Navier-Stokes equations.

Both the VOF and the level set methods bear limitations when implemented alone. Despite maintaining volume (or phase) conservation, the VOF requires reconstruction of interface locations from less precise cell averages. The interface is effectively smeared over some distance specific to the numerical scheme and then reconstructed to be sharp base on geometric approximations. Level set methods do not require interface reconstruction; however, they lack a discrete conservation property; phases enclosed by the interface can suffer from a significant loss of volume due to the accumulation of discrete conservation errors over time.

There are hybrid methods combining ideas of the VOF, level set, and particle methods. See for instance [41, 23, 13]. In this work we consider a new hybrid level set and volume-of-fluid that builds on the weak variational form of phase conservation proposed in [24], the optimal control formulation of [6], and the recent reformulation and extension of these approaches to a monolithic scheme in [38]. These methods simultaneously maintain the convenient and precise signed distance represenation of the phase geometry and a phase conservation property.

Despite the numerical convenience of the signed distance field, it can be difficult for model users to provide initial conditions, since it generally requires solution of the nonlinear Eikonal equation from some other description of the interface. In this work we provide a method that only requires input of the initial phase indicator field, which we pre-process to obtain the SDF. To do this, we solve the non-linear Eikonal equation with a stabilization term that blends ideas from [31], the elliptic redistancing by [6] and [38]. In addition, we modify the method in [38] by simplifying the numerical discretization and solving (at every time step) a redistancing pre-stage based on an extensions of the stabilized Hamilton-Jacobi formulation of redistancing from [24], which ensures convergence of the elliptic redistancing approach to the correct SDF.

The velocity field is obtained via the second order projection scheme for the incompressible Navier-Stokes equations with variable material parameters by [17]. We use Taylor-Hood finite elements. Additionally, we incorporate artificial viscosity based on [28, 1, 5, 16] and surface tension as proposed in [21]. The artificial viscosity has no tunable parameters and is algebraic, which makes it robust and suitable for unstructured meshes and arbitrary refinement levels.

The rest of this work is organized as follows. In §2 we describe the finite element spatial discretization. Afterwards, in §3 we describe the overall algorithm. Later, in §4 we review the conservative level set method by [38]. We consider only the continuous model and propose a simpler discretization in time. Furthermore, we propose a redistancing pre-stage that can be used to obtain the SDF from a discontinuous indicator field. In §5 we describe the full Navier-Stokes discretization. In this section we propose an algebraic and robust stabilization for the momentum equations. Section 6 is devoted to the numerical examples. Finally, we close with some discussion in §7.

2 Finite element spatial discretization

We use continuous Galerkin finite elements to discretize all equations in space. Let $d=\{1,2,3\}$ be the number of spatial dimensions and $\Omega\subset\mathbb{R}^{d}$ be a bounded domain with boundary $\partial\Omega\subset\mathbb{R}^{d-1}$ which we decompose into $\partial\Omega^{-}=\{\mathbf{x}\in\partial\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \mathbf{u}\cdot\mathbf{n}<0\}$ and $\partial\Omega^{+}=\{\mathbf{x}\in\partial\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \mathbf{u}\cdot\mathbf{n}\geq 0\}$. Given an SDF $\phi$ we define $\Gamma(t)=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \phi(\mathbf{x},t)=0\}$ to be the material interface. Time-dependent variables are defined on the time interval $t\in[0,T]$, where $T>0$. Given a computational mesh $\mathcal{T}_{h}$, we consider the finite element space $X_{h}^{p}=\{w\in\mathcal{C}^{0}(\Omega)\leavevmode\nobreak\ |\leavevmode\nobreak\ w|_{K}\in\mathbb{P}_{\text{p}},\forall K\in\mathcal{T}_{h}\}$ with $p=\{1,2\}$; i.e., we use only continuous piecewise linear and quadratic spaces. The spaces are spanned by basis functions $\{w_{1},\dots,w_{\text{dim}\left(X_{h}^{p}\right)}\}$ which possess the partition of unity property; i.e., $\sum_{j}w_{j}(\mathbf{x})=1$. The degrees of freedom associated with these basis functions are denoted by uppercase letters. The finite element solution $u_{h}(\mathbf{x})\in X_{h}$ is given by $u_{h}(\mathbf{x})=\sum_{j\in\mathcal{I}(\Omega_{i})}U_{j}w_{j}(\mathbf{x})$, where $\Omega_{i}$ is the patch of elements containing node $i$. Here, and in the rest of this paper, the notation $\mathcal{I}(z)$ is used for the index set containing the numbers of all basis functions whose support on $z$ is of nonzero measure.

3 Two-phase flow algorithm based on operator splitting

The overall two-phase flow algorithm is driven by a simple first-order operator splitting technique, as in [30, 40, 24]. In particular, operator splitting is used to decouple the level set and the Navier-Stokes stages. More sophisticated alternatives are possible as well, see for instance [27, 35, 33]. Following this approach one can discretize first the level set equation and subsequently the Navier-Stokes equations or vice versa; we choose the former, which is analogous to the splitting presented for variable-density flows in [17]. In the rest of this section, for simplicity of exposition, we assume a fixed time step; however, we apply this algorithm to a variable time stepping scheme, see §5.1. Let $\phi$, $\mathbf{u}$ and $p$ denote the level set, velocity and pressure fields respectively. For any given time step assume we know the solution at time $t^{n}$ and $t^{n-1}$; in particular, that we know $\phi^{n}$, $\mathbf{u}^{n}$, $\mathbf{u}^{n-1}$, $p^{n}$ and $p^{n-1}$. We proceed as follows:

Obtain a second order extrapolation of the velocity field: $\mathbf{u}^{*}=2\mathbf{u}^{n}-\mathbf{u}^{n-1}\approx\mathbf{u}^{n+1}$. 2. 2.

Using $\mathbf{u}^{*}$, solve the level set equation described in §4, to obtain $\phi^{n+1}$. 3. 3.

Given the level set solution at $t^{n+1}$, define the fluid density ($\rho^{n+1}$) and dynamic viscosity ($\mu^{n+1}$) via

[TABLE]

where $A$ and $W$ denote the air and water phases respectively, and $H_{\epsilon}$ is a regularized Heaviside function defined in §4. 4. 4.

Given $\rho^{n+1}$ and $\mu^{n+1}$, solve the variable-coefficient Navier-Stokes equations via a projection scheme, see §5.1, to obtain an updated velocity field $\mathbf{u}^{n+1}$. 5. 5.

Repeat until the final time is reached.

4 Monolithic conservative level set method

First note that a weak formulation of phase volume conservation for phase $\alpha$, with phase indicator function $\chi_{\alpha}$, and boundary described by the zero level set of an SDF $\phi$, can be written as

[TABLE]

where $W(\Omega^{*})$ is any suitable test space containing the constant function on subsets $\Omega^{*}$ of $\Omega$. For example, global conservation is ensured if $\Omega^{*}=\Omega$ or local conservation if $\Omega^{*}=K$. Equation (1) is ill-posed as an equation for $\phi$, but, if $W(\Omega^{*})$ has the partition of unity property, then the conservation property is maintained even after addition of some regularization $G$ [24]:

[TABLE]

In this work we consider a time-dependent form of equation (2) from [38] for solenoidal velocity field $\mathbf{v}$, which is given in strong form by

[TABLE]

where $\phi\in\mathbb{R}$ denotes the SDF level set function, $S_{\epsilon}$ is a smoothed sign function, $\lambda>0$ is a user defined parameter and $\delta>0$ is a small regularization parameter. We use $\delta=1\times 10^{-10}$ in all simulations. The smoothed sign function makes conservation symmetric with respect to the phases and the regularization is based on an elliptic form of the Eikonal equation[6]. Equation (3) combines level set evolution, signed distance property, and VOF evolution into a single equation as described further below. We follow [38] and consider

[TABLE]

Let $h(\mathbf{x})$ denote a characterization of the local mesh size, then $\epsilon=\alpha h(\mathbf{x})$ defines the thickness of the regularization of $H_{\epsilon}$ and $S_{\epsilon}$. We use $\alpha=\frac{3}{2}$ in all simulations.

Equation (3a) is a conservation law for $S_{\epsilon}(\phi)$. Under appropriate boundary conditions and a conservative finite element method (e.g. the partition of unity property), this implies that

[TABLE]

and since $\int_{\Omega}H_{\epsilon}(\phi)d\mathbf{x}$ represents the (regularized) volume of one of the phases, the method is volume conservative. The terms $\partial_{t}S_{\epsilon}(\phi)+\nabla\cdot\mathbf{v}S_{\epsilon}(\phi)=0$ correspond to a volume of fluid like model that impose not only conservation but are also responsible for advecting the interface to the correct position. This is true since the velocity field is assumed to be solenoidal ($\nabla\cdot\mathbf{v}=0$) and $\partial_{t}S_{\epsilon}(\phi)+\mathbf{v}\cdot\nabla S_{\epsilon}(\phi)=0\implies\partial_{t}\phi+\mathbf{v}\cdot\nabla\phi=0$. Therefore, these terms are important for consistency with the volume of fluid and level set methods. Henceforth, we refer to them as consistency terms of (3a). Since the Jacobian of $\partial_{t}S_{\epsilon}(\phi)+\nabla\cdot\mathbf{v}S_{\epsilon}(\phi)=0$ is singular, the terms $-\nabla\cdot\lambda(\nabla\phi-\mathbf{q})$ are added to the model; i.e., they act as a regularization. Note that $\mathbf{q}\approx\frac{\nabla\phi}{|\nabla\phi|}$; hence, the regularization terms also penalize deviations of the level set from the distance function. The result is a model for the SDF level set that is volume conservative and contains a term that penalizes deviations from the distance function.

Remark 4.0.1** (Boundary conditions).**

Note that technically one has to impose boundary conditions for the level set $\phi$ in the inflow boundary $\partial\Omega^{-}$. However, this information is only applied through the smoothed sign function $S_{\epsilon}$. Therefore, all we need to know is $S_{\epsilon}\left(\phi^{\text{BC}}\right)$; i.e., the phase (e.g., water or air) at that particular inflow boundary. The second boundary condition $(\nabla\phi-\mathbf{q})\cdot\mathbf{n}=0,\leavevmode\nobreak\ \forall\mathbf{x}\in\partial\Omega$, is imposed to guarantee conservation. Since the regularization term is integrated by parts during the finite element spatial discretization, applying this boundary condition is trivial.

Remark 4.0.2** (About the parameter $\lambda$).**

The parameter $\lambda$ controls the amount of regularization and penalization introduced to the model. For any given problem, it is important to guarantee consistency of (3a) w.r.t. the volume of fluid equation. Therefore, the authors in [38] propose to scale $\lambda$ by $h(\mathbf{x})$. We follow their definition and use

[TABLE]

where $\tilde{\lambda}=\mathcal{O}(1)$ is a dimensionless user defined parameter, $v^{\max}=\max_{\mathbf{x}}|\mathbf{v}|$, $C<1$ is the Courant number and $\overline{|\phi|}=\frac{1}{|\Omega|}\int_{\Omega}|\phi|d\mathbf{x}$.

We identify this parameter as the most important in the overall methodology. In all the problems in §6 we obtain satisfactory results using $\tilde{\lambda}=10$. However, it is essential to realize that the size of $|||\phi|-\overline{|\phi|}||_{L^{\infty}(\Omega)}$ changes depending on the domain. In particular, if the domain allows the maximum value of $|\phi|$ to grow from one problem to another, then the effective value of $\lambda$ is reduced. If $\lambda$ is too small we observe perturbations in the free surface. And if $\lambda$ is too large more dissipation is added. In general, for any given problem, we want to choose $\tilde{\lambda}$ to be as small as possible without introducing large perturbations to the free surface. In [38], the authors propose to use automated control theory to adjust this parameter at every time step. See for instance [6].

4.1 Discontinuous initialization of the distance function level set

The conservative level set model (3) requires an initial condition $\phi(\mathbf{x},t=0)$ to be given by an SDF. From a practical point of view this might be inconvenient in some situations. Instead, we propose to start the algorithm by considering an initial configuration given by a discontinuous function $\hat{\phi}\in\{-1,0,1\}$ representing the two phases and re-distance it to obtain an SDF. Due to the fact that the underlying equations for the SDF are nonlinear and the proposed initial conditions are far from the root, we use a robust multi-stage approach to ensure convergence and accuracy. This somewhat complex approach is used only at initialization and requires no user input beyond the initial phase configuration. The main objectives is a practical and robust SDF calculation on unstructured meshes. See figure 1a where we consider $\Omega=(0,1)^{2}$ and

[TABLE]

where $\mathbf{x}_{0}=(0.5,0.5)$ and $R=0.25$. First we project $\hat{\phi}$ onto the finite element space. Afterwards, the main algorithm is based on solving the Eikonal equation several times. The first time we re-distance the solution away from the interface. Then we concentrate on the cells containing the interface. And finally we perform a global redistancing to polish the result. We now provide details about this process.

Step 1: Lumped $L^{2}$-projection. We start by doing a lumped $L^{2}$-projection to obtain $\hat{\phi}_{h}=\sum_{j}\hat{\Phi}_{j}w_{j}(\mathbf{x})$ where, upon defining $m_{i}=\int_{\Omega}w_{i}d\mathbf{x}$, $\hat{\Phi}_{i}=\frac{h_{e}}{m_{i}}\int_{\Omega}\hat{\phi}w_{i}d\mathbf{x}$. Here $h_{e}$ is a characteristic element size; e.g., $h_{e}=\frac{1}{2}\left[\max_{\mathbf{x}\in\Omega}h(\mathbf{x})+\min_{\mathbf{x}\in\Omega}h(\mathbf{x})\right]$. Note that due to the scaling by $h_{e}$, $\hat{\phi}_{h}$ has units of distance. This step is important to introduce some dissipation into the initial condition. In figure 1c we show $\hat{\phi}_{h}$.

Step 2: Redistancing away from the interface. Given $\hat{\phi}_{h}\in X_{h}^{1}$, we solve the (viscous) Eikonal equation away from the interface to obtain $\phi_{h}^{*}\in X_{h}^{1}$ as follows:

[TABLE]

During this process, we freeze the DOFs associated with the interface $\hat{\Gamma}=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \hat{\phi}_{h}=0\}$ by imposing strongly $\Phi_{i}^{*}=\hat{\Phi}_{i},\leavevmode\nobreak\ \forall i\in\mathcal{I}(\hat{\Gamma})$. The first term in (5) is the consistency term w.r.t. the Eikonal equation. In the second term $c=\mathcal{O}(1)$, we use $c=0.1$, controls the amount of background dissipation. We introduce this dissipation to aim the solution to converge to the viscous solution of the Eikonal equation. The third term acts as non-linear stabilization and penalization from the distance function, see remark 4.1.1. In §4.2.1 we provide details about the discretization of $\mathbf{q}_{h}(\cdot)$. In figure 1d we show the result of this step.

Step 3: Redistancing close to the interface. Now we re-distance the solution close to the interface. Given $\phi^{*}_{h}\in X_{h}^{1}$, we solve the Eikonal equation to obtain $\phi^{**}_{h}\in X_{h}^{1}$ as follows:

[TABLE]

During this process, we impose strongly $\Phi_{i}^{**}=\Phi_{i}^{*},\leavevmode\nobreak\ \forall i\in\mathcal{I}(\Omega)\setminus\mathcal{I}(\hat{\Gamma})$; i.e., we allow changes only on the DOFs associated to the interface. By doing this, we aim to limit drastic movement of the interface. We perform only one Newton iteration of this step. In figure 1e we plot $\phi_{h}^{**}$.

Step 4: Global redistancing. Finally we obtain $\phi_{h}^{0}\in X_{h}^{1}$ by solving

[TABLE]

The first term, inspired by [6], is a penalization term to prevent large changes in the interface. Here $\alpha\in\mathbb{R}$ is a penalization constant, we use $\alpha=10^{9}$, and $m_{i}^{\Gamma^{**}}=\int_{\Omega}\delta_{\epsilon}(\phi^{**}_{h})w_{i}d\mathbf{x}\approx\int_{\Gamma^{**}}w_{i}d\mathbf{s}$, where $\Gamma^{**}=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \phi_{h}^{**}=0\}$. We perform only one Newton iteration of this step. We show $\phi_{h}^{0}$ in figure 1f.

To solve equations (5)-(7) we use a quasi-Newton method. Here we show only the approximate Jacobian of equation (7). The approximate Jacobian at the $k$-th Newton iteration is given by

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function. Note that, for simplicity, we drop the part related to $\mathbf{q}_{h}(\phi_{h}^{0})$.

Remark 4.1.1** (About the stabilization term).**

The last term in equations (5)-(7) acts as stabilization. Indeed $\nabla\phi-\mathbf{q}(\phi)\approx\left(1-\frac{1}{|\nabla\phi|}\right)\nabla\phi$. Therefore, this term behaves like the weak discretization of a Laplacian term with coefficient given by the residual of the Eikonal equation. This idea is similar to the high-order stabilization proposed in [31, §5] and the nonlinear residual-based variational multiscale stabilization used in [24]. The parameter $\tau$ controls the strength of the stabilization. We use $\tau=\frac{h(\mathbf{x})}{2}|{\bf V}|$ with $|{\bf V}|=1$ which is reasonable since ${\bf V}=S_{\epsilon}(\phi)\frac{\nabla\phi}{|\nabla\phi|}$ is the redistancing velocity. Note that this stabilization term is similar to the regularization and penalization terms in (3a). Moreover, the terms in (3a) can be interpreted as nonlinear stabilization where the coefficient is given not by the residual of the level set equation but by the residual of the Eikonal equation. See [38, Remark 3.2.2].

Remark 4.1.2** (About the sign of the stabilization term).**

Note that \left(1-\frac{1}{|\nabla\phi|}\right)<0\leavevmode\nobreak\ if $\leavevmode\nobreak\ |\nabla\phi|<1$. In this situation anti-diffusion is applied which helps in the convergence to $|\nabla\phi|=1$. This is similar to parabolic redistancing by [10] and elliptic redistancing by [6]. Moreover, equations (5)-(7) can be seen as a blend between hyperbolic and elliptic redistancing. By choosing $\tau=\mathcal{O}(h)$ we favor hyperbolic redistancing.

We consider again $\hat{\phi}$ to be given by (4) and perform a convergence test. The results are shown in table 1.

4.2 Full discretization of the conservative level set method

4.2.1 $\mathcal{C}^{0}$ normal reconstruction

The vector field $\mathbf{q}_{h}(\phi_{h})$ in (3b) approximates a normal field to the interface. We follow [38] and consider a weighted lumped $L^{2}$-projection given by

[TABLE]

4.2.2 Redistancing pre-stage

We begin each time step with a pre-redistancing stage using the algorithm in §4.1. The motivation behind this step is to introduce some hyperbolicity into the redistancing process. Our aim is prompt the redistancing to emanate from the interface. In our experience, this is not always possible with the elliptic redistancing, and thus not always possible with the penalization embedded in (3a). Given the solution $\phi^{n}_{h}$ at time $t^{n}$ we find $\phi^{*}_{h}\in X_{h}^{1}$ such that

[TABLE]

where $m_{i}^{\Gamma^{n}}=\int_{\Omega}\delta_{\epsilon}(\phi^{n}_{h})w_{i}d\mathbf{x}$, $\alpha=10^{9}$ and $\tau$ is given as in remark 4.1.1. To freeze the interface $\Gamma^{n}=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \phi_{h}^{n}=0\}$, we can set $\alpha=0$ and impose strongly $\Phi_{i}^{*}=\Phi_{i}^{n},\leavevmode\nobreak\ \forall i\in\mathcal{I}(\Gamma^{n})$. In this case, contrary to §4.1, we keep $\mathbf{q}_{h}$ independent of the solution $\phi^{*}_{h}$. This is done to avoid extra computational effort during this pre-stage.

4.2.3 Discretization of the level set equation

In [38] it is noted that one can use linear continuous Galerkin finite elements with no extra stabilization provided that the advection term is treated implicitly. We follow this idea and solve model (3a) via

[TABLE]

where $\phi_{h}^{*}$ is the result of the pre-stage given by (9). Note that we discretize the consistency terms in (3a) via a second order Crank-Nicolson discretization in time and use a first-order, implicit-explicit discretization for the regularization and penalization terms respectively. This approach keeps the method simple and efficient (compared to the two stage method proposed in [38]). Indeed, by doing a linearization of (3a) around the interface, the regularization and penalization terms are expected to be $\mathcal{O}(h^{2})$, see [38, Remark 3.2.2]. Therefore, one can expect no harm from being permissive with respect to the order of approximation of such terms. We solve (10) via Newton’s method. The Jacobian corresponding to $R(\phi^{n+1}_{h},w_{i})$ at the $k$-th Newton iteration is given by

[TABLE]

We close this section by solving a benchmark in the literature of level sets. The problem is known as periodic vortex, see [39]. The domain is given by $\Omega=(0,1)^{2}$. The initial condition and velocity field are given by

[TABLE]

where $\Gamma_{0}:=\{(x,y)\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ (x-x_{c})^{2}+(y-y_{c})^{2}=r^{2}\}$ is a circle of radius $r=0.15$ centered at $(x_{c},y_{c})=(0.5,0.75)$. We select the positive distance in (11a) if $\mathbf{x}=(x,y)$ is inside the circle $\Gamma_{0}$ and the negative distance otherwise. In figure 2 we show the solution at different times, and in table 2 we show the errors and convergence rates in the $L^{1}$- and $L^{2}$-norms for different refinements. Additionally, we compute some of the metrics in [38, §5] given by

[TABLE]

where $L$ is the $(d-1)$-dimensional measure of the zero level set $\Gamma(0)$, and $H(\cdot)$ is the sharp Heaviside function. The quantities $\text{I}_{\text{err}}(\phi_{h})$ and $\text{V}_{\text{err}}(\phi_{h})$ measure the extent of interface displacements and area/volume difference with respect to the exact SDF and exact Heaviside function, respectively, see [14, 36], while $\text{D}_{\text{err}}(\phi_{h})$ measures the deviation of $\phi_{h}$ from a distance function. We present the results in table 3 for five levels of refinement.

5 Incompressible Navier-Stokes solver

The incompressible, non-conservative form of the Navier-Stokes equations with variable material parameters is given by

[TABLE]

where $\varepsilon(\mathbf{u})=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)$ is the symmetric gradient, $\rho$ and $\mu$ are the density and viscosity respectively and $\mathbf{u}$, $p$ and $\mathbf{f}$ are the velocity, pressure and force fields respectively.

5.1 Time discretization via a projection scheme

We consider the second order projection scheme for the incompressible Navier-Stokes equations with variable density by [17] and adapt it for variable time step sizes. Let $r:=\frac{\Delta t^{n}}{\Delta t^{n-1}}$. Upon defining

[TABLE]

See §6 for details about common initial and boundary conditions.

Remark 5.1.1** (About the linearization of the momentum equations).**

Note that we follow [17] and use $\mathbf{u}^{*}$, a second order extrapolation of the velocity field, to linearize the momentum equations. Moreover, we use $\mathbf{u}^{*}$ in the next section to incorporate artificial viscosity to stabilize the advective term. By doing this, the equation remains linear with respect to the solution at time $t^{n+1}$.

Remark 5.1.2** (About accuracy).**

The expected accuracy of this projection scheme is second order in the $L^{2}$- and $H^{1}$-norms. We refer the reader to [17] for a convergence study.

5.2 Artificial viscosity

In this section, we add artificial viscosity to stabilize the advective term in the momentum equations. To do this we consider [28, 1, 5, 16], where artificial dissipative operators are introduced to enforce Discrete Maximum Principles (DMP). We do not have a DMP requirement; instead, we are interested in adding enough artificial viscosity to have a robust and well behaved discretization. This extra viscosity must vanish as the mesh is refined, must behave properly for unstructured meshes and for all refinements and should not add significant computational cost to the method. In particular, we are interested in keeping equation (14d) linear. Finally, we aim to have no tunable parameters and to preserve, to the best of our ability, the accuracy properties of the underlying method.

Componentwise smoothness indicator

We concentrate first on one component of the velocity vector; e.g., in the $x$-component $u_{h}$, whose degrees of freedom are denoted by $U$, see §2. The first step is to neglect the force, pressure and viscosity terms in (14d); as a result, we get a hyperbolic system, which is more suitable for the theory developed in the references above. In the rest of this section we consider backward Euler time stepping; nevertheless, we apply these results to the second order method in §5.1. By doing this the full discretization becomes

[TABLE]

where $M$ is the mass matrix and $K$ is the linearized advective matrix, whose components are given by $M_{ij}=\int_{\Omega}w_{i}(\mathbf{x})w_{j}(\mathbf{x})d\mathbf{x}$ and $K_{ij}=\int_{\Omega}\left[\mathbf{u}_{h}^{*}\cdot\nabla w_{j}(\mathbf{x})\right]w_{i}(\mathbf{x})d\mathbf{x}$ respectively. Now we introduce an artificial dissipative matrix, see e.g. [28] and references therein, whose components are given by

[TABLE]

The objective behind this idea is to add sufficient conditions to construct $U^{n+1}$ as a convex combination of $U^{n}$. This, along with using a positive lumped mass matrix, guarantees the DMP property. By doing this, however, the accuracy of the method is reduced to first order. To recover the second order accuracy expected from the projection scheme we use a smoothness based indicator, see [1, 5, 16], given by

[TABLE]

The indicator $0\leq\beta_{i}^{x}\leq 1$, which is associated with the $x$-component of the velocity field, is (close to) one when the solution is smooth and (close to) zero when the solution is non-smooth and at local extrema.

Isotropic artificial dissipation

The smoothness indicators for the other velocity components are computed similarly. Assume the problem is three dimensional. Then, we define

[TABLE]

which is subsequently used to reduce the amount of artificial dissipation from (16). This can be done in different ways. In this work, we follow [1] and define

[TABLE]

Finally, we apply this discrete operator to each component of the velocity and reincorporate the pressure, force and viscosity terms. The full discretization for the $x$-th component is given by

[TABLE]

and similarly for the other components of the velocity. Here $L$ is the stiffness matrix and $F^{x}$ is a vector that accounts for the pressure and force terms. Let $f^{x}$ and $\mathbf{n}^{x}$ denote the $x$-th component of the force field and the $x$-th component of the normal unit vector to the boundary respectively. Then, the entries of $L$ and $F^{x}$ are given by

[TABLE]

Note that we integrate by parts for the pressure term.

Remark 5.2.1** (About the DMP property).**

When the smoothness indicator $\beta^{x}$ is constructed based on $U^{n+1}$, the mass matrix is lumped such that $m_{i}=\int_{\Omega}w_{i}d\mathbf{x}>0$ and the corresponding $D_{ij}$ matrix is applied to (15), the DMP property is guaranteed, see e.g. [1]. By using $U^{*}$ instead, the DMP property is likely lost. In addition, the viscous, pressure and force terms require extra care. Also note that, to avoid loss of accuracy, we do not lump the mass matrix. Nevertheless, as mentioned before, we do not aim to satisfy a DMP. Instead, our objective is to add robust and reliable artificial viscosity to dissipate small spatial scales not resolved by the computational mesh in order to obtain a well-behaved method for any mesh and any refinement level while nevertheless converging rapidly.

Remark 5.2.2** (About accuracy).**

It is known that the discrete operator $D_{ij}^{\text{Low}}$ decreases the order of convergence to first order even for smooth solutions, see for instance [28]. Second order convergence is recovered except around local extrema when the smoothness indicator $\beta$ is used. The degeneracy of accuracy around local extrema occurs since $\beta=0$ at local extrema. To overcome this problem, one can consider an extra smoothness sensor based on second derivatives, which is suitable for the velocity space $X_{h}^{2}$, see [31]. The objective is to identify smooth local extrema to deactivate the artificial dissipation. In [29] this idea is employed for different scalar equations. By doing this, the authors demonstrate that the high-order accuracy of the underlying method is recovered. We do not explore this idea further in this work.

5.3 Surface tension

In this section, we summarize the work by [21] to incorporate surface tension effects. Given $\Gamma\subset\Omega$ to be the interface between the two fluids, the goal is to impose $[\mathbf{u}]|_{\Gamma}=0$ and $-[-p{\bf I}+2\mu\varepsilon(\mathbf{u})]|_{\Gamma}\cdot\mathbf{n}^{\Gamma}=\sigma\kappa\mathbf{n}^{\Gamma}$, where $\mathbf{n}^{\Gamma}$ is the unit normal vector to the interface, $\sigma$ is the surface tension coefficient and $\kappa$ is the curvature of the interface. These conditions can be imposed by incorporating

[TABLE]

into the left hand side of (13). Let $\underline{\nabla}z(\mathbf{x})=\nabla z(\mathbf{x})-[\mathbf{n}^{\Gamma}\cdot\nabla z(x)]\mathbf{n}^{\Gamma},\leavevmode\nobreak\ \forall\mathbf{x}\in\Gamma$ be the tangential gradient with respect to $\Gamma$ and $\underline{\Delta}z(\mathbf{x})=\underline{\nabla}\cdot\underline{\nabla}z(\mathbf{x}),\leavevmode\nobreak\ \forall\mathbf{x}\in\Gamma$ be the Laplace-Beltrami operator. Using the fact that $\kappa\mathbf{n}^{\Gamma}=\underline{\Delta}\mathbf{x}|_{\Gamma}$, see for instance [15], the weak form of (18), after multiplying $\mathbf{f}_{\text{st}}$ by a test function $\mathbf{w}(\mathbf{x})$ and integrating by parts, is given by

[TABLE]

where $\text{Id}_{\Gamma}(\mathbf{x})$ is the identity map on $\Gamma$. The boundary integral vanishes for the cases of interest in this work. Once the force integral related to the surface tension is defined, the authors in [21] propose to treat (19) semi-implicitly based on $\text{Id}_{\Gamma}^{n+1}\approx\text{Id}_{\Gamma}^{n}+\Delta t^{n+1}\mathbf{u}^{n+1}$ and to approximate the surface integrals by volume integrals via the regularized delta function of the level set. Doing this yields

[TABLE]

Assuming the acute angle condition is satisfied, see for instance [9], it is remarked in [21] that since the second term in (20) is positive, it contributes to the stability properties of the momentum equations. In other words, it adds viscosity to the velocity at the interface $\Gamma$. See more details about this discretization in the next section.

5.4 Full discretization

Let $\mathbf{u}_{h},\mathbf{w}\in[X_{h}^{2}]^{d}$ and $p_{h},\delta\psi_{h},\theta\in X_{h}^{\text{1}}$. The full discretization of (14) is given by

[TABLE]

Remark 5.4.1** (Velocity correction).**

After the pressure update, it is possible to correct the velocity field to obtain a weakly divergence free velocity field. To do this we can rearrange the terms in (21c) to obtain

[TABLE]

*and redefine the velocity such that $\mathbf{u}_{h}^{n+1}\rightarrow\mathbf{u}_{h}^{n+1}-\frac{1}{a}\nabla\delta\psi_{h}^{n+1}$ in the interior of the domain. By doing this it is clear that *

[TABLE]

It is remarked in [18, §3.5], and references therein, that not doing this correction does not affect the accuracy properties of the method. However, as explained in §4, the level set model assumes the velocity to be (at least weakly) divergence free. Therefore, we correct the velocity field.

6 Numerical experiments

A total of eight test problems were selected to demonstrate the behavior of the method we propose. The results were compared qualitatively and quantitatively versus other results in the literature and versus experimental measurements. These two-phase flow experiments represent a wide range of applications ranging from high-viscosity buckling fluids to free surface flows around obstacles. Our main objective is to test the robustness of the method for a wide range of physical parameters, surface tension, different type of boundary conditions, external forces and arbitrary mesh refinements. In particular, we set the numerical parameters to be the same for all problems; in our opinion, this demonstrates the robustness of the numerical method and suitability of the computational model for engineering applications.

Parameters

The physical parameters are part of the definition of the problem. They are given by the density ($\rho$) and viscosity ($\mu$) of each phase, the magnitude of the gravity ($g$) and the surface tension coefficient ($\sigma$). Unless otherwise stated, we consider the metric system and set

[TABLE]

where $W$ and $A$ refer to water and air, respectively. In the rest of this work we omit the units. The parameter $\tilde{\lambda}$ in the conservative level set method is the main numerical parameter, see remark 4.0.2. In all simulations we use $\tilde{\lambda}=10$. A less important parameter is $\epsilon$, which is used to regularize Heaviside and Dirac delta functions, see §4. During the redistancing process we set $\epsilon=\frac{1}{3}$ and elsewhere we use $\epsilon=\frac{3}{2}$. In all problems we let $h_{e}$ denote a characteristic element size (of the unstructured mesh) and report the number of elements.

Commonly used initial and boundary conditions

In different problems we start ‘at rest’; i.e., the initial velocity and pressure fields are set to zero. Likewise, we use similar and common boundary conditions in different problems. We refer to them as slip, non-slip, open top and inflow boundary conditions. These common boundary conditions for the velocity are shown in table 4, where $(\cdot)_{\tau}$ denotes the tangential component of a given vector field. The level set does not require boundary conditions when the boundary is set to slip or non-slip. When the boundary is open and $\mathbf{u}\cdot\mathbf{n}<0$ we set $S_{\epsilon}(\phi)=1$; i.e., we let only air through the open boundary. Finally, when the boundary is of inflow type we set either $S_{\epsilon}(\phi)=0$ or $S_{\epsilon}(\phi)=1$ to let water or air respectively into the domain.

Computational framework

This work was developed using Proteus (https://proteustoolkit.org), a toolkit for computational methods and simulation released under the MIT open source license with source available at https://github.com/erdc/proteus. The numerical linear algebra is handled by PETSc, see [3, 4, 2], through petsc4py [12]. Inside Proteus 1.5.1, we created a set of files to facilitate the definition of the different problems. This framework is not part of Proteus; nevertheless, to facilitate its use, we provide a release with this framework incorporated, see release 1.5.1-mp-r1 (https://github.com/erdc/proteus/releases/tag/1.5.1-mp-r1). Our aim is to provide a solid, robust and easy to use open-source computational framework for users interested in two-phase flows. We do not assume the users have any knowledge on finite elements nor expect the need to adjust the numerical parameters, except for potentially $\tilde{\lambda}$. Instead, we expect the user to provide information only about the problem; in particular, the domain, initial and boundary conditions and physical parameters. We encourage the interested reader to install the Proteus release and run different problems.

6.1 Two-dimensional experiments

6.1.1 Rising bubbles

We consider first a bubble of a light fluid at rest immersed in a heavier fluid. Due to the action of gravity, the bubble rises through the heavier fluid. In this problem, the effect of surface tension is critical to maintain the correct shape of the bubble and to obtain the correct position and rising speed. We consider the two common test cases in [22]. The domain of interest is $\Omega=(0,1)\times(0,2)$. We impose non-slip boundary conditions on the bottom and top boundaries and slip boundary conditions on the left and right boundaries. The initial condition for level set is initialized as explained in §4.1 considering an interface given by $\Gamma=\{(x,t)\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ r^{2}=(x-0.5)^{2}+(y-0.5)^{2},r=0.25\}$. The material parameters for the different test cases are shown in table 5. We consider structured meshes with refinement levels $h_{e}=\frac{1}{80},\frac{1}{160}$ and $h_{e}=\frac{1}{320}$, which corresponds to 6480, 25760 and 102720 elements respectively. The zero contour plot of the level set and the evolution of the center of mass $c_{y}=\frac{1}{|\Omega_{B}|}\int_{\Omega_{B}}yd\mathbf{x}$ and rising velocity $r_{v}=\frac{1}{|\Omega_{B}|}\int_{\Omega_{B}}\mathbf{u}_{y}d\mathbf{x}$ are shown in figure 3. Here $\Omega_{B}=\{(x,y)\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ \phi(x,t)\geq 0\}$ denotes the domain occupied by the bubble.

6.1.2 Dam break problem

Here we consider a common two dimensional dam break problem. In [34], the authors considered a similar problem and performed a series of experiments placing pressure and water height gauges at different locations. Based on these results, different numerical studies have been performed to validate numerical methods and codes. In this work, we consider the setup by [11] and reproduce some of their results. The domain of interest is $\Omega=(0,3.22)\times(0,1.8)$. An initial column of water at rest is located in $W=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ x\leq 1.2,\leavevmode\nobreak\ y\leq 0.6\}$. The rest of the domain ($A=\Omega\setminus W$) is filled with air. At time $t=0$ the column of water starts to fall because of the action of gravity. For this problem, we follow [11] and consider non-slip boundary conditions everywhere except for the top boundary, which is left open. We report the results of only one pressure gauge located at $P_{1}=(3.22,0.12)$ and two water height gauges located at $H_{1}=(2.228,y)$ and $H_{2}=(2.724,y),\leavevmode\nobreak\ \forall y\in[0,1.8]$. Let $d=3.22$ denote the length of the domain. We consider three unstructured meshes with $h_{e}=0.01d,\leavevmode\nobreak\ 0.00625d$ and $h_{e}=0.003125d$, which correspond to 17752, 45188 and 181033 triangular elements respectively. In figure 4 we compare the numerical results (for the different refinements) versus the experimental measurements for all the gauges. Finally, in figure 5 we consider the same times as [11, figure 17] and show the water-air interface. Note that only in this figure, to facilitate comparisons, we follow the reference and scale the $x$-axis as $\tilde{x}=\frac{x}{H}$ and the time as $\tilde{t}=t\sqrt{g/H}$, where $H=0.6$ is the initial height of the column of water.

6.1.3 Buckling flow

The fluid buckling phenomenon occurs in situations involving thin streams of highly viscous flow that encounter a boundary or plate. Depending on the Reynold’s number and the cross sectional geometry of the viscous jet, the fluid may fold, coil, or oscillate antiperiodically. Here, the method’s ability to simulate buckling flow is evaluated against results from previous literature, see e.g., [7, 44, 43, 8]. The domain $\Omega=(0,1)^{2}$ is initially occupied by the air phase, which is at rest. A small inlet on the top boundary of the box, defined by $I=\{\mathbf{x}\in\bar{\Omega}\leavevmode\nobreak\ |\leavevmode\nobreak\ |x-0.5|<0.05\}$, is the inflow boundary, where the velocity is strongly set to $\mathbf{u}=(0,-1)$. The rest of the top boundary is left open. The non-slip boundary condition is applied in the bottom, left and right boundaries. The material properties are shown in table 6. We use an unstructured mesh with a maximum element size of $h_{e}=5\times 10^{-3}$, which correspond to 126,646 triangular elements. In figure 6 we plot the heavy fluid phase for different times.

6.1.4 Filling tank

In this section we consider the problem of filling an empty tank with a water-like fluid. For this problem, the domain of interest is given by $\Omega=(0,0.4)^{2}$. The initial data consists of water in $W=\{(x,y)\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ x<0.01,\leavevmode\nobreak\ |y-0.325|<0.025\}$ and air in $A=\Omega\setminus W$. Both phases start at rest. The boundary $I=\{(x,y)\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ x=0,\leavevmode\nobreak\ |y-0.325|<0.025\}$ is set to inflow boundary. At $I$ we set strongly $\mathbf{u}=(0.25,0)$. The rest of the left boundary is defined as slip boundary, the bottom and the right boundaries are non-slip and the top is left open. We consider two sets of material parameters. First we follow [19, §11.3] and reproduce qualitatively their results. In this case the material parameters are

[TABLE]

Then we consider more realistic parameters given by (23a), without surface tension effects; i.e., $\sigma=0$. We use an unstructured mesh with maximum element size of $h_{e}=2\times 10^{-3}$, which corresponds to 126,888 triangular elements. The water phase at different times are shown in figure 7. Note the difference in the behavior of the solution due to different gravity and viscosity parameters.

6.2 Three-dimensional experiments

6.2.1 Dam break problem with obstacle

The first three dimensional problem that we consider is a dam break with an obstacle located downstream. The computational domain is $\Omega=(0,3.22)\times(0,1)\times(0,1)\setminus O_{1}$, where $O_{1}=[2.39,2.55]\times[0.3,0.7]\times[0,0.16]$ represents the obstacle. An initial column of water at rest is located in the domain $W=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ x\leq 1.22,\leavevmode\nobreak\ z\leq 0.55\}$. The rest of the domain; i.e, $A=\Omega\setminus W$, is filled with air. At $t=0$ the column of water starts to fall because of the action of gravity. Non-slip boundary conditions are imposed everywhere except for the top, which is left open. Experiments for this problem were performed by the Maritime Research Institute Netherlands (MARIN), see [26, 25]. In particular, pressure and water height gauges were located in different points of interest. In this work we report the results for only four pressure gauges located at $P_{1}=(2.389,0.526,0.025)$, $P_{3}=(2.389,0.526,0.099)$, $P_{5}=(2.414,0.474,0.165)$ and $P_{7}=(2.487,0.474,0.165)$ and four water height gauges located at $H_{1}=(0.582,0.5,z),\leavevmode\nobreak\ H_{2}=(1.732,0.5,z),\leavevmode\nobreak\ H_{3}=(2.228,0.5,z)$ and $H_{4}=(2.724,0.5,z),\leavevmode\nobreak\ \forall z\in[0,1]$. See figure 8. We consider an unstructured mesh with characteristic size $h_{e}=0.025$, which corresponds to 2,552,090 tetrahedral elements. In figure 9, we compare the numerical results versus the experimental measurements for the pressure and water height gauges and in figure 10 we show the water phase at different times.

6.2.2 Moses flow

In this simulation we consider the domain $\Omega=(0,3.22)\times(0,1)\times(0,1)\setminus O_{1}\cup O_{2}$, where $O_{1}=[0.6635,0.8245]\times[0.2985,0.7015]\times[0,0.35]$ and $O_{2}=[2.3955,2.5565]\times[0.2985,0.7015]\times[0,0.35]$ represent two obstacles located at both ends of the domain. Initially, the tank is partially filled with water and air at rest in the domains $W=\{\mathbf{x}\in\Omega\leavevmode\nobreak\ |\leavevmode\nobreak\ z\leq 0.3\}$ and $A=\Omega\setminus W$, respectively. We apply an external force designed to split the water at the center of the tank. Our aim with this problem is to test the robustness of the method with respect to external (and aggressive) forces. The force is

[TABLE]

Note that the external force is scaled by the fluid density. All boundaries are considered to be non-slip except the top which is left open. The grid is unstructured with element size $h_{e}=0.025$, which corresponds to 2,522,647 tetrahedral elements. In figure 11 we show the water phase at different times. We remark that the solution appears to tend to a steady state.

7 Conclusions

In this work we presented a numerical model of incompressible, immiscible two-phase flow. The general algorithm is driven by an operator splitting scheme in which a conservative level set method is first solved for a given velocity field to evolve the interface. Afterwards, the interface (at the new location) is employed to obtain an updated velocity field.

We extended the conservative monolithic level set method presented in [38] to make it more robust and to allow simpler initial conditions. The original monolithic scheme contains a second-order term that penalizes deviations from the distance function, which was derived from parabolic and elliptic redistancing, see [10, 6]. In our experience with two-phase flow simulations using the method in [38], we found it is not always possible to force the redistancing to emanate from the interface, which results in incorrect solutions (of the viscous Eikonal equation [32]). As this was not an issue with the earlier multistage conservative level set method presented in [24] based on the classical Eikonal equation, we propose performing an extra pre-processing step to induce the redistancing to be propagated from the interface, effectively providing a better initial guess to the original monolithic scheme. This extra step is done by solving the Eikonal equation directly. It is important to remark that, thanks to the penalization embedded into the conservative level set, solving the Eikonal equation in our case is not a computationally expensive task. This is true since the initial guess in the non-linear iterative process already contains most of the features of the solution. We see this pre-processing step as a small but necessary improvement for the initial condition at each time step and consider a more robust monolithic scheme as future work.

We use a second order projection scheme proposed for for variable-density incompressible flows in [17] extended to solve the Navier-Stokes equations with variable density, viscosity, and time steps. For the applications of consideration in this work, extra artificial viscosity is commonly needed to stabilize the momentum equations at high Reynolds number. Special care is required when some of these techniques are used with unstructured meshes. In this work, we extend recent ideas on solving hyperbolic equations via Discrete Maximum Principle (DMP) preserving continuous Galerkin finite elements. These methods are algebraic and suitable for both coarse and highly refined unstructured meshes. We relax the DMP preserving features of these schemes, which allows simplification of the stabilization; however, in our numerical experiments we always obtained well behaved solutions. We remark that this stabilization technique is free of tunable parameters and capable of high-order accuracy with suitable changes to the underlying approximation spaces.

Once the individual methods are described, we present an extensive set of numerical examples in two- and three-dimensions. In all these problems, we use the same numerical parameters which, in our opinion, demonstrates the robustness of the method. Finally, we provide an open source and freely downloadable computational framework to test and try the numerical examples presented in this work and others that the reader might be interested in.

8 Acknowledgments

The work of Manuel Quezada de Luna was supported primarily by an appointment to the Postgraduate Research Participation Program at the U.S. Army Engineer Research and Development Center, Coastal and Hydraulics Laboratory (ERDC-CHL) administrated by the Oak Ridge Institute for Science and Education through an inter-agency agreement between the U.S. Department of Energy and ERDC. Haydel Collins and Chris Kees were supported by the ERDC University program and the ERDC Future Investment Fund. Permission was granted by the Chief of Engineers, US Army Corps of Engineers, to publish this information.