A novel diffuse-interface model and a fully-discrete maximum-principle-preserving energy-stable method for two-phase flow with surface tension and non-matching densities

Two well-established classes of the interface capturing models are the level-set and phase-field models. Level-set formulations satisfy the maximum principle for the density but are not energy-stable. On the other hand, phase-field models do satisfy the second law of thermodynamics but lack the maximum principle for the density. In this paper we derive a novel model for incompressible immiscible two-phase flow with non-matching densities and surface tension that is both energetically-stable and satisfies the maximum principle for the density. The model finds its place at the intersection of level-set and phase-field models. Its derivation is based on a diffusification of the incompressible two-phase Navier–Stokes equations with non-matching densities and surface tension and involves functional entropy variables. Additionally, we present an associated fully-discrete energy-stable method. Isogeometric analysis is used for the spatial discretization and the temporal-integration is performed with a new time-integration scheme that is a perturbation of the second-order midpoint scheme. The fully-discrete scheme is unconditionally energy-dissipative, pointwise divergence-free and satisfies the maximum principle for the density. Numerical examples in two and three dimensions verify the energetic-stability of the methodology.


Introduction
Free-surface flows appear in a large class of applications ranging from marine and offshore engineering, e.g. violent sloshing of LNG in tanks or wave impacts, to bubble dynamics. Classical numerical methods for such free-surface flow problems follow the moving free-surfaces with mesh-motion. These so-called interfacetracking methods become very complex in case of topological changes (e.g. break-up or coalescence). A good alternative is formed by the interface-capturing methods which introduce an extra field that naturally captures the topological changes. This strategy transforms the moving boundary problem into a system of PDEs on a fixed domain. This significantly simplifies the complexity of associated numerical solution strategies, although it introduces new challenges.
Interface capturing methods can roughly be divided into volume-of-fluid, level-set and phase-field methods, see Elgeti and Sauerland [1] for a more elaborate discussion. Volume-of-fluid methods are popular methods for both incompressible flow [2][3][4] and compressible fluid flow [5,6] and their main advantage is mass conservation. Although, the need for compression techniques or interface reconstruction algorithms may destroy this feature. In general volume-of-fluid methods do not provably satisfy the maximum principle for the density, exceptions occur when a CFL-like condition is fulfilled, see e.g. [7]. This is precisely the property of level-set methods [8][9][10][11][12][13], they satisfy the maximum principle for the density ab initio (at the cost of mass conservation). An additional feature of the level-set method is the ability to accurately compute the surface tension force which typically based on the continuum model of Brackbill et al. [14]. It is well-known, see e.g. [15,16], that surface tension effects are better represented when using the level-set approach as compared with the volume-of-fluid approach. We refer to Gross and Reusken [17] for the error analysis of the surface tension force in the level-set method. The level-set appears in many applications including structural optimization [18] and image analysis [19], and is particularly a suitable tool in case of high density jumps [20][21][22]. The phase-field models [23][24][25] are known for their rigorous thermodynamical structure. Since phase-field models typically contain higher-order derivatives, isogeometric analysis is a popular discretization method of these models [26][27][28][29]. Phase field models find applications in many problems in computational mechanics including boiling [30], fracture [31,32] and tumor growth [33][34][35]. The main issue is that numerical methods for phase-field models do not provably satisfy the maximum principle for the density [36].
The natural notion of stability for nonlinear (interface) problems is entropy stability. In nonlinear problems entropy norms replace the Sobolev norms employed in the stability analysis of linear problems. Indeed, one obtains physically relevant solutions when imposing the second law of thermodynamics. Note that thermodynamically stable resembles in the current isothermal case energetically stable as Clausius-Duhem inequality reduces to an energy-dissipative inequality.
In this paper we derive a novel thermodynamically-stable maximum-principle-preserving diffuse-interface model for incompressible immiscible flows with surface tension and non-matching densities and present a corresponding fully-discrete finite element method. The model complements existing phase-field fluid models, and is innovative in the sense that it connects to classical level-set methods. It thereby presents a new class of methods in between phase-field and level-set methods. The fully-discrete method, based on the Galerkin finite element method, inherits several properties of the underlying model. Namely, it is energy-stable, it is pointwise divergence-free, and it satisfies the maximum principle for the density.
The first goal of this work is to derive from the sharp-interface incompressible Navier-Stokes equations with surface tension the following thermodynamically consistent diffuse-interface model: Here u is the velocity, p is the pressure and φ is the phase variable. The density ρ and the dynamic viscosity µ depend on the diffuse-interface thickness parameter ε and are computed as where the phase densities ρ 1 and ρ 2 and viscosities µ 1 and µ 2 are assumed constant. ϱ ε = ρ ′ ε (φ) denotes the density change with respect to the phase-field φ. Moreover, H ε = H ε (φ) denotes a polynomial regularization of the Heaviside functional (for the precise definition we refer to Section 6.2), δ ε (φ) = H ′ ε (φ) is the regularized Dirac delta functional with derivative δ ′ ε (φ). The variable v is a new contribution that we will discuss in the next paragraph. Furthermore, the g = −gȷ where g is the gravitational acceleration and ȷ is the vertical unit vector and y is the vertical position, and τ ε = τ ε (u, φ) = 2µ ε (φ)∇ s u is the diffusion tensor with ∇ s u the symmetric velocity gradient. The constant surface tension coefficient is σ . Moreover, the mean curvature functional is κ(∇φ) = ∇ · ν(∇φ) where ν(∇φ) = ∇φ/|∇φ| is the interface normal. Here |∇φ| is standard 2-norm of ∇φ, ∥∇φ∥ 2 , which may be regularized to avoid singularities as ∥∇φ∥ ϵ,2 , defined as ∥∇φ∥ 2 ϵ,2 = ∇φ · ∇φ + ϵ 2 . In the regularized case the left-hand side of (1e) is augmented with −σ δ ′ (φ)ϵ 2 /∥∇φ∥ ϵ,2 .
The distinctive component of model (1) is the introduction of the variable v in the last member of (1b) with its definition in (1e). To the best knowledge of the authors, this has not appeared in level-set formulations or methods before. In contrast, it is well-known in the phase-field community. Namely, it closely resembles (up to kinetic and gravitational contributions) the chemical potential. This is a familiar quantity in thermodynamics representing the partial molar Gibbs free energy and can be traced back to the seminal work of Gibbs [37]. We derive the variable v via the concept of functional entropy variables, proposed by Liu and coworkers [27], which is basically the same as the introduction of the chemical potential known in phase-field literature. One explanation of this approach is that it circumvents the limitation caused by the functional spaces as the new variable v equals the unavailable weighting function required to show energy dissipation. This creates extra freedom and as a result the associated weak form is equipped with energetic stability for standard divergence-conforming function spaces. This concept is the natural alternative to entropy variables when the mathematical entropy associated with the system of equations is a functional (instead of a function) of the conservation variables. To clarify, the key feature of the model (1) is that it is energy-stable in a standard functional setting, as opposed to the standard level-set formulation.
To explain why our new model (1) can be placed in between level-set and phase-field models we note the following. Phase-field models are diffuse-interface models. One approach to derive phase-field models is to start from a moving boundary problem, referred to as the sharp-interface theory, and then apply diffusification to smear out the layer. The sharp-interface problem can often be recovered by letting the interface thickness tend to zero and using matched asymptotic expansions [38]. The second way to obtain a phase-field model is via thermodynamical arguments in which free energy functionals play a key role. Remark that these phase-field models do not need to be associated with a sharp-interface model. The derivation of our model (1) falls into the first category in the sense we start with a sharp-interface model and subsequently apply diffusification. We also note that our level-set convection equation often appears in the derivation of a phase-field model via diffusification (e.g. in the Stefan problem [39]). On the other hand, we make use of free-energy functionals and our model is closely linked to the coupled Navier-Stokes-Cahn-Hilliard (NSCH) model [40] which falls into the second category. Both model (1) and the NSCH model represent incompressible two-component fluids with surface tensions. However, the coupled NSCH model assumes a constant density of the mixture. Gurtin et al. [41] have placed the NSCH model in a continuum mechanism framework using thermodynamical arguments and have proven compatibility with the second law of thermodynamics. Moreover, they have extended the model to the non-matching densities case. This more general model is closely linked to our model with non-matching densities.
The derivation of the diffuse-interface involves three approximations, one related to the diffusification while the other two are related to the regularization of the 2-norm. As just mentioned, we will see that the diffusification of the standard level-set formulation of the incompressible Navier-Stokes equations (performed in Eq. (43)) does not provide a thermodynamically consistent model. This is probably not surprising for phase-field modelers, as it is well-known in the phase-field community that a general diffusification does not necessarily provide a thermodynamically consistent model [39]. To ensure energy dissipation, one needs to diffusify the energies, and then re-derive the corresponding PDE model (by computing the variational derivative of the energy). On the other hand, it is presumably new to level-set experts. Therefore we discuss this, which is the starting point of the novel diffuse-interface model, in great detail. The regularization of the 2-norm singularity is the second approximation and is presented in (39)- (42). It appears that an additional modification, the third approximation given in (92)-(93), is required. Also this modification may not shock phase-field modelers as it stems from the surface energy, and thus the proposed model may more quickly be obtained via directly taking this road. To clarify, in absence of the 2-norm regularization diffusification is the only approximation.
The fact that the model (1) is energy-stable in the standard functional setting does not imply that the spatial and temporal discretization procedures are trivial. Both the spatial and temporal discretization procedures impose challenges. These challenges emerge from the complexity of the system, namely the nonlinearity and the strong coupling of the equations, and the need for convection stabilization techniques on the spatial level. Concerning the spatial discretization, we note that a standard Galerkin method for the phase-field/level-set equation (1d) is not stable procedure. However, employing a standard stabilization mechanism is not possible as it upsets the energetic stability of the system. As a second remark on stabilization techniques, we note that the construction of an energy-dissipative multiscale stabilization method, allowing the computation of high Reynolds number flows, is not trivial and is beyond the scope of the current paper. Furthermore, the discrete approximation of the curvature is a challenge. A typical approach for lower-order methods is the use of a projection step which leads to inaccuracies. On that regard we note that in a recent paper [42] the authors propose to use a smooth higher-order NURBS-based isogeometric discretization [43] which significantly improves the accuracy of the curvature evaluation. Concerning the temporal discretization, the main challenge lies in the identification of the different energy evolution terms. This difficulty stems from the diffusification process.
We emphasize that level-set formulations are never equipped with a provably-stable algorithm. In fact, even on the semi-discrete level (both spatial discrete-time continuous and vice versa) stable numerical schemes for level-set formulations do not exist. It is considered one of the main discrepancies of the existing level-set methods and is of practical importance. Akkerman et al. [44] show that for a viscous air-water level-set simulation artificial energy may be created. This leads to a nonphysical prediction of the fluid behavior and renders the results useless.
The second objective is to propose an unconditionally energy-stable fully-discrete scheme for the proposed model (1). Additionally, the method satisfies the maximum principle for the density and is pointwise divergence-free. The formulation does not require the evaluation of the curvature and as such is not limited to higher-order discretizations. To the best knowledge of the authors this is the first energy-stable discretization scheme of a level-set type formulation. To inherit energetic stability in a spatial semi-discrete sense we employ a NURBS-based isogeometric analysis Galerkin-type discretization. Furthermore, we introduce an SUPG stabilization mechanisms that does not upset the energy-dissipative property of the method. Additionally, we augment the momentum equation with a residual-based discontinuity capturing term. An essential ingredient of the time-integration scheme, which can be understood as a perturbation of the second-order midpoint rule, is the dual-density-derivative approach in which different discretization strategies for the density derivative ϱ in (1b) and (1e) are utilized. Although the discretizations are very different, we note that a dual-density strategy has appeared in the design of a linear energy-stable method for a quasi-incompressible phase-field model [36].
The remainder of this paper is organized as follows. Section 2 presents and analyzes the energy behavior of the sharp-interface incompressible Navier-Stokes equations with surface tension. In Section 3 we use the sharp-interface model as a starting point to derive the diffuse-interface level-set model. Additionally, we provide a detailed analysis of the energy behavior. Phase-field modelers could jump over this analysis whereas level-set experts are encouraged to check it. We emphasize that the energy analysis does not extend to the standard functional setting. This is the motivation to derive our novel energy-dissipative formulation in Section 4. We do this with the aid of functional entropy variables and again note that it coincides with the introduction of the chemical potential (up to kinetic and gravitional contributions). Then, in Section 5 we present the semi-discrete finite element energetically-stable formulation. Next, in Section 6 we present our novel time-integration scheme to arrive at the fully-discrete energy-dissipative method. Section 7 shows the numerical experiments in two and three dimensions which verify the energy-dissipative property of the scheme. We draw conclusions in Section 8. Remark 1.1. To keep the work comprehensible we have intentionally not included multiscale stabilization mechanisms and the usage of techniques that control the interface-width. Incorporating these additional techniques in the currently proposed algorithm would allow the simulation of violent flows. These developments lie beyond the scope of this paper. △

Governing equations
Let Ω ⊂ R d , d = 2, 3, denote the spatial domain with boundary ∂Ω . We consider two immiscible incompressible fluids that occupy subdomains Ω i ⊂ Ω , i = 1, 2, in the senseΩ =Ω 1 ∪Ω 2 and Ω 1 ∩ Ω 2 = ∅. A time-dependent smooth interface Γ = ∂Ω 1 ∩ ∂Ω 2 separates the fluids. The problem under consideration consists of solving the incompressible Navier-Stokes equations with surface tension dictating the two-fluid flow: in Ω i (t) (3a) with u(x, 0) = u 0 (x) in Ω i (0) and Γ (0) = Γ 0 for the fluid velocity u : Ω → R d and the pressure p : Ω → R. The stress tensor is given by: with viscous stress tensor: The jump of a vector v is denoted as The problem is augmented with appropriate boundary conditions. We denote with x ∈ Ω the spatial parameter and with t ∈ T = (0, T ) the time with end time T > 0. Furthermore, we set g = −gȷ where g is the gravitational acceleration and ȷ is the vertical unit vector. The initial velocity is u 0 : Ω → R d . We use the standard convention for the various differential operators, i.e. the temporal derivative reads ∂ t and the symmetric gradient denotes . The constants µ i > 0 and ρ i > 0 denote the dynamic viscosity and density of fluid i respectively. The normal speed of Γ (t) is denoted as V , the normal of Γ (t), denoted ν, is pointing from Ω 2 (t) into Ω 1 (t) and the tangential vector is t. The curvature is κ = ∇ · ν, i.e. κ(x, t) is negative when Ω 1 (t) is convex in a neighborhood of x ∈ Γ (t). Furthermore, the outward-pointing normal of ∂Ω denotes n. We defined u n = u · n and u ν = u · ν as the normal velocity of ∂Ω and Γ (t), respectively. Eq. (3a) represents the balance of momentum while (3b) is the continuity equation. Next, (3c) states that the velocities are continuous across the separating interface. The fourth equation, (3d), stipulates that the discontinuity of the stresses at the interface is governed by surface tension. In absence of surface tension it reduces to an equilibrium of the stresses. Note that a direct consequence of (3d) is the continuity of tangential stress at the interface: We assume that the surface tension coefficient σ ≥ 0 is constant, i.e. Maragoni effects are precluded. Furthermore, we assume that line force terms vanish as a result of boundary conditions or additional conditions (see also [45]). We refer to [46] for some well-posed properties of the problem. We introduce the notation with indicator χ D of domain D. System (3) may now be written as: where τ (u) ≡ 2µ∇ s u and u(x, 0) = u 0 (x) in Ω i (0) and Γ (0) = Γ 0 . As we aim to develop an energy-stable model, we first study the energy behavior of the sharp-interface model associated with system (9). This is the purpose of the remainder of Section 2. After the energy analysis in Section 2.2 we present a standard weak formulation of (9) in Section 2.3.

Energy evolution
We consider the dissipation of the energy of the problem (9). The total energy consists of three contributions, namely kinetic (K ), gravitational (G) and surface energy (S): with y = x · ȷ the vertical coordinate. The subscript s refers to the sharp-interface model.
Theorem 2.1. Let u and p be smooth solutions of the incompressible Navier-Stokes equations with surface tension (9) The total energy E s , given in (10), satisfies the dissipation inequality: where B s contains the boundary contributions: Proof. To establish the dissipative property (11) we will first consider the evolution of each of the energy contributions (10) separately and subsequently substitute these in the strong form (9). We start off with the kinetic energy evolution. The following sequence of identities holds: where n 1 and n 2 denote the outward unit normal of Ω 1 (t) and Ω 2 (t), respectively. The first identity results from the Leibniz-Reynolds transport theorem. To obtain the second equality one adds a suitable partition of zero, subsequently applies the divergence theorem on both Ω 1 (t) and Ω 1 (t), and lastly uses the chain rule.
In a similar fashion we have the identities for the gravitational energy evolution: The first identity emanates from the Leibniz-Reynolds transport theorem and the second is a direct consequence of the divergence theorem. Finally, we consider the energetic contribution due to surface tension. We have from the Reynolds transport theorem in tangential calculus, see e.g. [47], the identity: where we recall that we do not account for Maragoni forces (σ is constant). Here ν ∂ is the unit-normal vector to ∂Γ (t), tangent to Γ (t). We refer to [48,49] for alternative insightful derivations of (15). We discard the last member of the right-hand side of (15) as it represents a line force.
We multiply the momentum equation by u and subsequently integrate over the domain: Considering the second expression in (16) in isolation we have the two identities: The first identity follows from adding a suitable partition of zero. For the second equality we perform integration by parts and make use of the jump (9d) where we note that on Γ (t) we have n 1 = −ν and n 2 = ν. We deduce from the continuity equation: Next, we collect the identities (13), (14), (15), (17) and (18), substitute these into (16). The first member in (16) cancels with the first term in the ultimate expression in (13). By virtue of (17) the second term in (16) drops out. The third member of (16) disappears due to (14). Some of the other terms in (13), (14) and (17) vanish due to (15) and (18). Gathering the expressions we eventually arrive at: This completes the proof with

Weak formulation
Recall that we suppress line force contributions as a result of boundary or auxiliary conditions. At this point we also assume homogeneous boundary conditions to increase readability of the remainder of the paper. Results can be easily extended to non-homogeneous boundary conditions. We define (·, ·) Ω as the L 2 (Ω ) inner product on the interior and (·, ·) Γ as the L 2 (Γ ) inner product on the boundary. We take zero-average pressures for all t ∈ T . The space-time velocity-pressure function-space satisfying homogeneous boundary condition u = 0 denotes S T and the corresponding weighting function space denotes S. The standard conservative weak formulation corresponding to strong form (9) reads: Find {u, p} ∈ S T such that for all {w, q} ∈ S: with interface speed V = u · ν. The weak formulation (21) is equivalent to the strong form (9) for smooth solutions and the associated energy evolution relation coincides with that of the strong form (9).

Remark 2.2.
To show the energy evolution for the case of non-homogeneous boundary conditions one may enforce boundary conditions with a Lagrange multiplier construction [50][51][52][53] and subsequently use (15) to identify the surface energy contribution. △ Remark 2.3. In order to avoid evaluating second-derivatives the alternative form +(∇w, σ P T ) Γ for the surface tension term in (21) with tangential projection P T = I − ν ⊗ ν may be used. In Appendix A.1 we provide the derivation of this alternative form. △

Standard diffuse-interface model
In this section we present the standard diffuse-interface model and analyze its energy behavior. To do so, in Section 3.1 we provide the sharp-interface level-set formulation of (9). Next, in Section 3.2 we present the corresponding non-dimensional form. Then in Section 3.3 we introduce two approximations: the diffusification and regularization to avoid singularities. We conclude with a detailed study of the energy behavior of this formulation in Section 3.4.

Sharp-interface level-set formulation
We employ the interface capturing level-set formulation to reformulate model (21). To this purpose we introduce the level-set function φ : Ω (t) → R to describe the evolution of the interface Γ (t). The sub-domains and interface are identified as: The motion of the interface Γ (t) is governed by pure convection: This results from taking the temporal derivative of the zero-level set. The domain indicator may now be written as: where H is the Heaviside function with the half-maximum convention: The resulting density and fluid viscosity are: In order to write the surface term in (21) in the level-set context we need expressions for the surface normal, the curvature and require to convert the surface integral into a domain integral. This is how we proceed. The surface normal is now continuously extended into the domain viâ The curvature results from taking the divergence of (27): .
By using a (double) variable transformation, we refer to Chang et al. [54], the surface integral of (21) converts into Here δ Γ = δ Γ (φ, ∇φ) denotes the Dirac delta concentrated on the interface Γ (t): which extends the integral over boundary Γ (t) to the domain Ω [55]. In (30) δ(φ) represents the standard Dirac delta distribution. Since (29) holds for all weights w we obtain the strong form equivalent to (3), (9) and weak form (21) in terms of the variables u, p and φ as: with u(x, 0) = u 0 (x) and φ(x, 0) = φ 0 (x) in Ω . Even though (31c) can be understood as both the level-set and phase-field equation, we refer to it as level-set convection equation throughout this paper. From this point onward we skip the hat symbols for simplicity.

Non-dimensionalization
We now perform the non-dimensionalization of the incompressible Navier-Stokes equations with surface tension. Here we re-scale the system (31) based on physical variables. The dimensionless variables are given by: where L 0 is a characteristic length scale and U 0 is a characteristic velocity. A direct consequence is where we have used the scaling property of the Dirac delta: The dimensionless system reads: The used dimensionless coefficients are the Reynolds number (Re) which expresses relative strength of inertial forces and viscous forces, the Weber number (We) measuring the ratio of inertia to surface tension and the Froude number (Fr) which quantifies inertia with respect to gravity. The expressions are given by: We suppress the star symbols in the remainder of this paper.

Diffusification
In the following we introduce the diffusification which is the first approximation of model. We diffusify the interface over an interface-width ε > 0 via replacing the (sharp) Heaviside function (25) by a diffusified differentiable Heaviside H ε (φ). We postpone the specific form of H ε (φ) to Section 6. The diffusified delta function is δ ε,Γ (φ, ∇φ) = δ ε (φ)∥∇φ∥ 2 with one-dimensional continuous diffusified delta function δ ε (φ) = H ′ ε (φ). We refer to [56] for details concerning the approximation of the Dirac delta. The density and fluid viscosity are computed as Our procedure to arrive at an energy-dissipative formulation, presented in Section 4, requires a conservative formulation of the momentum equation. Using the continuity and level-set convection equations, the diffuse-interface model follows straightforwardly: in Ω . At this point we have assumed a constant interface width ε.
Additionally we regularize the interface Dirac delta: The diffuse-interface model with regularization of the 2-norm reads: In the following we suppress the diffuse-interface width ε for the sake of notational simplicity.

Energy evolution
In the following we show the energy balance of the diffuse-interface formulation (43). First we introduce some notation, subsequently we present a few lemmas, and finally use these lemmas to establish the local and global energy balance of the model (43). We emphasize that the steps taken in this section are not valid in the standard functional setting. This is the motivation for the construction of the new diffuse-interface model in Section 4.
The kinetic, gravitational and surface energy associated with system (43) are: The total energy is the superposition of the separate energies: The local energy is given by: We present the local energy balance and subsequently the global balance. To that purpose we first need to introduce some notation and Lemmas associated with the surface energy. Let us define the normal projection operator: and the tangential projection operator: The associated gradient operators are the gradient along the direction normal to the interface: and the gradient tangent to the interface: Lemma 3.2. The term ∥∇φ∥ ϵ,2 evolves in time according to: Proof. This follows when evaluating the normal derivative of the level-set convection equation. Taking the gradient of the level-set convection equation and subsequently evaluating the inner product of the result with ν(∇φ) yields: Applying the gradient operator to each of the members provides + ∇u : The first term in (53) coincides with the first member in expression (51). For the second term in (53) we note that the term in brackets equals the gradient of ∥∇φ∥ ϵ,2 . Finally, one recognizes the normal projection operator in the latter term of (53). This delivers: Adding a suitable partition of zero completes the proof. □ Remark 3.3. The evolution (51) may be linked to the recently proposed variation entropy theory [57]. Variation entropy is local continuous generalization of the celebrated TVD (total variation diminishing) property derived from entropy principles. It serves as a derivation of discontinuity capturing mechanisms [58]. Using the continuity equation (43b) we obtain an alternative form of (51): with η(∇φ) = ∥∇φ∥ ϵ,2 and f(φ, u) = uφ. In the stationary case, i.e. when the term ∇ N (∂f/∂φ) is absent, relation (55) represents the evolution of variation entropy η(∇φ). This occurs when the velocity normal to the interface is constant. △ Lemma 3.4. The surface Dirac δ Γ ,ϵ (φ, ∇φ) evolves in time according to: Proof. Multiplying the level-set convection equation by δ ′ (φ) provides: The superposition of (51) multiplied by δ(φ) and (57) multiplied by ∥∇φ∥ ϵ,2 provides the result. In other words, the operator in which I denotes the identity operator, applied to the level-set convection equation delivers the evolution of the surface Dirac (56). □ To derive the local energy balance we introduce the following identity.
Proposition 3.5. It holds: Proof. See Appendix A.2. □ We now present the local energy balance.
Lemma 3.6. The local energy balance associated with system (43) takes the form: The divergence terms represent the redistribution of energy over the domain and the second to last term accounts for energy dissipation due to diffusion. The last term that emanates from the regularization is unwanted. We return to this issue in Section 4.
Proof. First we consider the local kinetic energy of the system (43). By straightforwardly applying the chain-rule we find: From the momentum and level-set convection equations, i.e. (43a) and (43c), we deduce: For the energetic contribution due the gravitational force, the chain-rule and the level-set equation (43c) convey that: And for the local surface energy evolution we invoke Lemma 3.4: Superposition of (62)-(64) yields: With the aim of simplifying (65) we introduce the identities: 1 We The first and the second identity follow from expanding the gradient and divergence operators. To obtain the third we note and apply Proposition 3.5 on the second term. Invoking (66) into (65) and adding a suitable partition of zero yields: With the aid of the continuity equation (43b) the latter member on the right-hand side of (68) vanishes. This completes the proof. □ Remark 3.7. The energy balance of Lemma 3.6 may also be written as: In this form we clearly see that the second divergence term represents the diffusion of kinetic energy density. △ We can now present the global energy evolution.
Theorem 3.8. Let u, p and φ be smooth solutions of the strong form (43). The associated total energy E , given in (45), satisfies the dissipation inequality: where B contains the boundary contributions: and where we have set ϵ = 0.
Proof. This follows from integrating the energy balance of Lemma 3.6 over Ω and using the divergence theorem: ∫ We discard the line force terms on the right-hand side and reorganize to get: Using the homogeneous boundary condition and setting ϵ = 0 finalizes the proof. □ The energy balance associated with the original model (9) and that of the level-set formulation (43) comply.
Proof. In the limit ε → 0 we may transform (73) back to get: To close this section we note that one may avoid evaluating second derivatives appearing in the surface tension term. This holds for the original model (9) which we have addressed with briefly in Remark 2.3. In the following Proposition we note that this alternative form directly converts to the regularized model (43). Proposition 3.10. We have the identity: We With the aid of Proposition 3.10 one can directly evaluate the surface tension term and does not require any additional procedure such as the one from [59].

The novel energy-dissipative diffuse-interface model
In this section we derive the novel diffuse-interface model. In Section 3 we have in great detail depicted the procedure to arrive at the energy dissipative statement for the standard diffuse-interface model. This procedure involves several steps that are not valid in the standard functional spaces. For instance the operator (58) associated with the surface energy is not permittable in a standard functional setting . Independently, the temporal discretization also gives rise to issues. Standard second-order semi time-discrete formulations of (43) are also not equipped with an energy-dissipative structure. We demonstrate this in Appendix B. Lastly, we note that the standard diffuse-interface model contains an unwanted term stemming from the regularization approximation.
The issues arise from the fact that the standard model is too restrictive with regard to the function spaces in order to establish energy-dissipation. Enlarging the standard function spaces introduces many complications and as such we do not further look into this strategy. The alternative is to modify the standard diffuse-interface model (43). This is the road we pursue. We employ the concept of functional entropy variables proposed by Liu et al. [27]. Liu and co-workers introduce the concept of functional entropy variables for the isothermal Navier-Stokes-Korteweg equations [27] and for the Navier-Stokes-Korteweg equations including the interstitial working flux term [30]. Here we apply the formalism to the level-set formulation of the incompressible Navier-Stokes equations with surface tension. This creates the extra space to resolve both discrepancies mentioned above. Additionally, the unwanted regularization term also vanishes. Remark that the usage of functional entropy variables closely resembles the introduction of the chemical potential (apart from kinetic and gravitational contributions) known by phase-field modelers.

Functional entropy variables
Energetic stability for the incompressible Navier-Stokes equations with surface tension coincides with stability with respect to a mathematical entropy function. Thus to construct an energy-dissipative formulation for the incompressible Navier-Stokes equations the natural approach seems to adopt entropy principles. For systems of conservation laws classical entropy variables are defined as the partial derivatives of an entropy with respect to the conservation variables. The Clausius-Duhem inequality plays the role of energetic stability and this results from pre-multiplication of the system of conservation laws by the entropy variables. The standard approach of constructing an entropy stable discretization as in Hughes et al. [60,61] is not applicable since the mathematical entropy is not an algebraic function of the conservation variables. In the situation of a general mathematical entropy functional the derivatives should be taken in the functional setting. The corresponding Clausius-Duhem inequality is then the result from the action of the entropy variables on the system of conservation laws.
In the current study we wish to inherit the notion of energetic stability for the incompressible model with surface tension. To this purpose we use as mathematical entropy functional the energy density (46) which we recall here: Following the approach described above, energetic stability results from the action of the entropy variables on the system of equations. In contrast to [27] and [30] the notion of conservation variables does not exist. Instead, the derivatives of H should here be taken with respect to the model variables U = (φ, ρu). Remark that (76) is a functional of the model variables U: Note that H contains a gradient term ∥∇φ∥ ϵ,2 which is non-local and thus the appropriate derivative is the functional derivative. We define the entropy variables as functional derivatives: The resulting functional derivatives are for test functions δv = [δv 1 , δv 2 , δv 3 , δv 4 ] T : We emphasize that it is essential to use the expression in terms of the model variables (77) to evaluate (79). The associated explicit form of (79) reads: We may use the functional entropy variables to systematically recover the energy balance (60).
Theorem 4.1. Applying the functional entropy variables to the incompressible two-phase Navier-Stokes equations with surface tension recovers the energy balance (60): Proof. Application of the functional entropy variables on the time-derivatives provides: Next we apply the entropy variables on the fluxes to get: Testing the entropy variables with the surface tension term gives: Testing the entropy variables with the viscous stress yields: And finally testing with the body force yields: Addition of (83)-(86) gives: Recognize the operator (58) on the fourth line of the right-hand side of (87). We may thus use Lemma 3.4 and write: 1 We Invoking the identities (66) and (88) the expression (87) collapses to a f a f dasd f a f se f segssegsgsgseg We merge the terms in (89) and use the continuity equation (43b) to cancel the terms containing the divergence of velocity. Taking the superposition of (82) and (89) while recognizing H on the right-hand side of (89) completes the proof. □

Modified formulation
Theorem 4.1 implies that an energy-dissipative relation may be recovered when the functional entropy variables are available as test functions. For standard test function spaces we cannot select the weight V 1 . We circumvent this issue, similar as in [27], by explicitly adding V 1 as a new unknown v to the system of equations. Thus we introduce the extra variable: where we use the notation ϱ = ϱ(φ) := ρ ′ (φ). The variable v equals the chemical potential (up to kinetic and gravitational contributions). Note that in a dimensional form the units of v do not match that of the chemical potential; the difference is a length unit stemming from the phase variable. The question arises how to couple the extra variable (90) to the diffuse-interface model (43). In this regard, note that a direct consequence of (90) is: Recall that the diffuse-interface model (43) is only associated with an energy-dissipative structure for ϵ = 0, see Theorem 3.8. This dissipative structure does not change when performing a consistent modification. Thus adding a suitable partition of zero based on (91) to the momentum equation (43a) keeps the same energy behavior. Instead, we suggest to replace the surface tension term in (43), i.e.
This is the third approximation in the model. In this way we eliminate the unwanted regularization term. We now obtain the final form of our new model in terms of the variables u, p, φ and v as: with u(0) = u 0 and φ(0) = φ 0 in Ω . This is model (1) in non-dimensional form.  (94) is equivalent to the standard diffuse-interface models (38) and (43).
The corresponding weak formulation reads: where we recall ϱ = ∂ρ/∂φ and have u(x, 0) = u 0 (x) and φ(x, 0) = φ 0 (x) in Ω . The solution space W T and corresponding test-function space W are divergence-compatible. We take W T := V T × Q 3 T and W := V × Q 3 where we refer to [62,63] for the precise definitions of V, V T and Q, Q T . 1. The formulation satisfies the maximum principle for the density, i.e. without loss of generality we assume that ρ 2 ≤ ρ 1 and then have: 2. The formulation is divergence-free as a distribution: 3. The formulation satisfies the dissipation inequality: Dissipation inequality (98) is not equipped with terms supported on the outer boundary ∂Ω since these vanish due to assumed boundary conditions. Proof. 1. This is a direct consequence of the definition of ρ = ρ(φ).
2. The divergence-conforming space allows to take q = ∇ · u in (95b) and hence we find: 3. Selection of the weights ψ = v in (95c) and ζ = −∂ t φ in (95d) yields: We add Eqs. (100) and find: Performing integration by parts yields: Recall that the line integral terms vanish due to auxiliary boundary conditions. Noting that we arrive at: Next we take w = u in (95a) to get: From the identities (66), the continuity equation (97), homogeneous boundary conditions and integration by parts we extract the identities: ) Ω = − (∇u, τ (u, φ)) Ω + (u, v∇φ) Ω .

Energy-dissipative spatial discretization
In this section we present the spatial discretization of the modified model (95). First we introduce some notation, then discuss the stabilization mechanisms and subsequently provide the semi-discrete formulation.

Notation
We employ an isogeometric analysis discretization. To provide the appropriate setting, we introduce the parametric domain denoted asΩ := (−1, 1) d ⊂ R d with corresponding mesh M. The element size h Q = diag(Q) of an element Q in M is its diagonal length. The physical domain Ω ⊂ R d follows as usual via the continuously differentiable geometrical map (with continuously differentiable inverse) F :Ω → Ω and the corresponding physical mesh reads: The Jacobian mapping is J = ∂ x/∂ξ . The physical mesh size h K is given by with the subscript F referring to the Frobenius norm. Note that on a Cartesian mesh it reduces to the diagonal-length of an element. The element metric tensor reads with inverse Using the metric tensor we see that the Frobenius norm is objective: where Tr denotes the trace operator. We define approximation spaces W h T ⊂ W T , W h ⊂ W spanned by finite element or NURBS basis functions. The div-conforming solution space is W h T := V h T × (Q h T ) 3 and the corresponding test-function space is W 0,h := V 0,h × (Q 0,h ) 3 . We refer to [63,64] for the precise definitions. Furthermore, we use the conventional notation superscript h to indicate the discretized (vector) field of the corresponding quantity.

Stabilization
It is well-known that a plain Galerkin discretization is prone to the development of numerical instabilities. This motivates the use of stabilization mechanisms. We employ the standard SUPG stabilization [65] for the level-set convection, i.e. we augment the discrete level-set equation with with residual We use the standard definition for stabilization parameter τ as also given in [52]. To ensure that the stabilization term does not upset the energetic stability property we balance it with the term: in the momentum equation.
Remark 5.1. In the current paper we focus on an energy-dissipative method without multiscale stabilization contributions in the momentum equation such as [66]. Standard stabilized methods are not directly associated with an energy dissipative property and thus specific techniques are required to establish such a property, see e.g. [53,67,68]. We note that these methods are developed for the single-fluid case. An extension to the current two-fluid case may be the topic of another paper. △ A popular method to stabilize the momentum equation is to use discontinuity capturing devices. We follow this road and augment the momentum equation with the discontinuity capturing term: The discontinuity capturing viscosity is given by: with conservative momentum residual and C a user-defined constant. The term clearly dissipates energy.
Remark 5.2. In order to avoid evaluating second derivatives in the surface tension contribution, one may project the residual onto the mesh and subsequently use Proposition 3.5. △ Remark 5.3. Even though we present the stabilization and discontinuity capturing terms in an ad hoc fashion, we wish to emphasize that these may be derived with the aid of the multiscale framework. The natural derivation for discontinuity capturing terms can be found in [58]. △

Semi-discrete formulation
The semi-discrete approximation of (95) is stated as follows: where u h (0) = u h 0 and φ h (0) = φ h 0 in Ω and we recall ϱ h (0) = ϱ(φ h (0)). The initial fields u h 0 and φ h 0 are obtained via standard L 2 -projections of respectively u 0 (x) and φ 0 (x) onto the mesh. The density and fluid viscosity are computed as The discrete counterparts of the kinetic, gravitational and surface energy are: The total energy is the superposition of the separate energies: Similarly, the semi-discrete local energy reads Remark 5.4. One may use a skew-symmetric form for the convective terms. Via a partial integration step, and using the continuity equation (38b), we may replace the convective term in (38) by the first three terms on the right-hand side of (124). In the current situation the specific form of the convective terms (conservative or skewsymmetric) is not essential. This changes when the formulation is equipped with multiscale stabilization terms.
In the single-fluid case (in absence of surface tension) the well-known multiscale discretization that represents an energy-stable system is the skew-symmetric form, see e.g. [52,53,68]. In contrast to the current two-phase model, this property is for the single-fluid case directly inherited by the fully-discrete case when employing the mid-point rule for time integration. △ Remark 5.5. We note that additional dissipation mechanisms for the surface evolution can upset energy-stability of the system. Well-balanced dissipation, introduced for the Navier-Stokes-Korteweg equations [69], is a possible strategy to resolve this. △ The semi-discrete formulation (119) inherits to a large extend Theorem 4.4. The notable difference lies in the usage of stabilization terms.
Theorem 5.6. Let (u h , p h , φ h , v h ) be a smooth solution of the weak form of incompressible Navier-Stokes equations with surface tension (119). The formulation (119) has the properties: 1. The formulation satisfies the maximum principle for the density, i.e. without loss of generality we assume that ρ 2 ≤ ρ 1 and then have: 2. The formulation is divergence-free as a distribution: 3. The formulation satisfies the dissipation inequality: The proof of Theorem 5.6 goes along the same lines as that of Theorem 4.4.
Proof. 1 & 2. The first two properties are directly inherited from the continuous case. Note that the weighting function choice for the second property is in general not permitted. The specific NURBS function spaces proposed by Evans et al. [63,64] do allow this selection.

Selection of the weights ψ
Addition of Eqs. (128) results in: (129) By performing integration by parts we obtain: Recognize δH h δφ h on the left-hand side to arrive at: Next we take w h = u h in (95a) to get: Similar as in the continuous case, we have the identities: Noting that δH h δ(ρ h u h ) = (u h ) T and employing (133) we arrive at: The superposition of (131) and (134) yields:

Energy-dissipative temporal discretization
In this section we present the energy-stable time-integration methodology. We present a modified version of the mid-point time-discretization method. First we introduce some required notation in Section 6.1 and then explain the time-discretization of the terms that differ from the standard midpoint rule in Sections 6.2 and 6.3. The eventual method is presented in Section 6.4.
The simplest fully-discrete algorithm would be to start from the semi-discrete version of (119) and then discretize in time using the second-order mid-point time-discretization. An important observation is that this approach does not lead to a provably energy-dissipative formulation, see Appendix B. We note that this is in contrast to the single-fluid case (in absence of surface tension effects).
In the following we present our strategy to arrive at a provable energy-dissipative formulation. Our approach is to mirror the semi-discrete case as closely as possible. We first focus on the terms that are directly associated with temporal derivatives of the energies and then discuss the discretization of the remaining terms.

Notation
Let us divide the time-interval T into sub-intervals T n = (t n , t n+1 ) (with n = 0, 1, . . . , N ) and denote the size of interval T n as time-step ∆t n = t n+1 − t n . We use subscripts to indicate the time-level of the unknown quantities, i.e. the unknowns at time-level n are u h n , p h n , φ h n and v h n . Lastly, we denote the intermediate time-levels and associated time derivatives as: where ρ h n = ρ(φ h n ), and defined [[a h ]] n := a h n+1 − a h n for the jump of a certain quantity a h .

Identification energy evolution terms
In order to identify the energy evolution terms we wish to have the fully discrete version of Three issues arise: (i) the approximation of the internal energy density 1 2 ∥u h ∥ 2 2 in the additional equation (119d), (ii) the approximation of the interface density jump term ϱ h and (iii) the approximation of the surface tension contribution. In the following we discuss the considerations for their time-discretization.
(i) The first matter is resolved when taking a shift in the time-levels in the energy density, analogously as in Liu et al. [27], i.e. we take 1 2 u h n · u h n+1 in the additional equation. (ii) Concerning the second problem, we require a stable time-discretization of ϱ h such that the approximation of ϱ h ∂ t φ h equals that of ∂ t ρ h . This suggests to approximate ϱ h at the intermediate time level t n+1/2 as such that The approximation (138) is not defined when φ h n+1 = φ h n . If ϱ h is a polynomial function of φ h we may use truncated Taylor expansions around φ h n+1/2 to find: where M chosen such that latter terms in the sum vanish and where we use the notation h (m) (x) = d m h/dx m for the mth derivative of a function h = h(x). Remark that (140) is well-defined. This motivates to use a (piece-wise) higher-order polynomial for ϱ h . We define the regularized Heaviside as where H p = H p (ϕ) is the piece-wise polynomial regularization: This function is C 3 -continuous at ϕ = 0 and C 3 -continuous at ϕ = −1, ϕ = 1. Furthermore, we base the regularization of Dirac on the Heaviside, i.e. we have δ ε (φ h ) = H (1) ε (φ h ).
Remark 6.1. The regularized Dirac delta δ ε (φ h ) has area 1. △ Remark 6.2. If ϱ h is non-polynomial one may use perturbed trapezoidal rules. In case of positive higher-order derivatives this leads to a stable approximation for ϱ h . △ Remark 6.3. This regularization closely resembles the popular goniometric regularization: with Taylor series representation: Definition (144) satisfies condition (139): Remark 6.4. The approximation ϱ h F,n+1/2 (case 2), defined in (138), is well-behaved when φ h n+1 ≈ φ h n . To see this, we consider without loss of generality the case where φ h n < −ε and −ε < φ h n+1 < 0. A Taylor series representation of H ε around φ h n+1 = −ε reveals: for some ξ ∈ (−ε, φ h n+1 ). It is now convenient to express φ h n as a perturbation of −ε relative to the distance between φ h n and −ε. In other words, we write φ h n = −ε − υ(φ h n+1 + ε) for υ = −(φ h n + ε)/(φ h n+1 + ε) > 0. Substitution into (147) gives: Noting that |H (4) P (ϕ)| < 60 we obtain the bound: (iii) We now turn our focus to the surface tension contribution, which writes in semi-discrete form: Recall that in the semi-discrete form the surface energy evolution follows when substituting ζ h = −∂ t φ h : Here we have utilized following identities: • for the first term: • for the second term: • and for combining the terms: We wish to follow the same steps in the fully-discrete setting. However, these identities are not directly guaranteed in a fully discrete sense. In the following we describe the fully-discrete approximation of each of the three terms in (150), i.e. δ ′ (φ h ), δ(φ h ) and ∇φ h /∥∇φ h ∥ ϵ,2 , that complies with these identities. To that purpose we introduce the mid-point approximation of the time-derivative. Proposition 6.5. The mid-point approximation of the time-derivative satisfies the product-rule in the following sense: where a h and b h are scalar or vector fields.
(III) We start off with the last identity (152c). The fully-discrete version of the product rule in (152c) follows from Proposition 6.5: This implies that we require the approximation: We now aim to identify the first and the second term on the right-hand side of (154) with first and second term on the right-hand side of (151) respectively. (I) To identify the first term we require, in a similar fashion as for ϱ, the approximation ς h n+1/2 ≈ δ ′ (φ h )(t n+1/2 ) to satisfy: To this purpose we define with truncated series: and the fraction: (II) We take the approximation: such that (152b) is satisfied in a fully-discrete sense:

Discretization other terms
We discretize the continuity equation using the mid-point rule, i.e.
( q h , ∇ · u h n+1/2 ) which implies pointwise divergence-free solutions on a fully-discrete level. Next, we require the fully-discrete version of the identities: which make use of the pointwise divergence-free property. These identities are fulfilled when we have Applying the chain-rule implies that we can take as approximation in the momentum equation: where the subscript m refers to the momentum equation.

Fully-discrete energy-dissipative method
We are now ready to present the fully-discrete energy-dissipative method: Given u h n , p h n , φ h n and v h n , find u h n+1 , p h n+1 , φ h n+1 and v h n+1 such that for all (w h , q h , ψ h , ζ h ) ∈ W 0,h : Remark 6.7. Due to Proposition 6.5 the time-derivative in the momentum equation may be implemented as: Theorem 6.8. The algorithm (166) has the properties: 1. The scheme satisfies the maximum principle for the density, i.e. without loss of generality we assume that ρ 2 ≤ ρ 1 and then have: 2. The scheme is divergence-free as a distribution: 3. The scheme satisfies the dissipation inequality: Proof. 1 & 2. Analogously to the semi-discrete case.

Selection of the weights
We add Eqs. (171) and find: We Using (146), (154), (156) and (161) we get Next we take w h = u h n+1/2 in (166a) to get: By virtue of (163) and (169) we have the identities: ) These reduce (174) to Addition of (173) and (176) by using (167) gives: Using the identity we identify the sum of the first two terms on the left-hand side of (177) as the change of kinetic energy. Next, the third term on the left-hand side of (177) represents change in gravitational energy. The latter term on the left-hand side of (177) resembles the surface energy evolution. We are left with: Remark 6.9. Following Brackbill [14] we employ the time-step restriction ∆t n ≤ ∆t max with 3 We 2π whereρ = (ρ 1 + ρ 2 )/2. △

Numerical experiments
In this section we evaluate the proposed numerical methodology on several numerical examples in two and three dimensions. To test the formulation we use both a static and dynamic equilibrium problem and check the energy dissipative property of the method. We do not test the method on a 'violent' problem in order to avoid the usage of redistancing procedures. All problems are evaluated with NURBS basis functions that are mostly C 1 -quadratic but every velocity space is enriched to cubic C 2 in the associated direction [63,64].

Static spherical droplet
Here we test the surface tension component of the formulation by considering a spherical droplet in equilibrium [14,42,70,71]. Viscous and gravitational forces are absent and hence the surface tension forces are in balance with the pressure difference between the two fluids. The interface balance (3d) thus reduces to: which is also referred to as the Young-Laplace equation. The exact curvature is given by: where we recall d = 2, 3 as the number of spatial dimensions. The spherical droplet of radius r = 2 of fluid 1 with density ρ 1 = 1.0 is immersed in fluid 2 with density ρ 2 = 0.1. The surface tension coefficient is σ = 73 which implies that surface tension forces dominate. The computational domain is a cubic with a side length of 8 units and the spherical droplet is positioned in the center of it. On all surfaces a non-penetration boundary condition (u n = 0) is imposed. We employ three meshes with uniform elements: 20 × 20, 40 × 40 and 80 × 80. We take ε = 2h K for all simulations in this section. The time-step is taken as ∆t n = 10 −3 which satisfies (180) for each of the meshes. We exclude the discontinuity capturing mechanisms for this problem, i.e. we set C = 0. In Fig. 2 we display the pressure for the finest mesh.
In Fig. 3 we display the pressure contours for each of the meshes. The corresponding pressure jump is 37.97, 36.80 and 36.56 for the meshes 20 × 20, 40 × 40 and 80 × 80 respectively. This implies second-order convergence.
In Figs. 4 and 5 we depict the energy evolution and dissipation for each of the meshes. The theoretical value of the surface energy is 2πr σ ≈ 917.34 which is well represented on the finest mesh. We see that the total and surface energies are (virtually) constant and the kinetic energy grows but has an insignificant contribution to the   total energy. Remark that this is not in conflict with the energy-dissipative property of the numerical discretization; even though kinetic energy increases, the total energy dissipates.
In Fig. 6 we report the velocity magnitude. We use the typical way to do so for this test case, i.e. report the velocity magnitude after one time-step and after 50 time-steps. Note that this test-case represents a physically-stable situation and as such it is desirable that velocities and thus the kinetic energy vanish. However, we see that this is not the case. Very small velocity currents are present. This is a well-known problem in the level-set community, and these currents are known as parasitic currents. The magnitude is comparable to what other researcher have obtained for this test-case. Several techniques can be used to reduce parasitic currents. One possibility is for example to use a so-called balanced-force algorithm [15] which assumes that the curvature is determined analytically. However, it is important to realize that the occurrence of parasitic currents is a property of the model, not of the method. The reason is that the system is not in a total energy-stable state. We explain this in more detail in Remark 7.1. As a consequence, numerical strategies that reduce parasitic currents are in some sense 'nonphysical tricks'. Moreover, we see in Figs. 4 and 5 that mesh refinement does not effect the kinetic energy evolution and dissipation. This confirms the claim that parasitic currents are not a discrepancy of the discretization method.
Remark 7.1. The occurrence of parasitic currents is a property of the diffuse-interface model (both of the standard model and the new model), not of the discretization. To see this, note that spurious oscillations would only be absent in a stationary state of the numerical model. A stationary state is the state with the minimum energy, and as such it is a minimizer of in which we recall that δ Γ ,ϵ (φ, ∇φ) is diffusified: δ ε,Γ ,ϵ (φ, ∇φ) = δ ε (φ)∥∇φ∥ ϵ,2 . The minimizers are characterized by  augmented by the conservation of mass constraint, for some constant C ∈ R. In absence of gravity we arrive at: Assuming that near the interface φ is not constant, we can discard the second term by setting ϵ = 0, which yields the requirement: Since the regularization of the Dirac delta functional is to some extend free, this equation only holds for C = 0, which results in two situations: • κ(∇φ) = 0: a straight interface.
In other words, parasitic currents occur near a curved interface. This explains why spurious currents, even though possibly very small, will always be present in this model (except for straight interfaces). Note that we recover the circular and spherical interface in 2 and 3 dimensions when taking the limit ε → 0, i.e. κ =C, on Γ (184) withC = −We C. △ In Fig. 7 we plot the variable v h n+1 . Note that the maximum theoretical value is where the max x∈Ω δ ε (φ) = max x∈Ω . We see that the finest mesh is able to accurately represent v h n+1 whereas on the coarser meshes v h n+1 is smeared out significantly.

Droplet coalescence 2D
In this example, inspired by Gomez et al. [72], we simulate the merging of two droplets into a single one. Gravitational forces are absent. Due to pressure and capillarity forces the single droplet then develops to a circular shape. We take as computational domain the unit box Ω = [0, 1] d and apply no-penetration boundary conditions. The initial configuration consists of two droplet at rest (u 0 = 0) with centers at c 1 = (0.4, 0.5) and c 2 = (0.78, 0.5) and radii r 1 = 0.25 and r 2 = 0.1 respectively. The regularized interfaces of the droplets initially overlap on a small part of the domain. If this is not the case the droplets remain at their position and thus no merging would occur. In contrast with the Navier-Stokes Korteweg equations, in this situation the interface has a finite width, due to the definition of H ε (φ). The Navier-Stokes Korteweg equations have no absolute notion of interface width; its effect decays exponentially. The droplets have a larger density (ρ 1 = 100) than the surrounding fluid (ρ 2 = 1) while the viscosities are equal: µ 1 = µ 2 = 1. We take as surface tension the low value of σ = 0.1 which causes a slowly merging process. To initialize the level-set we split the domain into two parts (x ≤ 0.665 and x > 0.665), such that each contains one droplet, and apply the standard distance initialization to each subdomain. We use 50 × 50 elements, set the time-step as ∆t = 0.1 and take C = 0.4.
We show in the Figs. 8 to 13 a detailed view of the merging process. The colors patterns are set per snapshot such that difference are most apparent.
In Figs. 14 and 15a we show the energy evolution and dissipation. In this case the theoretical value of the initial surface energy is 2π (r 1 + r 2 )σ ≈ 0.2199. We observe that the total and surface energies monotonically decrease in time. The kinetic energy increases when the droplets move towards each other (t < 10) and decreases during the merging process and subsequently flattens out.   In order to test whether the equilibrium state has been reached we evaluate the circularity of the droplet. The circularity is defined as the fraction of the perimeter evaluated from the droplet volume and the perimeter itself: The circularity depicted in Fig. 15b confirms the equilibrium state as γ tends to 1.

Droplet coalescence 3D
Here we simulate the merging of two droplets in three dimensions. We use the same physical parameters as in the two-dimensional case. The centers of the droplets are at c 1 = (0.4, 0.5, 0.6) and c 2 = (0.75, 0.5, 0.5) and the radii remain the same: r 1 = 0.25 and r 2 = 0.1. Also here the regularized interfaces of the droplets initially overlap. Again, to initialize the level-set we partition the domain, see Fig. 16a and apply the standard distance initialization to each subdomain. The initial configuration is depicted in Fig. 16b. We use 50 × 50 × 50 elements, set the time-step as ∆t = 0.1 and take C = 0.1.
We show in Fig. 17 snapshots of the merging process. In Fig. 18 we visualize the energy evolution and dissipation. The theoretical value of the initial surface energy is 4π (r 2 1 + r 2 2 )σ ≈ 0.0911. The behavior of the various energies  is similar as in the two-dimensional case. Also in this case the energy-dissipative property of the numerical method is confirmed.

Conclusion
In this work we have proposed a novel diffuse-interface model for immiscible two-phase flow with non-matching densities and surface tension. The model is derived from the sharp-interface incompressible Navier-Stokes equations with surface tension using the concepts of diffusification and functional entropy variables. The model lies at the intersection of level-set and phase-field models. The main properties are its energetic stability in the standard functional setting and the satisfaction of the maximum-principle for the density. In addition we have presented a corresponding fully-discrete energy-stable method that satisfies the maximum principle for the density and is pointwise divergence-free. We use isogeometric analysis for the spatial discretization and present a new perturbation of the second-order midpoint scheme for the time-integration. We have presented numerical examples in two and three dimensions which confirm the energy-stability of the method.  The energy-dissipative property makes the current model and associated numerical method a promising starting point for several extensions. A first suggestion is to equip the developed method with multiscale stabilization mechanisms that are energetically stable. The natural approach is to extend stabilization mechanisms that are energetically stable for single fluid flow [53,68]. Another possible research direction entails the control of the interface-width. The current model does not contain a physical mechanism that controls the interface-width. To gain control of the interface-width one could use re-distancing procedures that are well-known in level-set methods. To construct an energy-dissipative method that involves a re-distancing procedure one could attempt a two-step numerical scheme in which the redistancing procedure is decoupled from the current algorithm. These two extensions would allow a provably-stable simulation of a wide range of problems with non-matching densities. These include problems in which surface tension plays a key role, e.g. rising bubbles and falling droplets, and more violent problems in which inertia dominates, e.g. a dam-break problem.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
where Π Γ (x) is defined as the normal projector of x onto the interface Γ . The surface gradients of these fields are now given by while the tangential divergence of v is the trace of the surface gradient: Note the slight abuse of notation; we use the same notation for the surface gradient as employed for the surface gradient in the regularized level-set model. Alternative expressions for the surface gradients are where P T denotes the tangential projection tensor: Hereν is the continuous extension of the outward unit normal pointing from Ω 1 into Ω 2 and I is identity matrix.
Using the above identities we have Lemma A.1. Buscaglia et al. [49]: For any tangentially differentiable vector field w we have: Using (A.6) and Lemma A.1 we may write the surface tension term as

A.2. Diffuse-interface level-set model
In the following we utilize index notation.