High-Order Finite-Element Framework for the Efficient Simulation of Multifluid Flows

Abstract

In this paper, we present a comprehensive framework for the simulation of Multifluid flows based on the implicit level-set representation of interfaces and on an efficient solving strategy of the Navier-Stokes equations. The mathematical framework relies on a modular coupling approach between the level-set advection and the fluid equations. The space discretization is performed with possibly high-order stable finite elements while the time discretization features implicit Backward Differentation Formulae of arbitrary order. This framework has been implemented within the Feel++ library, and features seamless distributed parallelism with fast assembly procedures for the algebraic systems and efficient preconditioning strategies for their resolution. We also present simulation results for a three-dimensional Multifluid benchmark, and highlight the importance of using high-order finite elements for the level-set discretization for problems involving the geometry of the interfaces, such as the curvature or its derivatives.

1. Introduction

Understanding and predicting the dynamics of systems consisting of multiple immiscible fluids in contact is a great challenge for numerical computations, as they involve bulk coupling of the fluids—related to the long-range features of the Navier-Stokes equations—and possible strong surface effects.

Such interfaces arise in a variety of applications in physics and engineering. For instance, the dynamics of drops immersed in another fluid [1], shocks in high-speed flows of aerospace interest [2], vortex sheets in wakes [3] or at the boundaries between immiscible fluids in engineering applications [4]. A lot of effort has been thus recently put in the development of efficient numerical methods to solve these strongly coupled fluid problems while tracking the interfaces to accurately resolve the surface effects.

One of the main difficulties of these Multifluid simulations is to keep track of the interfaces between the fluids as they evolve in time, as well as to characterize accurately their geometry. Two main approaches are typically used to account for the interface changes. In the first one, the interface is tracked explicitely with a moving mesh. This large class of “front-tracking” methods—which embrace, for example, the Arbitrary Lagrangian-Eulerian (ALE [5]), the Fat Boundary Method (FBM, [6])—feature a very accurate description of the contact surface, but at high computational costs. The second approach uses only a fixed mesh, and represents the interface implicitly with some additional field. For example, this is the approach used by the Volume Of Fluid (VOF [7]), the Phase-Field method [8,9] and the Level-Set method (LS [10]). These methods in general provide easier coupling formulations and smaller algebraic problems but with less accuracy regarding the description of the interface.

In this work, we present a level-set based framework for the simulation of generic Multifluid flows. It features efficient solving strategies for the fluid–incompressible Navier-Stokes (or Stokes)— equations and the level-set advection. This framework is implemented using the open-source Feel++ library—Finite Element Embedded Library in C++ (v0.104.1, Feel++ Consortium, Strasbourg, France)—[11,12,13] within the FluidMechanics, AdvectionDiffusionReaction, LevelSet and Multifluid toolboxes.

These toolboxes expose user-friendly interfaces and allow versatile parametrization of the problems and solving strategies while managing the parallel assembly of the finite-element algebraic systems and their resolution seamlessly.

We also present a 3D benchmark for our toolboxes. This numerical experiment is an extension of the classic rise of a drop in a viscous fluid in 2D. We use two different setups and compare our results with other approaches to validate our framework.

We stress that all the implementations and testcases presented below are freely available [14,15] and can be used for a large class of fluid-structure interaction and suspension problems.

2. Simulating Multifluid Flows

We consider a system consisiting of several non-miscible fluids with different physical properties occupying some domain $\mathsf{\Omega }$. We denote ${\mathsf{\Omega }}_i\subset \mathsf{\Omega }$ ($\mathsf{\Omega }={\bigcup }_i{\mathsf{\Omega }}_i$) the domain occupied by the ith fluid. We want to study the dynamics of such a system, which means both solving for the motion of the fluids and tracking each fluid subdomain as it evolves in time.

2.1. Fluid Equations

In this work, we assume that all the considered fluids are incompressible, and thus obey the incompressible Navier-Stokes equations

$$\begin{array}{cc}\hfill \rho [{\partial }_t\overrightarrow{u}+(\overrightarrow{u}·\nabla )\overrightarrow{u}]-\nabla ·\sigma (\overrightarrow{u},p)& ={\overrightarrow{f}}_b+{\overrightarrow{f}}_s\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \nabla ·\overrightarrow{u}& =0\hfill \end{array}$$

(2)

where $\overrightarrow{u}$ and p are respectively the velocity and pressure fields, $\rho$ is the fluid density, $\sigma$ is the fluid stress tensor and ${\overrightarrow{f}}_b$, ${\overrightarrow{f}}_s$ are the bulk and surface forces respectively. We consider in the following Newtonian fluids, so that

$$\sigma (\overrightarrow{u},p)=-p\mathbf{I}+2\mu \mathit{D}(\overrightarrow{u})$$

(3)

with $\mu$ the fluid viscosity and $\mathit{D}(\overrightarrow{u})=\frac{1}{2}(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T)$ is the rate of strain tensor of the fluid.

Note that these equations are satisfied for each fluid independently. In a Multifluid system, they must be supplemented with boundary conditions for the velocity and pressure at both the domain boundary and at the interface between two fluids. In this work, we assume that there is no slip at the interfaces, so that the velocity and pressure fields are continuous at all the interfaces. The boundary conditions at the domain boundaries can be Dirichlet or Neumann:

$$\begin{array}{ccc}\hfill \overrightarrow{u}& ={\overrightarrow{g}}_D\hfill & \hfill \mathrm{on}\partial {\mathsf{\Omega }}_D\\ \hfill \sigma (\overrightarrow{u},p)\overrightarrow{n}& ={\overrightarrow{g}}_N\hfill & \hfill \mathrm{on}\partial {\mathsf{\Omega }}_N\end{array}$$

(4)

2.2. Interface Tracking: The Level-Set Method

Let us consider two disjoint—not necessarily connected—fluid domains ${\mathsf{\Omega }}_i$, ${\mathsf{\Omega }}_j\subset \mathsf{\Omega }$ such that $\partial {\mathsf{\Omega }}_i\cap \partial {\mathsf{\Omega }}_j\ne \varnothing$, and denote ${\mathsf{\Gamma }}_{ij}$ the interface between them. In order to efficiently track this interface ${\mathsf{\Gamma }}_{ij}(t)$ as it evolves in time, we use the level set method [10,16,17] which provides a natural way to compute geometrical properties of the interface and to handle possible topological changes in a completely Eulerian framework.

It features a Lipschitz continuous scalar function $\varphi$ (the level set function) defined on the whole domain. This function is arbitrarily chosen to be positive in ${\mathsf{\Omega }}_i$, negative in ${\mathsf{\Omega }}_j$ and zero on ${\mathsf{\Gamma }}_{ij}$, so that the interface is implicitely represented by the 0-level of $\varphi$.

As ${\mathsf{\Omega }}_i$ and ${\mathsf{\Omega }}_j$ evolve following the Navier-Stokes Equation (2) dynamics, ${\mathsf{\Gamma }}_{ij}$ gets transported by the velocity field $\mathit{u}$ (uniquely defined at the interface by assumption and construction) and therefore obeys the advection equation:

$$\frac{\partial \varphi }{\partial t}+\mathit{u}·\nabla \varphi =0.$$

(5)

We choose as level set function the signed distance to the interface

$$\varphi (\mathit{x})=\left\{\begin{array}{cc}dist(\mathit{x},{\mathsf{\Gamma }}_{ij})& \mathit{x}\in {\mathsf{\Omega }}_i,\hfill \\ 0& \mathit{x}\in {\mathsf{\Gamma }}_{ij},\hfill \\ -dist(\mathit{x},{\mathsf{\Gamma }}_{ij})& \mathit{x}\in {\mathsf{\Omega }}_j,\hfill \end{array}\right.$$

(6)

as the intrinsic property $|\nabla \varphi |=1$—which is characteristic of the Euclidean distance for all points $\overrightarrow{x}$ such that ${\arg \min }_{{\overrightarrow{x}}_0\in {\mathsf{\Gamma }}_{ij}}|\overrightarrow{x}-\overrightarrow{x_0}|$ is a unit set–eases the numerical resolution of the advection equation by ensuring that no steep gradients are present. Furthermore, the regularity of the distance function allows us to use $\varphi$ as a support for the regularized interface Dirac and Heaviside functions (see below).

However, the advection Equation (5) does not preserve the property $|\nabla \varphi |=1$, and can lead to the accumulation or rarefaction of $\varphi$ iso-levels depending on the velocity. It is therefore necessary to enforce the condition $|\nabla \varphi |=1$ externally, and to reset $\varphi (t)$ to a distance function without moving the interface, Refs. [18,19,20]. To redistantiate $\varphi (t)$, we can either solve a Hamilton-Jacobi equation, which “transports” iteratively the isolines of $\varphi$ depending on the sign of $|\nabla \varphi |-1$ until the steady “distance” solution is attained, or use the fast marching method, which computes the actual distance to the interface for each degree of freedom using an upwind $\mathcal{O}(NlogN)$ scheme starting from the interface. We provide further details about our numerical implementation below (c.f. Section 3.5.3).

As mentioned above, the signed distance to the interface $\varphi$ allows us to easily define the regularized interface-related Dirac and Heaviside functions which can be used to compute integrals on the interface [21,22] or on a fluid subdomain

$${\delta }_{\epsilon }(\varphi )=\left\{\begin{array}{cc}0,& \varphi \le -\epsilon ,\\ \frac{1}{2\epsilon }\left[1+cos\left(\frac{\pi \varphi }{\epsilon }\right)\right],& -\epsilon \le \varphi \le \epsilon ,\\ 0,& \varphi \ge \epsilon .\end{array}\right.$$

(7)

$$H_{\epsilon }(\varphi )=\left\{\begin{array}{cc}0,& \varphi \le -\epsilon ,\\ \frac{1}{2}\left[1+\frac{\varphi }{\epsilon }+\frac{sin(\frac{\pi \varphi }{\epsilon })}{\pi }\right],& -\epsilon \le \varphi \le \epsilon ,\\ 1,& \varphi \ge \epsilon .\end{array}\right.$$

(8)

where $\epsilon$ is a parameter controlling the “numerical thickness” of the interface. We note in these definitions that enforcing $|\nabla \varphi (t)|=1$ is critical to ensure that the interfacial support of ${\delta }_{\epsilon }$ and $H_{\epsilon }$ is kept constant and larger than the mesh size. Typically, we choose $\epsilon$ to scale with the mesh size as $\epsilon \sim 1.5h$, where h the average mesh size.

The Heaviside function is used to define physical quantities which have different values on each subdomain. For example, we define the density and viscosity of two-fluid flows as $\rho ={\rho }_j+({\rho }_i-{\rho }_j)H_{\epsilon }(\varphi )$ and $\mu ={\mu }_j+({\mu }_i-{\mu }_j)H_{\epsilon }(\varphi )$. The delta function allows to define quantities on the interface using non-conform meshes. In particular in the variational formulations we can replace integrals over the interface $\mathsf{\Gamma }$ with integrals over the entire domain $\mathsf{\Omega }$ using the regularized delta function: if $\varphi$ is a signed distance function (i.e., $|\nabla \varphi |=1$), we have ${\int }_{\mathsf{\Gamma }}1\simeq {\int }_{\mathsf{\Omega }}{\delta }_{\epsilon }(\varphi )$.

In the following, we consider without loss of generality only the two-fluid case, with only one level set function tracking the interface between fluids 1 and 2.

2.3. Finite Element Formulation

We use finite elements methods to solve Equations (2) and (5), and work with a continous Galerkin formulation. As mentioned above, this can be done by smoothing the discontinuities of the fluid parameters (e.g., the fluid density and viscosity) at the interfaces using the regularized Dirac and Heaviside functions Equations (7) and (8). The continuity of the velocity and pressure fields is then imposed strongly in the formulation, and we can work with function spaces defined on the whole domain $\mathsf{\Omega }$.

We introduce $L^2(\mathsf{\Omega })$ the usual function space of square integrable functions, $H^1(\mathsf{\Omega })$ the function space of square integrable functions which gradients are also square integrable and the vectorial Sobolev spaces ${\mathbf{H}}^1(\mathsf{\Omega })=\{\overrightarrow{v}\in {\left[H^1(\mathsf{\Omega })\right]}^d\}$, ${\mathbf{H}}_{g_D}^1(\mathsf{\Omega })=\{\overrightarrow{v}\in {\mathbf{H}}^1(\mathsf{\Omega }),\overrightarrow{v}{|}_{\partial {\mathsf{\Omega }}_D}={\overrightarrow{g}}_D\}$ and ${\mathbf{H}}_{0_D}^1(\mathsf{\Omega })=\{\overrightarrow{v}\in {\mathbf{H}}^1(\mathsf{\Omega }),\overrightarrow{v}{|}_{\partial {\mathsf{\Omega }}_D}=\overrightarrow{0}\}$.

The variational formulation of the two-fluid coupling problem Equations (2) and (5) then reads Find $(\overrightarrow{u},p,\varphi )\in {\mathbf{H}}_{g_D}^1(\mathsf{\Omega })\times L^2(\mathsf{\Omega })\times H^1(\mathsf{\Omega })$ s.t. $\forall (\overrightarrow{v},q,\psi )\in {\mathbf{H}}_{0_D}^1(\mathsf{\Omega })\times L^2(\mathsf{\Omega })\times H^1(\mathsf{\Omega })$,

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\rho (\varphi )\left[{\partial }_t\overrightarrow{u}+(\overrightarrow{u}·\nabla )\overrightarrow{u}\right]·\overrightarrow{v}+{\int }_{\mathsf{\Omega }}\sigma (\overrightarrow{u},p;\mu (\varphi )):\nabla \overrightarrow{v}& ={\int }_{\mathsf{\Omega }}{\overrightarrow{f}}_b·\overrightarrow{v}+{\int }_{\mathsf{\Omega }}{\overrightarrow{f}}_s(\varphi )·\overrightarrow{v}+{\int }_{\partial {\mathsf{\Omega }}_N}{\overrightarrow{g}}_N·\overrightarrow{v}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\left(\nabla ·\overrightarrow{u}\right)q& =0\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\left({\partial }_t\varphi +\overrightarrow{u}·\nabla \varphi \right)\psi & =0\hfill \end{array}$$

(11)

where we have integrated by parts the stress tensor term and used the Neumann boundary condition in Equation (4). The surface force term ${\int }_{\mathsf{\Omega }}{\overrightarrow{f}}_s·\overrightarrow{v}$ is here generically written as a bulk integral, the surfacic aspect being hidden in the force expression ${\overrightarrow{f}}_s(\varphi )$ which in general contains some Dirac distribution $\delta (\varphi )$. In our case, the surface forces actually account for the interfacial forces between the two fluids (i.e., the surface tension), so that the surface integral can be evaluated as a bulk integral using the regularized level set delta function.

3. Numerical Setup

3.1. The Feel++ Toolboxes

The numerical implementation is performed using the Feel++—finite element C++ library [11,12,13], which provides a comprehensive finite-element framework for performing HPC resolution of generic problems. Feel++ allows to use a very wide range of Galerkin methods and advanced numerical methods such as domain decomposition methods including mortar and three fields methods, fictitious domain methods or certified reduced basis. It features a very expressive domain specific embedded language, and powerful localization techniques used for interpolation and mixed spaces integration. It also provides automatic and efficient assembly of the algebraic systems, as well as seamless parallelization. The Feel++ library has been used in many contexts ranging from the development and numerical verification of innovative mathematical methods to the simulation of large multi-physics systems [23,24,25].

The Feel++ mathematical kernel can be used to solve linear and nonlinear partial differential equation using arbitrary order Galerkin methods (fem, sem, cg, dg, hdg, crb). It naturally handles 1d, 2d, 3d manifolds using simplices and hypercubes meshes [11,12,13,26] and provides in particular:

i

a polynomial library allowing for a wide range of polynomial expansions including $H_{\mathrm{div}}$ and $H_{\mathrm{curl}}$ elements,

ii

a lightweight interface to Boost.UBlas, Eigen3 and PETSc/SLEPc

iii

a scalable in-house solution strategy, in particular with specialized preconditioners which can easily be tuned from the configuration files

iv

a natural language for Galerkin methods allowing the definition of function spaces, (bi)linear forms, operators, functionals and integrals,

v

a framework that allows user codes to scale seamlessly from single core computation to thousands of cores and enables hybrid computing.

Feel++ provides also an environment for modeling and solving various kinds of scientific and engineering problems. The framework, implemented in several Feel++ toolboxes, provides a language to describe these models. Based on JSON and INI(we use the .cfg extension) file formats, it is possible to configure and simulate a large class of models by defining the relevant physical quantities—such as the material properties, the boundary conditions, the sources, the couplings and the solvers. In this paper, we focus on the Multifluid toolbox which allows to setup all the necessary ingredients for a Multifluid flow simulation. This toolbox, which mainly manages the coupling between the fluids present in the system, is built on others toolboxes inclined toward monophysics problem resolution, namely the FluidMechanics, LevelSet and AdvectionDiffusionReaction toolboxes.

3.2. Fluid-Interface Coupling

The coupling between the fluid and the level-set in Equations (9)–(11) is highly nonlinear and solving these equations monolithically would require specific nonlinear solvers adapted to each particular coupling force ${\overrightarrow{f}}_s(\varphi )$. In order to ease the implementation and development processes, and to benefit from efficient solving strategies, we choose an explicit—non-monolithic—coupling in our numerical approach.

At each time step, we first solve the fluid equations using the physical parameters and the surface forces computed from the last-step level-set function. We then use the obtained fluid velocity to advect the level-set and get the new interface position. We can then apply some redistantiation procedure depending on the chosen strategy before proceeding to the next iteration.

The fluid-level-set coupling algorithm then writes

Note that the successive resolution of the fluid and level-set equations can also be iterated within one time step, until a fix point of the system of equations is reached. In practice however, for reasonably small time steps, the fix-point solution is already obtained after the first iteration.

3.3. Space/Time Discretization

We introduce ${\mathcal{T}}_h\equiv \left\{K_e,1\le e\le N_{elt}\right\}$ a compatible tessellation of the domain $\mathsf{\Omega }$, and denote ${\mathsf{\Omega }}_h={\bigcup }_{e=1}^{N_{elt}}{K}_e$ the discrete–unstructured–mesh associated with average mesh size h. We work within the continuous Galerkin variational formulation framework, and use Lagrange finite elements to spatially discretize and solve the equations governing the evolutions of the fluid and the level-set. We thus introduce ${\mathcal{P}}_h^k\equiv {\mathcal{P}}_h^k\left({\mathsf{\Omega }}_h\right)$ the discrete (h-dependent) finite element space spanned by Lagrange polynomials of order k.

Then, from the time interval $[0,T]$, we select $M+1$ equidistributed times: $\left\{t_i=i\delta t,0\le i\le M\right\}$, where $\delta t$ is the time step. In the Navier-Stokes and level-set advection equations, we apply a fully implicit time discretization by using a backward differentiation formula of arbitrary order N named ${\mathrm{BDF}}_N$. Let $\varphi$ be a function and denote ${\varphi }^{(i)}$ the function at time $t_i$. The time derivative of this function is then discretized as

$${\partial }_t\varphi =\frac{1}{\delta t}\left[{\alpha }_0{\varphi }^{(n+1)}-\sum _{i=1}^N{\alpha }_i{\varphi }^{(n+1-i)}\right]+\mathcal{O}(\delta t^{N+1})$$

(12)

where the ${\alpha }_i$ are ${\mathrm{BDF}}_N$ coefficients. In the numerical benchmarks reported below, we used the ${\mathrm{BDF}}_1$ and ${\mathrm{BDF}}_2$ schemes, which write respectively

$$\begin{array}{cc}\hfill {\mathrm{BDF}}_1:{\partial }_t\varphi & \approx \frac{1}{\delta t}\left[{\varphi }^{(n+1)}-{\varphi }^{(n)}\right]\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathrm{BDF}}_2:{\partial }_t\varphi & \approx \frac{1}{\delta t}\left[\frac{3}{2}{\varphi }^{(n+1)}-2{\varphi }^{(n)}+\frac{1}{2}{\varphi }^{(n-1)}\right]\hfill \end{array}$$

(14)

3.4. Solving the Incompressible Navier-Stokes Equations

The spatial discretization of the Navier-Stokes equations is handled via a inf-sup stable finite element (Taylor-Hood) ${[{\mathcal{P}}_h^{k+1}]}^d$, see e.g., [27]. We define $V_h$ and $Q_h$ the discrete function spaces where we search the velocity and the pressure solutions respectively. They are given by:

• $V_h=\left\{\overrightarrow{v}\in {\mathbf{H}}_{g_D}^1\cap {[{\mathcal{P}}_h^{k+1}]}^d\right\}$
• $Q_h=\left\{q\in {\mathcal{P}}_h^k\right\}$

We need also to introduce the test function space on the velocity with $W_h=\left\{\overrightarrow{v}\in {\mathbf{H}}_{0_D}^1({\mathsf{\Omega }}_h)\cap {[{\mathcal{P}}_h^{k+1}]}^d\right\}$. Hereafter we consider the case $k=1$. The discrete weak formulation associated to Equations (9) and (10) reads:

$$\begin{array}{c}\mathit{Find}({\overrightarrow{u}}_h^{(n+1)},p_h)\in v_h\times Q_h\mathit{such} \mathit{that}\forall ({\overrightarrow{v}}_h,p_h)\in W_h\times Q_h:\hfill \\ \begin{array}{cc}\hfill {\int }_{{\mathsf{\Omega }}_h}\rho (\varphi )& \left[\frac{{\alpha }_0}{\delta t}{\overrightarrow{u}}_h^{(n+1)}+\left({\overrightarrow{u}}_h^{(n+1)}·\nabla \right){\overrightarrow{u}}_h^{(n+1)}\right]·{\overrightarrow{v}}_h+{\int }_{{\mathsf{\Omega }}_h}\sigma ({\overrightarrow{u}}_h^{n+1},p_h^{n+1};\mu (\varphi )):\nabla {\overrightarrow{v}}_h+{\int }_{{\mathsf{\Omega }}_h}\left[\nabla ·{\overrightarrow{u}}_h^{(n+1)}\right]q_h\hfill \\ & ={\int }_{{\mathsf{\Omega }}_h}\sum _{i=1}^N\frac{{\alpha }_i}{\delta t}{\overrightarrow{u}}^{(n+1-k)}·{\overrightarrow{v}}_h+{\int }_{{\mathsf{\Omega }}_h}{\overrightarrow{f}}_b·{\overrightarrow{v}}_h+{\int }_{{\mathsf{\Omega }}_h}{\overrightarrow{f}}_s(\varphi )·{\overrightarrow{v}}_h+{\int }_{\partial {\mathsf{\Omega }}_N}{\overrightarrow{g}}_N·\overrightarrow{v}\hfill \end{array}\hfill \end{array}$$

The nonlinear problem presented by Equation (15) is solved monolithically with Newton’s method. From an algebraic point of view, at each nonlinear iteration, the following classical saddle-point system is inverted

$$\left(\begin{array}{cc}A& B\\ C& 0\end{array}\right)\left(\begin{array}{c}U\\ P\end{array}\right)=\left(\begin{array}{c}F\\ 0\end{array}\right)$$

(15)

where A corresponds to the velocity block and B, C to the velocity/pressure coupling. This system is solved with a gmres algorithm using the simple preconditioner, see [27].

In the benchmark presented in the next section, we have only Dirichlet boundary conditions defined on the whole boundary of the domain (i.e., $\partial {\mathsf{\Omega }}_N\equiv \varnothing$). In this case, the problem represented by the weak formulation (15) is not well-posed as the pressure is defined up to a constant. To solve this issue, our strategy consists to

• add the information to the Krylov subspace method (gmres) that the system has a null space, i.e., the pressure constant.
• rescale the pressure solution after each iteration of the Newton algorithm by imposing a mean pressure equal to 0.

Another solution available in the FluidMechanics toolbox is to add a Lagrange multiplier in order to impose the mean value of the pressure. The disadvantages of this approach are to increase the stencil of the matrix and to make its block structure much denser. This implies a reduction in the efficiency of the solver.

3.5. The Level-Set Framework

3.5.1. Level-Set Advection

As mentioned above, the transport of the level-set by the fluid is accounted by the advection equation Equation (11) in variational form. This equation is discretized in time using a backward differentiation formula ($BD{F}_N$) and space discretization is performed using the ${\mathcal{P}}_h^k$ discrete spaces introduced above. However, as is well-known for central differencing schemes in hyperbolic partial differential equations, the naive Galerkin discretization of Equation (11) on ${\mathcal{P}}_h^k$ can lead to spurious oscillatory instabilities. To circumvent this well-known problem, we stabilize the discrete advection equation. Four different stabilization methods are supported in our framework: the Streamline Upwind Petrov-Galerkin (SUPG, [28]), the Galerkin Least Squares (GLS, [29,30]), the Sub-Grid Scale (SGS, [31]) and the Continuous Interior Penalty (CIP, [32]) methods. Detailed description of these methods can be found in [33,34,35]. It should be noted that the CIP stabilization is much more costly than the three others, as it densifies the corresponding algebraic system, and requires larger stencils and thus larger connectivity tables.

In the benchmarks, we used the GLS stabilization method, introducing the bilinear form

$$S_{GLS}=\sum _{K\in {\mathcal{T}}_h}{\int }_K\tau \mathcal{L}[\psi ]\left(\mathcal{L}[\varphi ]-\sum _{i=1}^N{\alpha }_i{\varphi }^{(n+1-i)}\right)$$

(16)

with $\mathcal{L}[\varphi ]=\frac{{\alpha }_0}{\delta t}\varphi +\overrightarrow{u}·\nabla \varphi$, $\tau$ a coefficient chosen to adjust the stabilization to the local advection strength. The choice for the parameter $\tau$ was extensively discussed, in particular for the case of advection-diffusion-reaction (c.f. [31,36,37]), and we provide several of them in our AdvectionDiffusionReaction toolbox. However, for the specific case of pure transient advection, which can also be seen as an advection-reaction equation after time discretization, we use the simple expression

$$\tau =\frac{1}{\frac{2|\overrightarrow{u}|}{h}+\frac{2{\alpha }_0}{\delta t}}$$

(17)

In summary, the discrete FEM we solve for the level-set at time $t+\delta t$ given the values at previous times is

$$\begin{array}{c}\mathit{Find}{\varphi }^{(n+1)}\in {\mathcal{P}}_h^{k}s.t.\forall \psi \in {\mathcal{P}}_h^k,\hfill \\ {\int }_{{\mathsf{\Omega }}_h}\left[\frac{{\alpha }_0}{\delta t}{\varphi }^{(n+1)}+{\overrightarrow{u}}^{(n+1)}·\nabla {\varphi }^{(n+1)}\right]\psi +S\left({\varphi }^{(n+1)},\psi ;{\overrightarrow{u}}^{(n+1)},{\{{\varphi }^{(i)}\}}_{i\le n}\right)={\int }_{{\mathsf{\Omega }}_h}\sum _{i=1}^N\frac{{\alpha }_i}{\delta t}{\varphi }^{(n+1-i)}\psi \hfill \end{array}$$

(18)

where $S\left({\varphi }^{(n+1)},\psi ;{\overrightarrow{u}}^{(n+1)},{\{{\varphi }^{(i)}\}}_{i\le n}\right)$ is the SUPG, GLS, SGS or CIP stabilization bilinear, which vanishes as $h\to 0$ and in the three first cases involves the previous time steps to ensure consistency of the stabilized equation with its continuous version.

The assembly is performed in parallel with automatic choice of the appropriate quadrature order and seamless inter-process communication management.

The linear system is then solved using the PETSc library [38] with a GMRES solver preconditioned with an additive Schwartz Method (GASM [39]).

3.5.2. Geometrical Quantities

In many physical situations, the immersed structures play an active role in the dynamics of the whole system through the surface forces exerted by their “membranes” on the fluid. These forces can most of the time be computed from the geometry of the surface as they can be related to deformations of the interface. The level-set description of interfaces is very convenient when it comes to the computation of the geometrical parameters required in the interface forces. As an example, the normal and curvature of the surface can naturally be obtained from $\varphi$ as

$$\begin{array}{c}\hfill \overrightarrow{n}=\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\end{array}$$

(19)

$$\begin{array}{c}\hfill \kappa =\nabla ·\overrightarrow{n}\end{array}$$

(20)

These quantities are defined on the whole domain in an Eulerian way, and can therefore naturally be used in surface integrals as introduced in Equation (9).

To compute such derivatives of fields, we use projection operators in order to work with fields still in ${\mathcal{P}}_h^k$. To this end, we introduce the $L^2$ projection operator ${\mathsf{\Pi }}_{L^2}^k$ and the smoothing projection operator ${\mathsf{\Pi }}_{sm}^k$, defined respectively as

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{L^2}^k:& V_h\to {\mathcal{P}}_h^k\hfill \\ & u\mapsto u_{L^2}\equiv {\mathsf{\Pi }}_{L^2}^k(u)\mathrm{s}.\mathrm{t}.\forall v\in {\mathcal{P}}_h^k\hfill \end{array}$$

$${\int }_{{\mathsf{\Omega }}_h}{u}_{L^2}v={\int }_{{\mathsf{\Omega }}_h}uv$$

(21)

and

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{sm}^k:& V_h\to {\mathcal{P}}_h^k\hfill \\ & u\mapsto u_{sm}\equiv {\mathsf{\Pi }}_{sm}^k(u)\mathrm{s}.\mathrm{t}.\forall v\in {\mathcal{P}}_h^k\hfill \end{array}$$

$${\int }_{{\mathsf{\Omega }}_h}{u}_{sm}v+\epsilon {\int }_{{\mathsf{\Omega }}_h}\nabla u_{sm}:\nabla v-\epsilon {\int }_{\partial {\mathsf{\Omega }}_h}\left(\overrightarrow{N}·\nabla u_{sm}\right)v={\int }_{{\mathsf{\Omega }}_h}uv$$

(22)

with $\epsilon$ a—small—smoothing parameter typically chosen as $\epsilon \approx 0.03h/k$. These projection operators can be defined for both scalar and vector fields in the same way by using the appropriate contractions.

Note that the smoothing projection operator introduces some artificial diffusion which is controlled by $\epsilon$. This diffusion is reponsible for the smoothing of the projected field, but can also introduce artefacts in the computation, so that $\epsilon$ needs to be carefully chosen depending on the simulation. In practice, the smoothing operator ${\mathsf{\Pi }}_{sm}^k$ is used to compute derivatives of order $\ge k+1$ of fields in ${\mathcal{P}}_h^k$, as such order of derivation are subject to noise.

The projection operators are implemented by the Projector class, which optimizes the assembly of the algebraic systems corresponding to Equations (21) and (22) by storing the constant terms—such as the mass matrix for example—to prevent unnecessary computations at each projection. We usually also store the preconditioner of the system for reuse.

The projection linear systems are solved using PETSc’s GMRES solver with a GASM preconditioning method.

In the benchmarks below, we used $\overrightarrow{n}={\mathsf{\Pi }}_{L^2}^1\left(\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$ and $\kappa ={\mathsf{\Pi }}_{sm}^1(\nabla ·\overrightarrow{n})$ for ${\mathcal{P}}_h^1$ simulations, and $\overrightarrow{n}={\mathsf{\Pi }}_{L^2}^1\left(\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$ and $\kappa ={\mathsf{\Pi }}_{L^2}^1(\nabla ·\overrightarrow{n})$ for higher order ones.

3.5.3. Redistantiation

As mentioned above, as the level-set is advected by the fluid, it loses its “distance” property, i.e., $|\nabla \varphi |=1$ is not in general preserved by the advection equation. Therefore, to ensure numerical stability and prevent accumulation or rarefaction of the level-set iso-lines which support the regularized Dirac and Heaviside functions, one needs to redistantiate the level-set, i.e., recover a signed distance function to $\{\varphi =0\}$. Too main approaches can be used: the first one relies on a efficient direct computation of the distance using the well-known “fast-marching algorithm” [20], while the second consists in solving an Hamilton-Jacobi equation [18,19] which steady state enforces $|\nabla \varphi |=1$. Both methods are implemented in our framework; we provide some details in the following.

The fast-marching method is an efficient algorithm to solve the Eikonal equation in general, with an algorithmic complexity of $\mathcal{O}\left(N_{dof}log(N_{dof})\right)$. It uses the upwind nature of this equation to “march” away from the 0-iso-level to the rest of the domain, hence the Dijkstra-like complexity. Our implementation of the fast-marching algorithm is based on [35], and therefore strongly relies on $\varphi$ being a ${\mathcal{P}}_h^1$ field.

If the level-set is of order $k>1$, we can still perform the fast-marching redistantiation using a ${\mathcal{P}}_{iso}^1$ space as proposed in [40]. This order-1 Lagrange function space is constructed on a mesh obtained by replacing all the ${\mathcal{P}}_h^k$ degrees of freedom by nodes. The ${\mathcal{P}}_{iso}^1$ hence has the same number of degrees of freedom than the original ${\mathcal{P}}_h^k$ one. We can then project ${\varphi }^k\in {\mathcal{P}}_h^k$ onto ${\varphi }_{iso}^1\in {\mathcal{P}}_{iso}^1$, perform the fast-marching on ${\varphi }_{iso}^1$, and project back.

The Hamilton-Jacobi method follows a different approach, and works directly within the finite-element framework. It consists in solving an Hamilton-Jacobi-like equation which steady states is a signed distance function, namely

$${\partial }_{\tau }\varphi +sign(\varphi )\left(1-|\nabla \varphi |\right)=0$$

(23)

The $sign(\varphi )$ term anchors the position of $\{\varphi =0\}$, and ensures that the redistantiation front gets transported from this 0-iso-level, inward or outward depending on the initial sign of $\varphi$. In practice, this equation is discretized using the advection framework introduced above as

$$\begin{array}{c}\mathit{Find}{\varphi }^{(n+1)}\in {\mathcal{P}}_h^{k}s.t.\forall \psi \in {\mathcal{P}}_h^k\hfill \\ {\int }_{{\mathsf{\Omega }}_h}\left[\frac{{\alpha }_0}{\delta \tau }{\varphi }^{(n+1)}+sign({\varphi }^{(0)})\frac{\nabla {\varphi }^{(n)}}{|\nabla {\varphi }^{(n)}|}·\nabla {\varphi }^{(n+1)}\right]\psi +S\left({\varphi }^{(n+1)}\right)={\int }_{{\mathsf{\Omega }}_h}\left[sign({\varphi }^{(0)})+\sum _{i=1}^N\frac{{\alpha }_i}{\delta \tau }{\varphi }^{(n+1-i)}\right]\psi \hfill \end{array}$$

(24)

with $sign({\varphi }^{(0)})\equiv 2H_{ϵ}({\varphi }^{(0)})-\frac{1}{2}$ and $S\left({\varphi }^{(n+1)}\right)$ a stabilization bilinear form. This equation is usually solved for a few iterations with $\delta \tau \sim h$. The main advantage of the Hamilton-Jacobi approach is that it can be used straigthforwardly for high-order level-set functions. However, it is in general slower than the fast-marching method, and numerical errors can quickly accumulate and lead to spurious motion of the interface resulting in strong loss of mass of the $\{\varphi \le 0\}$ domain. To avoid these pitfalls, the number of iterations and the pseudo time-step $\delta \tau$ must be carefully chosen, which is not an easy task in general.

4. 3D Rising Drops Benchmark

We now present a 3D benchmark of our numerical approach, using the Navier-Stokes solver developed with the Feel++ library described in [41]. This benchmark is a three-dimensional extension of the 2D benchmark introduced in [42] and realised using Feel++ in [43]. The setup for this benchmark was also used in [44] to compare several flow solvers.

4.1. Benchmark Problem

The benchmark consists in simulating the rise of a 3D drop in a Newtonian fluid. The equations solved are the aforementioned coupled incompressible Navier Stokes equations for the fluid and advection equation for the level set Equations (9)–(11) with ${\overrightarrow{f}}_b$ and ${\overrightarrow{f}}_s$ respectively the gravitational and surface tension forces, defined as:

$$\begin{array}{cc}\hfill {\overrightarrow{f}}_b& ={\rho }_{\varphi }\overrightarrow{g}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\overrightarrow{f}}_s& =\sigma \kappa \overrightarrow{n}{|}_{\mathsf{\Gamma }}\simeq \sigma \kappa \overrightarrow{n}{\delta }_{\epsilon }(\varphi )\hfill \end{array}$$

(26)

with $\overrightarrow{g}\equiv -0.98{\overrightarrow{e}}_z$ the gravity acceleration and $\sigma$ the surface tension.

We consider $\mathsf{\Omega }$ a cylinder with radius $R=0.5$ and height $H=2$, filled with a fluid and containing a droplet of another imiscible fluid. We denote ${\mathsf{\Omega }}_1=\{\overrightarrow{x},s.t.\varphi (\overrightarrow{x})>0\}$ the domain outside the droplet, ${\mathsf{\Omega }}_2=\{\overrightarrow{x},s.t.\varphi (\overrightarrow{x})<0\}$ the domain inside the drop and $\mathsf{\Gamma }=\{\overrightarrow{x},s.t.\varphi (\overrightarrow{x})=0\}$ the interface. We impose no-slip boundary conditions $\overrightarrow{u}{|}_{\partial \mathsf{\Omega }}=0$ on $\mathsf{\Omega }$ walls. The simulation is run from $t=0$ to 3.

Initially, the drop is spherical with radius $r_0=0.25$ and is centered on the point $(0.5,0.5,0.5)$ assuming that the bottom disk of the $\mathsf{\Omega }$ cylinder is centered at the origin. Figure 1 shows this initial setup.

We denote with indices 1 and 2 the quantities relative to the fluid in respectively ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$. The parameters of the benchmark are then ${\rho }_1$, ${\rho }_2$, ${\mu }_1$, ${\mu }_2$ and $\sigma$. We also define two dimensionless numbers to characterize the flow: the Reynolds number which is the ratio between inertial and viscous terms and is defined as

$$\mathrm{Re}={\displaystyle \frac{{\rho }_1\sqrt{|\overrightarrow{g}|{(2r_0)}^3}}{{\nu }_1}},$$

and the Eötvös number represents the ratio between the gravity force and the surface tension

$${\mathrm{E}}_0={\displaystyle \frac{4{\rho }_1|\overrightarrow{g}|r_0^2}{\sigma }}.$$

Table 1. Numerical parameters taken for the benchmarks.

Tests	${\mathit{\rho }}_1$	${\mathit{\rho }}_2$	${\mathit{\nu }}_1$	${\mathit{\nu }}_2$	$\mathit{\sigma }$	Re	${\mathbf{E}}_0$
Case 1 (ellipsoidal drop)	1000	100	10	1	24.5	35	10
Case 2 (skirted drop)	1000	1	10	0.1	1.96	35	125

reports the values of the parameters used for two different test cases proposed in [44]. At $t=3$, the first one leads to an ellipsoidal-shaped drop while the second one gives a skirted shape due to the larger density and viscosity contrasts between the inner and outer fluids.

To quantify our simulation results, we use three quantities characterizing the shape of the drop at each time-step: the center-of-mass

$${\overrightarrow{x}}_c=\frac{1}{|{\mathsf{\Omega }}_2|}{\int }_{{\mathsf{\Omega }}_2}\overrightarrow{x},$$

the rising velocity—focusing on the vertical component

$${\overrightarrow{u}}_c=\frac{1}{|{\mathsf{\Omega }}_2|}{\int }_{{\mathsf{\Omega }}_2}\overrightarrow{u},$$

and the sphericity—defined as the ratio between the area of a sphere with same volume and the area of the drop—

$$\mathsf{\Psi }={\displaystyle \frac{4\pi {\left({\displaystyle \frac{3}{4\pi }}|{\mathsf{\Omega }}_2|\right)}^{\frac{2}{3}}}{|\mathsf{\Gamma }|}}.$$

Note that in the previous formulae, we have used the usual “mass” and area of the drop, respectively defined as $|{\mathsf{\Omega }}_2|={\int }_{{\mathsf{\Omega }}_2}1$ and $|\mathsf{\Gamma }|={\int }_{\mathsf{\Gamma }}1$.

4.2. Simulation Setup

The simulations have been performed on the supercomputer of the Grenoble CIMENT HPC center up to 192 processors. To control the convergence of our numerical schemes, the simulations have been run with several unstructured meshes, which characteristics are summarized in

Table 2. Mesh properties and degrees of freedom: mesh characteristic size, number of tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom and total number of degrees of freedom of the simulation.

Mesh Properties			Finite-Element DOF
$\mathit{h}$	Tetrahedra	Points	Order 1 DOF	Order 2 DOF	#DOF
0.025	380,125	62,546	62,546	490,300	1,595,992
0.02	842,865	136,932	136,932	1,092,644	3,551,796
0.0175	1,148,581	186,136	186,136	1,489,729	4,841,459
0.015	1,858,603	299,595	299,595	2,415,170	7,844,700
0.0125	2,983,291	479,167	479,167	3,881,639	12,603,251

.

We run the simulations looking for solutions in finite element spaces spanned by Lagrange polynomials of order $(2,1,k)$ for respectively the velocity, the pressure and the level set. The corresponding numbers of degrees of freedom for each mesh size are reported in Table 2 and the numerical parameters used for the simulations are provided in

Table 3. Numerical parameters used for simulations and resulting simulation times for each test case.

Numerical Parameters			Total Time (h)
$\mathit{h}$	#proc	$\mathbf{\Delta }\mathit{t}$	Case 1	Case 2
0.025	64	$1\times {10}^{-2}$	3.5	3.6
0.02	128	$9\times {10}^{-3}$	4.8	5.1
0.0175	128	$8\times {10}^{-3}$	8.9	9.5
0.015	192	$7\times {10}^{-3}$	12.3	13.5
0.0125	192	$6\times {10}^{-3}$	33.8	39.6

.

The Navier-Stokes equations are solved using Newton’s method and the resulting linear system is solved with a preconditioned flexible Krylov gmres method using the simple preconditioner introduced in [45]. The inner “inversions” of the velocity block matrix are performed using a block Jacobi preconditioner, with an algebraic multigrid (gamg) preconditioner for each velocity component block. Note that these velocity block inversions are only preconditioned—without any KSP iteration. The nested Schur complement “inversion” required for simple is solved with a few iterations of gmres preconditioned with an algebraic multigrid (gamg).

The linear advection equation is solved with a Krylov gmres method, preconditioned with an Additive Schwarz Method (gasm) using a direct $LU$ method as sub-preconditionner.

5. Results

In this section, we analyze the simulations of the rising drop for the two cases with low-order—$k=1$—Section 5.1 and Section 5.2 and high-order—$k=2$—Section 5.3 level set discretization space. Except for the comparison between the fast-marching and the Hamilton-Jacobi methods Section 5.1.1, the level set redistantiations were performed with the fast-marching method every 10 time-steps.

5.1. Case 1: The Ellipsoidal Drop

a shows the shape of the drop in the x-z plane at the final $t=3$ time step for the different aforementioned mesh sizes. The shapes are similar and seem to converge when the mesh size is decreasing. The drop reaches a stationary circularity as shown in Figure 2d, and its topology does not change. The velocity increases until it attains a constant value. Figure 2c shows the results obtained for the different mesh sizes. The evolution of the mass of the drop versus time is shown in Figure 2e. It highlights the rather good mass conservation property of our simulation setup, as about $3\%$ of the mass is at most lost for the coarsest mesh, while the finest one succeeds in keeping the loss in mass below $0.7\%$.

We also note that our simulation perfectly respects the symmetry of the problem and results in a axially symmetric final shape of the drop, as shown in Figure 3.

5.1.1. Comparison between Hamilton-Jacobi and Fast-Marching Reinitialization

As mentioned in Section 3.5.3, two reinitialization procedures can be used to overcome the “deformation” of the level set which becomes more and more different from the distance to the interface function as it is advected with the fluid velocity. The fast-marching method resets the values of $\varphi$ on the degrees of freedom away from the interface to match the corresponding distance. The Hamilton-Jacobi method consists in solving an advection equation which steady solution is the wanted distance function.

We have run the $h=0.0175$ simulation with both reinitialization methods to evaluate the properties of each one, and compare them using the monitored quantities. Figure 4 shows the respective shapes of the drop at final time ($t=3$), and gives the obtained results.

The first observation is that the mass loss (see e) is considerably reduced when using the FM method. It goes from about 18% mass lost between $t=0$ and $t=3$ for the Hamilton-Jacobi method to less than 2% for the fast-marching method. This resulting difference of size can be noticed in Figure 4. The other main difference is the sphericity of the drop. Figure 5d shows that when using the fast-marching method, the sphericity decreases really quickly and stabilises to a much lower value than the one obtained with the Hamilton-Jacobi method. This difference can be explained by the fact that the fast-marching method does not smooth the interface. The shape can then contain some small irregularities leading to a bad sphericity. Even so, with both methods the sphericity stays quite constant after the first second of the simulation. The rising velocity and the vertical position do not show any significant difference between the two reinitialization methods.

5.1.2. Comparison with Previous Results

Figure 6 shows a plot of our results compared to the ones presented in [44]. In this paper, the authors perform simulations on the same setup and with the same test cases as considered here. To ensure consistency of their results, they use three different flow solvers (hence three different space discretization methods) coupled with two different interface capturing methods: the DROPS and NaSt3D solvers coupled to a level set approach, and the OpenFOAM solver which uses a volume-of-fluid method.

To evaluate the effect of the characteristic mesh size, we plot the results we obtained for the simulations run with both $h=0.025$ and $h=0.0125$ along with the results from [44].

We can observe an overall good agreement between our results and the benchmark performed in [44].

5.2. Case 2: The Skirted Drop

In the second test case, the drop gets more deformed because of the lower surface tension and the higher viscosity and density contrasts. Figure 7 displays the monitored quantities for this test case. We observe that the shape of the “skirt” of the drop at $t=3$ is quite strongly mesh dependent, but converges as the mesh is refined. The other characteristics of the drop are not so dependent on the mesh refinement, even for the geometrically related ones, such as the drop mass, which shows a really small estimation error (only $2\%$ difference between the coarsest and finest meshes), and displays the really good conservation properties of our simulations. We again also note in Figure 8 the symmetry of the final shape of the drop, which highlights the really good symmetry conservation properties of our approach.

5.2.1. Comparison between Hamilton-Jacobi and Fast-Marching Reinitialization

As for the test case 1, we provide a comparison of the results for the test case 2 obtained using either the fast-marching or the Hamilton-Jacobi reinitialization method. The comparison is performed with an average mesh size ($h=0.0175$). The obtained shapes are shown in Figure 9 and the results are given in Figure 10. As before, they highlight noticeable differences between the two methods for geometrically related quantities such as mass loss, sphericity and final shape. We can even observe a non-negligible difference for the latter in the region of the “skirt”. This difference, mainly related to the diffusive properties of the Hamilton-Jacobi method, can also be observed on the 3D shapes in Figure 9. The good agreement of the results obtained using the fast-marching method tends to suggest that the Hamilton-Jacobi method is not accurate enough—or would require more careful and costly adjustment of its parameters—for this kind of three-dimensional simulation.

5.2.2. Comparison with Previous Results

As in Section 5.1.2, we compare our results to the benchmark [44], and show the relevant quantities in Figure 11.

We also observe a good agreement between our simulations and the ones from the benchmark. We however note that the final shape of the skirted drop is very sensitive to the mesh and none of the groups agree on the exact shape which can explain the differences that we see on the parameters in Figure 11 at time $t>2$.

5.3. High-Order Simulations

As already mentioned, our framework naturally allows the use of high-order Galerkin discretization spaces. As an illustration, we present here benchmark simulation results performed using finite element spaces spanned by Lagrange polynomials of order $(2,1,2)$ and $(2,1,3)$ for each test case. The mesh size considered here is $h=0.02$, and the results are shown in Figure 12 and Figure 13 for test cases 1 and 2 respectively. We expect the increase in order of the level-set field to improve the overall accuracy.

We can indeed observe that the final shapes of high-order simulations look smoother in both cases, as confirmed by the sphericity plots. The effect is highly noticeable on the “skirt” which appears for the second test-case, which looks even smoother than the one obtained with the finest ($h=0.0125$) $(2,1,1)$ simulation.

Let us also highlight that the small differences observed with the $(2,1,3)$ simulations, in particular for the final shapes, are most likely related to the absence of articifial diffusion error in the computation of geometrical quantities (especially for the curvature), which suggest more robust and realistic results for these simulations.

We can also notice that more “physically” controlled quantities, such as the position of the center-of-mass and the vertical velocity are less impacted by the polynomial order of the level-set component, which is not so surprising, as these quantities are mainly determined by the (level-set-dependent) fluid equations, which discretization orders where kept constant for this analysis.

6. Conclusions and Outlooks

We have presented in this paper a comprehensive numerical framework for the simulation of Multifluid flows. This framework is based on level-set methods solved by a (possibly high-order) finite element method. The explicit coupling between the level-set and the fluid has proven to be efficient and has allowed us to take advantage of reliable and efficient preconditioning strategies to solve the fluid equations.

The framework has been implemented within the Feel++ toolboxes and leverages the efficiency of the library to run on large numbers of processors in parallel. It also features user-friendly interfaces, and allows for easy model setup and parametrization using the json and cfg standard formats. The use of state-of-the-art metaprogramming techniques allows to seamlessly perform simulations in two or three dimension, and to increase the polynomial order of the finite elements. We highlight again that both the implementation and the benchmarks configuration files are open-source and available online [14,15].

The presented Multifluid framework has been validated using a three-dimensional two-fluid numerical benchmark and achieved results in agreement with the simulations performed with other methods.

We have also compared two different level-set reinitialization procedures (the fast-marching and the Hamilton-Jacobi methods) and observed significantlty different behaviors, in particular the former is much better at mass conservation than the latter.

High-order simulations were performed to highlight the increased smoothness of the computed interfaces. High-order discretizations of the level-set function also greatly helps for the computation of geometrical quantities, such as the curvature of the interface. It avoids the need for artificial diffusion in the computation of such derivatives, which can be prove essential for accurate accounts of surface effects. In particular the $(2,1,3)$ simulations which feature complete diffusion-free geometrical quantities suggest that increasing the order of the level-set discretization can be of great interest when seeking highly accurate results regarding shapes, or when physical forces involving high derivatives of the level-set field are present.

Further improving the accuracy of the level-set and related quantities using higher order and/or hybrid methods is still ongoing.

The framework presented and validated here provides the building blocks for the simulation of complex fluids in complex geometries. The use of a high-order level-set method to track interfaces allows in particular the simulation of a large number of immersed vesicles or cells with implicit account of the surface forces at constant cost with respect to the number of particles. Its versatility shall be used in a near future to better understand the flow of blood cells in realistic vascular systems, or the dynamics of swimming droplets interacting with external surfactants.