An energy-preserving level set method for multiphase flows
AAn energy-preserving level set method for multiphase flows
N. Valle a , F.X. Trias a , J. Castro a a Heat and Mass Transfer Technological Centre (CTTC), Universitat Politècnica de Catalunya -BarcelonaTech (UPC), ESEIAAT, Carrer Colom 11, 08222 Terrassa (Barcelona)
Abstract
The computation of multiphase flows presents a subtle energetic equilibrium betweenpotential (i.e., surface) and kinetic energies. The use of traditional interface-capturingschemes provides no control over such a dynamic balance. In the spirit of the well-known symmetry-preserving and mimetic schemes, whose physics-compatible dis-cretizations rely upon preserving the underlying mathematical structures of the space,we identify the corresponding structure and propose a new discretization strategy forcurvature. The new scheme ensures conservation of mechanical energy (i.e., surfaceplus kinetic) up to temporal integration. Inviscid numerical simulations are performedto show the robustness of such a method.
Keywords:
Multiphase Flow, Symmetry-preserving, Mimetic, Conservative LevelSet, Energy-preserving
1. Introduction
Multiphase flows are ubiquitous in industrial applications. They are present in arich variety of physical phenomena such as vaporization [1], atomization [2], electro-hydrodynamics [3] or boiling [4], among others [5, 6].The use of interface-capturing schemes is widespread for the computation of mul-tiphase flows due to its computational e ffi ciency. The Volume-Of-Fluid (VOF) by Hirtand Nichols [7], the Level Set method developed by Osher and Sethian [8] and mostrecently phase field methods, introduced by Anderson et al. [9], are the most popular ∗ Corresponding author
Email addresses: [email protected] (N. Valle), [email protected] (F.X. Trias), [email protected] (J. Castro)
Preprint submitted to Elsevier September 4, 2019 a r X i v : . [ phy s i c s . c o m p - ph ] A ug nterface capturing schemes for multiphase flows. An overview of these can be foundin [10] and references therein. Despite the pros and cons that each method presents,we made our development concrete on the Conservative Level Set (CLS) initially de-veloped by Olsson and Kreiss [11] and Olsson et al. [12] due to its good conservationproperties, curvature accuracy and ease of handling topological changes. This wasextended to unstructured collocated meshes in [13].Of particular interest are the incompressible Navier-Stokes equations, ρ (cid:32) ∂(cid:126) u ∂ t + ( (cid:126) u · ∇ ) (cid:126) u (cid:33) = ∇ · σ ∇ · (cid:126) u = σ is composed of the hydrostatic and the deviatoric ones ( σ = − p I + τ ). In turn, τ is defined by Stokes constitutive equation τ = µ S , while the straintensor is given by S = / (cid:16) ∇ (cid:126) u + ( ∇ (cid:126) u ) T (cid:17) .The proper solution of equations (1) requires an appropriate decoupling of pressureand velocity. In this regard, the Fractional Step Method (FSM) [14] is an excellent toolwhich properly enforces the incompressibility constraint. However, the FSM results ina Poisson equation which needs to be solved, which takes most of the computationaltime in a typical simulation.The construction of discrete di ff erential operators in the seminal work of Verstap-pen and Veldman [15, 16] aims at preserving physical quantities of interest, namelymomentum and kinetic energy, by preserving several mathematical properties at thediscrete level. This merges with the conception of mimetic finite di ff erence methods[17], where the discretization is performed to satisfy the inherent mathematical struc-ture of the space, naturally producing a physics-compatible discretization. The presentwork is motivated by such an appealing idea. This has been named mimetic [17] ordiscrete vector calculus [18], among others [19, 20]. Mimetic methods delve into theconstruction of discrete di ff erential operators by producing discrete counterparts ofmore fundamental mathematical concepts, making extensive use of exterior calculus.This approach results in the algebraic concatenation of elementary operators, namelymatrices and vectors. Such an approach can be seen as the mathematical dual of thephysics-motivated work on symmetry-preserving schemes and provides with a di ff er-ent point of view which fortifies the analysis of this family of methods, which have2een used in both academic [21] and industrial problems [22, 23], among others. How-ever, to our knowledge, there is no a straightforward extension of these ideas into themultiphase flow community yet.While Direct Numerical Simulation (DNS) of single-phase flows has reached sub-stantial maturity, multiphase flows lag behind due to its increased complexity, namelydue to two main issues: i) surface tension and ii) di ff erences in physical properties. Theformer results in a dynamic equilibrium between kinetic and potential energies, whichare exchanged through the capillary term. Indeed this is the reason why symmetry-preserving discretizations [15, 16], despite conserving flawlessly kinetic (and thus to-tal) energy in single-phase flows, do not su ffi ce to preserve mechanical energy in mul-tiphase flows, as this transfer needs to be taken into account explicitly. The later poseschallenges regarding how interpolations need to be done without breaking physicallaws. In the framework of VOF, Fuster [24] proposed a discretization that preservesthe (skew-)symmetries of the momentum equation, preserving kinetic energy up tosurface tension, which is regarded as an energy source. However, as far as surface ten-sion is not included into the analysis, this is a necessary, but not a su ffi cient conditionfor preserving mechanical energy. Regarding the viscous term, the work of Sussman etal. [25] provided with a conservative discretization. The interested reader is referred to[26] and references therein for a comparison between di ff erent discretization strategiesfor the viscous term. Despite the impressive progress done so far, none of the abovehave included surface tension, and thus potential energy, in the analysis of conservationof energy. It is well-known, however, that the imbalances in the surface tension termmay lead to spurious currents and, eventually, unstable solutions [27]. In the frame-work of phase field methods, the impact of surface tension on the energy balance hasbeen included in the works of Jacqmin [28] for the Cahn-Hilliard equation, and Jametet al. [29] and Jamet and Misbah [30] for the Allen-Cahn formulation. This paper aimsto dig into a discretization including surface tension which, without recompression,preserves mechanical energy for level set schemes.The rest of the paper is arranged as follows: in Section 2 a glimpse of algebraictopology is provided. This sets the foundations to review the well-known symmetry-preserving discretization for single-phase flows in Section 3 and, inspired by this, de-3elop a new energy-preserving scheme for multiphase flows in Section 4. Comparativeresults between current techniques and the newly developed methods are presented inSection 5. Finally, conclusions and future insights are outlined in Section 6.
2. Topological model
Any numerical approach requires a finite-dimensional representation of the spacesunder consideration. This implies a discrete representation of the domains involvedin the setup of the problem. Single-phase flows fit well into a fixed frame, typicallydiscretized on a fixed grid. On the other hand, multiphase flows require to account fora moving interface which splits the domain at question into two regions. This inter-face needs to be properly discretized in order to preserve several inherent topologicalproperties. The way this is accomplished has led to a diversity of multiphase methods[10].
Let Ω be the working domain bounded by ∂ Ω and assume that M is a partition of Ω into a non-overlapping mesh. An illustrative example is given in Figure 1. Althoughwe stick to structured meshes for computational simplicity, the formulation presentedhere is independent of the mesh structure and thus can be extended into unstructuredmeshes. Incidence matrices are used to account for the connectivity within geometricentities. An example for T FC , the incidence matrix relating faces with cells accordingto the orientation of the mesh given in Figure 1 is shown next T FC = f f f f f f f f f f f f c − − + + c − − + + − c − − + + c − − + + (2)They replace usual neighboring relations such as φ c = (cid:80) f ∈ c φ f for the sum of facevalues related to cell c . In addition, its transpose provides with an explicit form for4 ∂ Ω c c c c ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n v v v v v v v v v Figure 1: Left: Domain Ω and its boundary ∂ Ω . Right: Mesh M . c i corresponds with the ith cell, ˆ n i corresponds with the normal vector to the jth face (i.e., f j ) and v k corresponds with the kth vertex. Itsincidence matrix is stated in equation (2). ∆ f φ = φ + − φ − for the di ff erence across face f , among others. Basic geometric proper-ties such as edge lengths ( W E ), face surfaces ( A F ) and cell volumes ( V C ) are arrangedas diagonal matrices. This matrix perspective presents several advantages: i) mesh in-dependence , ii) computational simplicity and iii) readily accessible algebraic analysis.While we restrain ourselves from digging into the first two, the later is useful both forreviewing the classical symmetry-preserving scheme and the development of the noveltechnique described here. Hereafter, lowercase letters correspond with vectors, whosesubscript indicates the geometric entity to which they are linked (e.g., p c corresponds topressure located at cells). Capital letters correspond with matrices, whose subscript(s)identify rows and (if di ff erent) columns (e.g. T FC is the face-to-cell incidence matrix). Interfaces imply a moving topology along the working domain, which implies a La-grangian frame of reference. Interface tracking schemes track such a frame explicitly,at the expenses of numerical complexity [31]. On the other hand, interface capturingschemes preserve a fully Eulerian approach, by mapping quantities expressed in theLagrangian frame back into the Eulerian one [7, 32, 33]. This results in a simpler im-plementation of the interface at the cost of an implicit representation. At this point5e split the presentation between the techniques used to actually capture the evolutionof the interface and the ones used to obtain explicit geometric information out of theimplicit form.
Let’s assume now that the domain Ω presents an interface at Γ , which splits Ω into Ω + and Ω − . We note that the volume of a single phase Ω + can be defined as (cid:90) Ω + dV = (cid:90) Ω H ( r ) dV (3)where r corresponds with the signed shorter distance from an arbitrary point to the in-terface, as can be seen in Figure 2, while H ( r ) is its corresponding Heaviside step func-tion, which is valued 1 at phase Ω + and 0 otherwise. Note that this function is the key tomap a Lagrangian frame ( Ω + ) back into an Eulerian one ( Ω ). Specific tracking of sucha quantity is the basis of the Volume of Fluid (VOF) method [7], which yields to theconcept of volume fraction or, more generally, marker function. Despite being formallyneat, the implementation of specific convection schemes is required, eventually requir-ing full geometric reconstruction, resulting in an intricate implementation. A di ff erentapproach is to capture the interface with a CLS [11, 12]. This captures the interface asthe isosurface of a continuous and smooth function θ . The level set marker function, θ ,results in a smoothed Heaviside step function that preserves (cid:82) Ω θ dV = (cid:82) Ω H ( r ) dV . It isconstructed as the convolution of the distance function r as follows H ( r ) ≈ θ ( r ) = (cid:18) tanh (cid:18) r (cid:15) (cid:19) + (cid:19) (4)where (cid:15) corresponds with a smoothing factor. Note that θ ( r ) → H ( r ) as (cid:15) →
0. Furtherdetails can be found in [13, 34].By imposing the conservation of the marker function, we can advect such a markerin an incompressible flow as [11, 12, 13]: ∂θ∂ t + ( (cid:126) u · ∇ ) θ = r c + r c − c + c − Figure 2: Distances r c ± are defined as the shorter distances of the interface to the cell. These are then normalto the interface and correspond with the minimum radius of the tangent sphere. Surface reconstruction may start by defining the interface normal. It is computedas [11] ˆ η i = ∇ θ |∇ θ | (6)Which implies, by definition, that the gradient of the marker function θ is parallelto the normal. On the other hand, curvature is computed as κ = −∇ · ˆ η i (7)Now, the surface area of Γ can be computed in any of the following forms A = (cid:90) Γ dA = (cid:90) Ω δ ( r ) dV = (cid:90) Ω ∇ H ( r ) · ˆ η i dV (8)where δ ( r ) is Dirac’s delta function, which formally is the distributional derivative ofthe Heaviside step function. This is the basis of the celebrated continuum surface forceof Brackbill et al. [35] for surface tension and, in general, of other smoothed interfacemethods. Regardless of the reconstruction method of choice, surface needs to satisfyfirst variation of area formula, which relates surface and volume variations throughcurvature and velocity as ddt (cid:90) Γ dA = − (cid:90) Γ κ(cid:126) u · ˆ η i dA (9)This is a fundamental identity, and the ultimate responsible of the correct conver-sion between kinetic and surface energy, as it will be shown in Section 4.1. A detailed7roof of this can be found in chapter 8.4 of Frankel [36]. As an illustrative example,let us consider the surface variation of a spherical surface. If we analyze how A = π r evolves under the action of the normal velocity, ˙ r , we obtain that dAdt = π r ˙ r , which canbe rearranged as dAdt = r A ˙ r , where we identify κ = r , the mean curvature of sphere.We are now going to prove that the use of a smooth marker function as in equation (4)leads to a consistent modeling of the interface by reconstructing surface area with itssmooth counterpart ˜ A as A (8) = (cid:90) Ω ∇ H ( r ) · ˆ η i dV (4) ≈ (cid:90) Ω ∇ θ · ˆ η i dV (6) = (cid:90) Ω |∇ θ | dV = ˜ A (10)We now show that equation (10) is a compatible approximation of A . In particular,we prove that equation (9) is still valid when we replace A by ˜ A , which is defined overthe volumes and thus much more convenient to compute. First, as a preliminary stage,we take the gradient of the transport equation (5) in the pursue of a relation betweenthe marker function and the smoothed surface ∂ ∇ θ∂ t + ∇ (cid:0) ( (cid:126) u · ∇ ) θ (cid:1) = · , · ), which simplifies bi-linear integrals as (cid:82) f gdS = ( f , g ). In addition,concepts such as orthogonality, duality or (skew-)symmetry are naturally expressed inthis framework. Further details can be found in Appendix A. With this in mind, wecan proceed to approximate the left hand side of equation (9) via equation (10), to yieldthe following ddt (cid:90) Γ dA (10) ≈ d ˜ Adt = ddt ( ∇ θ, ˆ η i ) = (cid:32) d ∇ θ dt , ˆ η i (cid:33) + (cid:32) ∇ θ, d ˆ η i dt (cid:33) = (cid:32) ∂ ∇ θ∂ t , ˆ η i (cid:33) + (cid:40)(cid:40)(cid:40)(cid:40)(cid:40)(cid:40) (cid:0) ( (cid:126) u · ∇ ) ∇ θ, ˆ η i (cid:1) + (cid:32) ∇ θ, ∂ ˆ η i ∂ t (cid:33) + (cid:40)(cid:40)(cid:40)(cid:40)(cid:40)(cid:40) (cid:0) ∇ θ, ( (cid:126) u · ∇ )ˆ η i (cid:1) = (cid:32) ∂ ∇ θ∂ t , ˆ η i (cid:33) + (cid:26)(cid:26)(cid:26)(cid:26)(cid:26) (cid:32) ∇ θ, ∂ ˆ η i ∂ t (cid:33) (11) = − (cid:0) ∇ (cid:0)(cid:0) (cid:126) u · ∇ (cid:1) θ (cid:1) , ˆ η i (cid:1) = (cid:0) ∇ · ˆ η i , (cid:126) u · ∇ θ (cid:1) (12)where we exploit the skew-symmetry of the convective operator in the second row ofeq. (12), benefit form ∂ ˆ η i ∂ t ⊥ ∇ θ in the third one and the duality between gradient anddivergence in the last one. 8egarding the approximation of the right hand side of equation (9), we can proceedby including equation (7) and then using equation (10) to move from surface to volumeintegrals as − (cid:90) Γ κ(cid:126) u · ˆ η i dA (7) = (cid:90) Γ ( ∇ · ˆ η i ) (cid:126) u · ˆ η i dA (10) ≈ (cid:0) ∇ · ˆ η i , (cid:126) u · ∇ θ (cid:1) (13)We finally obtain: ddt (cid:90) Γ d ˜ A = (cid:0)(cid:0) (cid:126) u · ∇ (cid:1) θ, ∇ · ˆ η i (cid:1) = − (cid:90) Γ κ(cid:126) u · ˆ η i dA (14)We conclude that, from a continuum point of view, the use of a marker function togetherwith the surface reconstruction strategy stated in equation (10) results in a consistentcapture of both the interface and its geometric features, namely the first variation ofarea equation (9). Notice that this analysis is not exclusive to level sets, but ratherextensible to other interface capturing schemes as far as the surface can be cast into apotential form as in equation (10).
3. Symmetry-preserving discretization of single-phase flows
In an incompressible flow, in the absence of external forces, the net balance of me-chanical energy is due to the viscous term of the Navier-Stokes equation solely. Thisis a relevant property for the simulation of turbulent flows, particularly for the com-putation of DNS. In this section, the well-known finite volume, staggered, symmetry-preserving discretization of Verstappen and Veldman [15, 16] is briefly reviewed. Thissets the ground of the newly developed energy-preserving scheme presented in section4. Assuming constant physical properties, Navier-Stokes equations (1) can be rear-ranged to yield ρ (cid:32) ∂(cid:126) u ∂ t + (cid:0) (cid:126) u · ∇ (cid:1) (cid:126) u (cid:33) = −∇ p + µ ∇ (cid:126) u ∇ · (cid:126) u = The evolution of kinetic energy, E k = (cid:0) (cid:126) u , ρ(cid:126) u (cid:1) , in a single-phase flow is obtainedby taking the inner product of (cid:126) u and equation (15), which, in the absence of external9orces and without contributions from the boundaries, yields: dE k dt = − ρ ( (cid:126) u , ( (cid:126) u · ∇ ) (cid:126) u ) − ( (cid:126) u , ∇ p ) + µ ( (cid:126) u , ∇ (cid:126) u ) = − µ (cid:13)(cid:13)(cid:13) ∇ (cid:126) u (cid:13)(cid:13)(cid:13) ≤ (cid:126) u · ∇ ) = − ( (cid:126) u · ∇ ) ∗ ), thecontribution of this term to kinetic energy is null. Duality of the gradient operator withdivergence (i.e., ∇ ∗ = −∇· ) together with the incompressible constrain of the velocity( ∇ · (cid:126) u =
0) results in a null contribution of the pressure term to kinetic energy [16].Finally, by exploiting again the duality between gradient and divergence in the viscousterm, this results in a negative-definite operator, µ ( (cid:126) u , ∇ (cid:126) u ) = − µ ( ∇ (cid:126) u , ∇ (cid:126) u ) = − µ (cid:13)(cid:13)(cid:13) ∇ (cid:126) u (cid:13)(cid:13)(cid:13) ,which, as expected, dumps kinetic energy. We are now going to review the well-know symmetry-preserving, second-order,staggered, finite volume discretization introduced by Verstappen and Veldman in [15],which was subsequently extended to fourth order in [16], from the algebraic perspectiveby means of the tools introduced in Section 2.1. This lays the foundation of the newlydevelopments introduced in Section 4. The discretization in a staggered grid startsby defining the discrete divergence operator, D , directly from the Gauss-Ostrogradskytheorem (cid:90) Ω ∇ · (cid:126) udV = (cid:90) ∂ Ω (cid:126) u · (cid:126) ndS ≈ − T FC A F u f = M C Du f , (17)where M C ∈ R | C |×| C | stands for the metric of the cells, which is a diagonal matrixcontaining cells volume. T FC ∈ R | C |×| F | takes care of the appropriate sum of fluxes overthe faces and A F ∈ R | F |×| F | is the diagonal matrix containing the surface of all faces.Finally u f stands for the staggered velocities. We can rearrange equation (17) to yield D = − M C − T FC A F , (18)this leads to D ∈ R | C |×| F | , as expected for a staggered grid arrangement. Next, thediscrete gradient operator, G , is constructed to preserve duality( u f , Gp c ) F = − ( Du f , p c ) C , (19)10here ( a c , b c ) C = a Tc M C b c stands for the weighted inner product in the cell space, C .It can be defined conversely for the face space, F . Further details on inner products canbe found in Appendix A.In the context of Navier-Stokes equations, preserving this duality at the discretelevel results into a null contribution of the pressure term to kinetic energy [16] G = − M F − D T M C (20)where M F ∈ R | F |×| F | corresponds to the metric of the face space, and thus the definitionof such a metric induces the proper construction of G . This is nothing but the definitionof the staggered control volume. Notice that G ∈ R | F |×| C | . Again, this locates gradientsat faces, as expected for a staggered discretization. M F is defined as M F = ∆ x F A F (21)Note that ∆ x F ∈ R | F |×| F | is the diagonal arrangement of the distance between cellcenters across the face, while A F ∈ R | F |×| F | is also diagonal and contains face surfaces.The final form of G results in G = ( ∆ x F ) − T CF (22)where the standard second-order approximation of the gradient arises naturally fromthe definition of the staggered control volume (i.e., one induces the other).By concatenation, the discretization of the scalar Laplacian operator, L , the essen-tial element of the FSM, can be defined as follows (cid:90) Ω ∇ pdV ≈ M C Lp c = M C DGp c (23)As expected, L ∈ R | C |×| C | . We can expand the final integrated form of the discreteLaplacian as M C L = − T FC A F ( ∆ x F ) − T CF (24)where such a discretization results in a negative-definite operator since A F and ∆ x F arepositive-definite, and T FC = T CFT . This is the ultimate responsible of the di ff usivecharacter of viscosity in the context of Navier-Stokes equations.Finally, the convective term can proceed as in Hicken et al. [37] in order to con-struct a skew-symmetric discretization. Even when dedicated convective operators may11e constructed for Cartesian meshes, this provides with a more flexible approach. Theidea is to construct proper face-to-cell and cell-to-face shift operators in order to exploitthe collocated convective operator as C ( u f ) F = Γ f → c ( I d ⊗ C ( u f ) C ) Γ c → f (25)which guarantees skew-symmetry as far as the vector-valued shift operators are trans-pose (i.e., Γ f → c = Γ T c → f ), as it is the case. Further details can be found in [15, 16, 37, 38]and references therein. By preserving (skew-)symmetries of the operators, as it was described above, theconservation of kinetic energy is guaranteed in the semi-discrete setup (i.e., up to tem-poral integration [39]), mimicking then the continuous behavior of the system. In par-ticular, the semi-discretized energy balance equation reads d E k dt = − ( u f , C ( u f ) F u f ) F − ( u f , Gp c ) F + µ ( u f , L F u f ) F ≤ , (26)which is the discrete counterpart of equation (16). As expected, the only term con-tributing to kinetic energy is the viscous one, i.e., µ ( u f , L F u f ) F , where L F is the standard,negative-definite, staggered di ff usive operator [16]. Note that this holds thanks to thespecific construction of the operators involved and if the incompressibility constrain ofvelocity is satisfied at the discrete level as well (i.e., Du f = c ).
4. Energy-preserving discretization of multiphase flows
Multiphase flows present discontinuities at the interface due to the di ff erence ofphysical properties and the existence of interfacial phenomena, namely, surface ten-sion. This section develops, on top of the symmetry-preserving scheme reviewed inthe previous section, a novel energy-preserving scheme for the discretization of curva-ture. Curvature plays a key role in the development of discontinuities, [ · ], across theinterface as [ σ ] ˆ η i = − γκ ˆ η i (27)12here γ states for the surface tension coe ffi cient, which we assume constant. Thisconfigures the resulting surface tension force, which acts at the interface by imposinga jump condition into the stress tensor which “pulls” the interface towards a lower freeenergy state. The original governing equations (1) can then be reformulated as ρ (cid:32) ∂(cid:126) u ∂ t + ( (cid:126) u · ∇ ) (cid:126) u (cid:33) = −∇ p + ∇ · τ ∇ · (cid:126) u = (cid:2) p (cid:3) = ˆ η Ti [ τ ] ˆ η i − γκ (29)where σ = − p I + τ is split into hydrostatic and deviatoric. The discussion about thediscretization of τ is out of the scope of this work, so the interested reader is refereed toLalanne et al. [26] for a thoughtful discussion on this topic. At this point, it is assumedthat τ presents a prescribed discontinuity at the interface. Regardless of viscous dissipation, incompressible, multiphase flows, do not pre-serve kinetic energy. Instead, surface ( E p = γ A ) and kinetic ( E k = (cid:0) (cid:126) u , ρ(cid:126) u (cid:1) ) energyare exchanged through the pressure term [28] such that, except for the aforementionedviscous term, mechanical energy is conserved. Interface deformation results then in atransfer, through the pressure jump, between surface and kinetic energy. In order toanalyze such a transfer, we start by analyzing the evolution of kinetic energy for mul-tiphase flows, which is obtained by taking the inner product of (cid:126) u and, this time, thegeneral formulation of an incompressible, Newtonian fluid given in equation (28) dE k dt = − (cid:0) (cid:126) u , (cid:0) ρ(cid:126) u · ∇ (cid:1) (cid:126) u (cid:1) − (cid:0) (cid:126) u , ∇ p (cid:1) + (cid:0) (cid:126) u , ∇ · τ (cid:1) (30)As stated in Section 3.1, the skew-symmetry of the convective term results in a nullcontribution to kinetic energy, while the stress term includes an extra contribution dueto the discontinuity at the interface stated in equation (27). − (cid:0) (cid:126) u , ∇ p (cid:1) + (cid:0) (cid:126) u , ∇ · τ (cid:1) = (cid:0) ∇ · (cid:126) u , p (cid:1) − (cid:0) ∇ (cid:126) u , τ (cid:1) − (cid:90) Γ (cid:2) p (cid:3) (cid:126) u · ˆ η i dS + (cid:90) Γ [ τ ] (cid:126) u · ˆ η i dS (31)Further details on the treatment of discontinuities within the inner product can be foundin Appendix A. 13ext, by considering an incompressible flow ( ∇ · (cid:126) u = ∇ (cid:126) u into symmetric ( S ) and skew-symmetric ( W )parts we obtain dE k dt = − (cid:0) (cid:126) u , ∇ p (cid:1) + (cid:0) (cid:126) u , ∇ · τ (cid:1) = − (cid:0) ∇ (cid:126) u , τ (cid:1) − (cid:90) Γ (cid:126) u (cid:2) p (cid:3) ˆ η i dA + (cid:90) Γ (cid:126) u [ τ ] ˆ η i dA = − (cid:0) ∇ (cid:126) u , τ (cid:1) + γ (cid:90) Γ κ(cid:126) u · ˆ η i dA = − ( S + W , µ S ) + γ (cid:90) Γ κ(cid:126) u · ˆ η i dA = − S , µ S ) + γ (cid:90) Γ κ(cid:126) u · ˆ η i dA = − (cid:13)(cid:13)(cid:13) √ µ S (cid:13)(cid:13)(cid:13) + γ (cid:90) Γ κ(cid:126) u · ˆ η i dA (32)As expected, viscosity results in a negative contribution to kinetic energy, whereassurface tension can take any sign depending on whether the interface is expanding orcontracting.On the other hand, the evolution of surface energy is related to the evolution of theinterfacial area. By considering Helmholtz’s free energy [40] dF = γ dA , (33)and plugging it in equation (9), we state that dE p dt = (cid:90) Γ ddt dF = γ (cid:90) Γ ddt dA = − γ (cid:90) Γ κ(cid:126) u · ˆ η i dA (34)Finally, performing a global balance of energy by combining equations (32) and (34),we obtain dE m dt = dE k dt + dE p dt = − (cid:13)(cid:13)(cid:13) √ µ S (cid:13)(cid:13)(cid:13) (35)As expected, surface tension does not play a role in the dissipation of energy, butrather produce a dynamic exchange between kinetic and surface ones of magnitude γ (cid:82) Γ κ(cid:126) u · ˆ η i dA . In the same spirit that symmetry-preserving methods aim at ensuring a null contri-bution of both pressure and convective terms in equation (16) at the discrete level, the14ask in a multiphase flow simulation adds to the requirements to preserve the propertransfer between kinetic and potential energies as d E m dt = d E k dt + d E p dt (36)Namely, if symmetry-preserving schemes were constructed to satisfy at a discrete levelequation (19), energy-preserving methods also satisfy the discrete version of equation(9) in order to properly capture energetic exchanges between kinetic and surface ener-gies. This requires the reformulation of the convective term for variable density flowsto preserve skew-symmetry, as proposed by Rozema et al. [41], which however is outof the scope of this work. Nonetheless, the transfer between kinetic and surface energyoccurs trough the surface tension term as ddt ( G θ c , ˆn f ) F = − ( UG θ c , k f ) F (37)where U = diag ( u f ) ∈ R | F |×| F | is the diagonal arrangement of face velocities, θ c ∈ R | C | is the cell-centered marker function vector and k f ∈ R | F | is the staggered curvaturevector. We consider the advection of the marker function in terms of the discretizedequation (5) d θ c dt = − C ( u f ) C θ c (38)where C ( u f ) C ∈ R | C |×| C | stands for the convective term of the marker function. It mayusually include a high-resolution scheme, as we shall see later, but so far we consider itas a single operator. We disregard the role of recompression stages in time derivativesbut rather consider them as correction steps, which is discussed later on this section.As previously exposed for the continuum case, we can proceed by constructing thediscrete counterpart of equation (14) as (cid:32) G d θ c dt , ˆn f (cid:33) F = − ( UG θ c , Υ D ˆn f ) F (39)where a new shift operator, Υ ∈ R | F |×| C | , is introduced in order to map the curvaturefrom cells to faces. Exploiting the duality of the discrete gradient and divergence oper-ators, equation (19), we obtain − (cid:32) d θ c dt , D ˆn f (cid:33) C = − ( UG θ c , Υ D ˆn f ) F (40)15y subsequently expanding the inner products, we obtain − (cid:32) d θ c dt (cid:33) T M C D ˆn f = − UG θ cT M F Υ D ˆn f ∀ ˆn f (41)which must hold regardless of the interface normal, ˆn f , and consequently independentlyof the cell-centered curvature, D ˆn f . This implies that − (cid:32) d θ c dt (cid:33) T M C = − ( UG θ c ) T M F Υ (42)must hold at any time, while releasing a degree of freedom regarding the definition ofthe normal. We can now plug equation (38) in for the time derivative and expand thetranspose terms − (cid:32) d θ c dt (cid:33) T M C = ( C ( u f ) C θ c ) T M C = θ cT C ( u f ) TC M C = − θ cT G T U T M F Υ ∀ θ c (43)which should hold for any θ c . This leads to C ( u f ) TC M C = − G T U T M F Υ (44)where, exploiting the diagonal arrangement of both U and M F to cast G T U T M F into G T M F U , we can use equation (20) to obtain the final condition as − ( M C C ( u f ) C ) T = M C DU Υ (45)From where we can infer that the convective scheme of the marker function deter-mines the curvature shift operator. This identity guarantees that energy transfers arebalanced and thus total mechanical energy, E m , is preserved up to temporal integra-tion, in the same way that kinetic energy, E k , is preserved in the symmetry-preservingdiscretization presented in Section 3.2 for the single-phase case.Regarding the construction of C ( u f ) C , any high-resolution scheme can be embeddedinto the algebraic form C ( u f ) C = DU Ψ , where Ψ ∈ R | F |×| C | is the actual high-resolutioncell-to-face interpolator. For the CLS, this typically corresponds with SUPERBEE[11]. We can split Ψ as Ψ = Π + Λ [42], to produce C ( u f ) C = DU ( Π + Λ ) (46)16his represents the symmetric ( DU Π ) and skew-symmetric ( DU Λ ) components of C ( u f ) C . The extension to VOF schemes, nicely summarized by Patel et al. [43],requires a previous casting of the advection scheme into the same framework intro-duced in [42]. Plugging equation (46) into equation (45) results in the final form of thededicated cell-to-face interpolation for curvature Υ = Π − Λ (47)which guarantees a proper potential and kinetic energy transfer. An illustrative examplecan be seen in Figure 3. In short, any upwind-like component used for the advectionof θ c turns into a downwind-like component for the interpolation of k f . This can becompared with the second-order midpoint rule used by Olsson and Kreiss where Υ =Π [11]. Ψ( u f ) = c c fu f +1 = Π c c fu f +1 / / c c fu f +1 / − / u f ) = c c fu f +1 = c c fu f +1 / / − c c fu f +1 / − / Figure 3: Example of a particular high-resolution scheme Ψ for the advection of θ c (in this example, thewell-known upwind scheme) and the corresponding dedicated curvature interpolator, Υ . In this case, theinterpolation scheme for curvature is downwind. By mimicking equations (30) and (32) we obtain the discrete counterpart of kineticenergy as d E k dt = γ ( UG θ c , Υ k c ) + µ ( u f , L F u f ) , (48)which assumes a proper discretization of all other terms described in Section 3. Weproceed similarly for potential energy by mimicking equation (34) to define discrete17otential energy as d E p dt = γ (cid:32) G d θ c dt , ˆn f (cid:33) (49)We obtain the semi-discretized total energy equation by combining equations (48)and (49), which, in combination with equation (39) yields d E m dt = d E k dt + d E p dt = γ ( UG θ c , Υ k c ) + µ ( u f , L F u f ) + γ (cid:32) G d θ c dt , ˆn f (cid:33) = µ ( u f , L F u f ) ≤ Du f = c ).The role of interface recompression deserves a special remark. Customarily in-cluded in the level set literature [32, 11], its role is to recover the interface sharp-ness that may have been deteriorated by the convective schemes by taking additional,correcting steps after an initial advection stage. Nonetheless, even when performedconserving mass, as in [11, 12], the nature of recompression results in a non-null con-tribution to potential energy, which violates the conservation of mechanical energy. Forthis reason, the energy-preserving method presented here disregards recompression tofocus on the physical coupling between marker advection and momentum transport.Similarly, other interface capturing schemes may consider additional steps aimed atrecovering interface quality and / or mass conservation [44]. While the results presentedhere allow to adopt this formulation into the momentum equation, including additionalcorrecting steps require an individualized analysis.
5. Results
Equipped with the discretization described in Section 4, we assess its performancefor canonical tests for multiphase flow systems. We focus on inviscid simulations inorder to isolate the performance of our newly developed discretization. Equations (1)18nd (5) are discretized according to the above-mentioned discretization. These read as d u f dt = − C ( u f ) F u f − Gp c + γ K F G θ c (51) d θ c dt = − C ( u f ) C θ c (52)where K F = diag ( Υ k c ) is the diagonal arrangement of the staggered curvature.Density ratio has been fixed to 1 in order to isolate the surface tension term, sim-plify the discretization of the convective term and facilitate the solution of the pressure-velocity decoupling. Nonetheless, as far as the convective term preserves skew-symmetry and the Poisson equation is solved exactly, ratios di ff erent than 1 may beincluded flawlessly. Surface tension forces are included as mentioned in Section 4The system is integrated in time with a second-order Adams-Bashforth schemewhile the pressure-velocity decoupling is achieved with a classical FSM [14]. An ef-ficient FFT decomposition in the periodic direction coupled with a Cholesky solver isused to ensure divergence-free velocity fields to machine accuracy.All simulations are carried on a Ω = [2 H × H ] square domain, where H is boththe semi-width and semi-height of the cavity. Top and bottom faces present periodicboundary conditions, while at the sides no-flux boundary conditions is imposed for themarker function (i.e., ∇ θ · ˆ n wall =
0) while free slip is set for velocity (i.e., (cid:126) u · ˆ n wall = d θ c d τ + D Γ f → c N C ( I C − Θ C ) θ c = DE F G θ c , (53)where τ stands for pseudo-time, Γ f → c is vector-valued shift operator, N C ∈ R | dC |×| C | maps scalars to vector fields aligned with the interface normal, while Θ C = diag ( θ c )and E F = diag ( (cid:15) f ) are the diagonal arrangements of, respectively, θ c and (cid:15) ; where (cid:15) is the face-centered smoothing factor defined in Section 2.2.1. Further details can befound in Olsson and Kreiss [11] for the CLS and in Trias et al. [38] for the constructionof the operators. The classical setup of a zero gravity cylindrical column of liquid is tested in orderto show the impact of the newly proposed method into spurious currents. The sectionof the column is located at the center of the domain and is given a radius of R = . H .Velocity is initially stagnant and that is how it should remain throughout the simulation;however, spurious currents are expected to appear due to errors in the calculation ofcurvature [27]. The initial setup is depicted in Figure 4.Linear perturbation theory provides with the time period for an initially cylindri-cal interface perturbed as r ( θ ) = R + r p cos ( s θ ), where s = , , , . . . correspond-ing to ellipsoidal, triangular of rectangular deformations, respectively [45]. Becauselinear theory predicts perfect equilibrium for both s = s =
1, we arbitrarilyassume an ellipsoidal perturbation (i.e., s =
2) in order to obtain a reference state.The oscillation period can be computed for any s as T = π/ (cid:113) ρ R /γ s ( s −
1) [46],while the characteristic length scale is L = π R , which leads a characteristic speedof c = L / T = (cid:113) γ s ( s − / ρ R , while pressure is referenced to ρ c . Integration iscarried over 5 T .Results in Figure 5 show how the newly proposed method (right column) resultsin an energy stable simulation by counterbalancing the numerical increase in kineticenergy with a decrease of potential energy. This yields to a stagnant situation in whichboth kinetic and potential energy restore their initial values (i.e., E k = E p = igure 4: Initial setup of the marker function for cylindrical column test case in a 128 ×
128 mesh. Contourlines are plotted for ∆ θ = . On the other hand, the standard midpoint rule interpolation for curvature (left column)results in an increase in total energy .In Figure 6 it can be seen how the newly developed curvature interpolation scheme(right) provides, first of all, an order of magnitude smaller oscillations that the standardone produces (left). In addition, there is a dramatic increase in the flow quality withinthe interface, extending the benefits of the high resolution advection scheme for themarker into the velocity field. On the other hand, the use of the standard midpoint rulefor updating curvature pollutes the flow within both phases.It is remarkable how, despite initializing the interface to a theoretical minimum en-ergy situation (i.e., cylindrical cross-section), numerical imbalances when computingcurvature does not reflect such a situation [27]. Nonetheless, the use of an energy-preserving scheme acts in order to keep energy constant, and so counter-balances suchan artificial movement by modifying the curvature accordingly. This results in a robustmethod which eventually is perturbation-proof.In comparison with Figure 5, Figure 7 shows the impact of recompression in bothschemes. As can be seen, the newly developed interpolation method can do little interms of energy, as the recompression stage increases the energy of the system. Actu-ally, we see how the increase in kinetic energy is even higher than in the previous case,with no recompression associated.On the other hand, Figure 8 shows how the velocity field is clearly distorted in21 -0.0004-0.0003-0.0002-0.000100.00010.00020.00030.0004 -4e-05-3e-05-2e-05-1e-0501e-052e-053e-054e-05 -3e-05-2e-05-2e-05-1e-05-5e-0605e-061e-051e-052e-052e-05 E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E k E p E m d E k dtd E p dtd E m dt Figure 5: Energy evolution of the cylindrical section for a pure advection case (i.e., no recompression) withthe standard interpolation of Olsson and Kreiss [11] for the curvature (left) and the newly proposed method(right) in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively.Figure 6: Velocity magnitude and interface location at t = T for a cylindrical column in a 128 ×
128 meshadvected without recompression. Left figure uses the standard midpoint rule while the right one uses thenewly developed energy-preserving one. Contour lines are plotted for ∆ θ = . -0.002-0.001-0.000500.00050.0010.0020.0020.003 -0.00500.0050.010.0150.020.025 -0.001-0.0008-0.0006-0.0004-0.000200.00020.00040.00060.00080.001 E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E k E p E m d E k dtd E p dtd E m dt Figure 7: Energy evolution of cylindrical column with the complete Olsson and Kreiss method [11] with asingle recompression stage (left), and the same method including the modified curvature interpolation (right)in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively. both cases, degrading the solution with respect to the pure advection algorithm oneand two orders of magnitude with respect to the midpoint and the energy-preservinginterpolation schemes, respectively. Noticeably, we still retain, even by including therecompression scheme, a higher quality of the velocity field within the bounded regionfor the newly developed interpolation scheme. The impact of recompression in theoverall quality of the solution is discussed in Section 6.23 igure 8: Velocity magnitude and interface location at t = T for a cylindrical column in a 128 ×
128 meshadvected with a single recompression step. Left figure uses the standard 2nd order shift operator while theright one uses the newly developed energy-preserving one. Contour lines are plotted for ∆ θ = . In order to stretch the previous result to a dynamic equilibrium situation, an el-lipsoidal section is set by distorting the initially cylindrical case. As in the cylindricalwater column, spurious currents may appear, while this time they accompany legitimatecurrents as a result of regions with a moderate non-constant curvature. The ellipse iscentered in the domain and is defined by x = . cos ( α ) and y = . sin ( α ), where α ∈ [0 , π ). Velocity field is initialized at rest and should follow to the oscillation ofthe ellipsoid throughout the simulation. The initial setup is depicted in Figure 9.In the same fashion that in the cylindrical water column described above, linearperturbation theory is employed in order to obtain a reference state. Characteristiclength is set to L = π R , where R = .
3. Time, velocity and pressure scales used arethe same than those for the cylindrical section case.Figure 10 shows how, while the standard midpoint interpolation (left) clearly in-creases the mechanical energy of the system, the newly proposed energy-preservinginterpolation scheme for curvature (right) preserves mechanical energy, which yieldsphysically consistent results and numerically stable simulations. There is, however,both positive and negative o ff sets for kinetic and potential energies. While kinetic and24 igure 9: Initial setup of the marker function for the oscillating ellipse test case in a 128 ×
128 mesh. Contourlines are plotted for ∆ θ = . potential energy are supposed to oscillate between 0 and its maximum or minimum foran ideal harmonic oscillator, we observe that this is not the case. This is explained byan imbalance in the momentum equation, which provides an artificial acceleration inthe fluid, resulting in an increase of the kinetic energy base state [47]. By virtue of theenergy-preserving scheme the oscillation gap for potential energy is reduced accord-ingly, resulting in a decrease of the elongation amplitude. This plays a relevant rolein the next case presented, the capillary wave, which is further discussed in the nextsubsection. Despite this well-known issue results still show the expected oscillatory be-havior of the ellipsoid. This can be checked from the bottom row of Figure 10, wherethe magnitude of the energy transfers remains approximately constant throughout thesimulation. In terms of the oscillating behavior, the increase in mechanical energy forthe naive interpolation results not only in artificially higher values of kinetic energy,but also in a phase di ff erence with respect to the energy-preserving one.Figure 11 presents the marker and velocity fields after t = T with a pure advec-tion scheme. Results for the energy-preserving scheme (right) show a shift in phasewith respect to the midpoint interpolation scheme (left). Velocity is not only higherfor the naive approach, but also the shape of the interface provides with non-physicalcurvature, as it can be observed by the kink appearing along the horizontal centerlineof the ellipsoid (left), which can be compared with the smoother profile present in theenergy-preserving approach (right). 25 -0.06-0.04-0.0200.020.040.06 -0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05 -0.03-0.02-0.01-0.01-0.00500.0050.010.020.020.03 E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E k E p E m d E k dtd E p dtd E m dt Figure 10: Energy evolution of the ellipsoidal section for a pure advection case (i.e., no recompression) withthe standard interpolation of Olsson and Kreiss [11] for the curvature (left) and the newly proposed method(right) in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively. In summary, the use of the energy-preserving scheme provides a higher degree ofreliability, by preserving mechanical energy also in a dynamical equilibrium situation.Despite the numerical errors in which the discretization of momentum may occur, themethod is robust and still preserves mechanical energy.The results obtained by including a single recompression step into the algorithmare presented in Figure 12. They show how, irrespective of the use of an energy-preserving scheme into the advection scheme, the amount of energy included into thesystem in order to keep a sharper profile results in a small, but non-physical, increaseof mechanical energy. Compared with Figure 10, it can be seen how the di ff erence isnot as much in mechanical energy but rather in the nature of the oscillations. Whileresults without recompression still preserve to some extent the oscillating nature of thephysical system, recompression produces an enhanced smoothing, resulting in a flatprofile in terms of both kinetic and potential energy.The claim stated above can be clearly seen in Figure 13, where the initial ellipsoid,26 igure 11: Pressure field and interface location at t = T for the oscillating ellipse in a 128 ×
128 meshadvected without recompression. Left figure uses the standard 2nd order shift operator while the right oneuses the newly developed energy-preserving one. Contour lines are plotted for ∆ θ = . expected to present a dynamical equilibrium, results in a fully rounded shape. Besides,Figure 13 shows how the resulting fields, in both cases, are irrespective of the interpo-lation scheme for curvature used for the advection scheme. Further discussion on theimpact of recompression in the final result is discussed in Section 6.27 -0.03-0.02-0.02-0.01-0.01-0.00500.0050.010.020.020.03 -0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05 -0.03-0.02-0.01-0.01-0.00500.0050.010.020.020.03 E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E / E t/ π p ρR /s ( s − γ d EE / d t t/ π p ρR /s ( s − γ E k E p E m d E k dtd E p dtd E m dt Figure 12: Energy evolution of the oscillating ellipse with the complete Olsson and Kreiss method [11]with a single recompression stage (left), and the same method including the modified curvature interpolation(right) in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively.Figure 13: Pressure field and interface location at t = T for the oscillating ellipse in a 128 ×
128 meshadvected with a single recompression step. Left figure uses the standard 2nd order shift operator while theright one uses the newly developed energy-preserving one. Contour lines are plotted for ∆ θ = . .3. Capillary wave A pure capillary wave is set by originally locating the interface at x = . sin ( ky ),producing an initial wave along the vertical center line of wavelength 2 π/ k . We set k = π/ H , so that a single oscillating period is contained within the domain. Velocityis initially at rest. With the mentioned boundary and initial conditions, the wave isexpected to oscillate indefinitely, alternating states of maximum potential energy (i.e.,maximum elongation) and minimum kinetic energy (i.e., fluid at rest) and vice-versa.Initial setup is presented in Figure 14.As is well known from linear perturbation theory [45], the oscillation of the givensetup present a characteristic period of T = π (cid:112) ρ/γ k tanh ( kH ), which is used as thereference value for time. On the other hand, the characteristic length scale is L = π/ k ,the wavelength of the perturbation. This yields a characteristic velocity of c = L / T = (cid:112) γ k tanh ( kH ) / ρ , while pressure is referenced to ρ c , where ρ stands for the average.Integration in time is set to 2 T . Figure 14: Initial setup of the marker function in a 128 ×
128 mesh. Contour lines are plotted for ∆ θ = . Results in Figure 15 show how the energy preserving discretization proposed in thepresent work preserves mechanical energy (top row, solid line) by balancing the result-ing energy transfers (bottom row, solid line). While the standard midpoint interpolationof curvature results in a non-physical increase of mechanical energy, which ultimatelyleads to instabilities, the novel proposed method provides a stable discretization.Even when mechanical energy is conserved in the newly proposed method, both29 . . -0.03-0.02-0.0100.010.020.03 . . -0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 . . -0.03-0.02-0.0100.010.020.03 . . E / E t/ π p ρ/γk tanh ( kH ) d EE / d t t/ π p ρ/γk tanh ( kH ) E / E t/ π p ρ/γk tanh ( kH ) d EE / d t t/ π p ρ/γk tanh ( kH ) E k E p E m d E k dtd E p dtd E m dt Figure 15: Energy evolution of the capillary wave for a pure advection case (i.e., no recompression) withthe standard interpolation of Olsson and Kreiss [11] for the curvature (left) and the newly proposed method(right) in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively. the amplitude of kinetic and potential oscillations (Figure 15, top row, right) and themagnitude of the energy transfers (Figure 15, bottom row, right) exhibit a significantdamping. The reason behind such a damping is the non-null contribution of surface ten-sion to the momentum equation (the desired result for a closed surface) which producesan artificial acceleration of the fluid. The origin of such artificial acceleration lies in thediscretization of curvature, particularly the computation of normals, which is at the ori-gin of the errors that propagate to the momentum equation. This non-physical increasein kinetic energy manifests itself as an increase of the base level of kinetic energy ato ff -peaks, as can be seen in the top row of Figure 15. While naive interpolation tech-niques are unresponsive to such energy increments, the new energy-preserving methodadjusts the transfers between kinetic and potential energies through surface tension tokeep mechanical energy constant. As a result, the artificial and progressive increase inthe kinetic energy level leaves no room to capillary oscillations, driving the system toa stagnant, but stable, situation. 30 igure 16: Pressure field and interface location at t = T for a pure capillary wave in a 128 ×
128 meshadvected without recompression. Left figure uses the standard 2nd order shift operator while the right oneuses the newly developed energy-preserving one. Contour lines are plotted for ∆ θ = . Figure 17, on the other hand, includes a recompression step into the evolution ofthe wave. Results show clearly how, despite its known advantages [11], the resultingsolution does not preserve energetic balances but rather increase total energy of thesystem, leading to eventual instabilities. It can be seen how the gain in sharpnessintroduced by recompression schemes is at the expenses of destroying the advantagesof the energy-preserving discretization. Results in Figure 18 can be compared withthose of Figure 16, which shows how recompression increases the total energy of thesystem. Namely, the scale in Figure 18 shows how velocity magnitudes are clearlyhigher regardless of the advective step is energy-preserving or not. Among them, theenergy-preserving scheme shows milder velocity fields. This role of recompression isanalyzed in Section 6. 31 . . -0.03-0.02-0.0100.010.020.030.04 . . -0.1-0.0500.050.10.150.2 . . -0.2-0.2-0.1-0.0500.050.10.10.2 . . E / E t/ π p ρ/γk tanh ( kH ) d EE / d t t/ π p ρ/γk tanh ( kH ) E / E t/ π p ρ/γk tanh ( kH ) d EE / d t t/ π p ρ/γk tanh ( kH ) E k E p E m d E k dtd E p dtd E m dt Figure 17: Energy evolution of the capillary wave with the complete Olsson and Kreiss method [11] with asingle recompression stage (left), and the same method including the modified curvature interpolation (right)in a 128 ×
128 mesh. Top rows show the discrete values of kinetic ( E k ), potential ( E p ) and total ( E m )energy. Bottom rows show their semi-discretized time derivative according to equations (48), (49) and (50),respectively.Figure 18: Pressure field and interface location at t = T for a capillary wave in a 128 ×
128 mesh advectedwith a single recompression step. Left figure uses the standard 2nd order shift operator while the right oneuses the newly developed energy-preserving one. Contour lines are plotted for ∆ θ = . . Concluding Remarks By incorporating the first variation of area, equation (9), into the continuum formu-lation we have explicitly imposed a novel condition to the system. Equation (14) showsthat the use of a smooth marker function is compatible with such a condition. Thiscondition is implicitly incorporated into the discretized system by means of a newlydeveloped curvature cell-to-face shift operator, Υ , defined in equation (47). Analyticaland numerical assessments provide evidence that, in the absence of recompression, thenovel interpolation scheme preserves mechanical energy up to temporal integration bybalancing kinetic and potential energy transfers to machine accuracy.The exact value of both kinetic and potential energy is not achieved due to the lackof conservation of linear momentum. This implies that, while the transfers betweensurface and kinetic energy are equal and of opposite sign, its magnitude is not neces-sarily the correct one..In this regard, the adoption of a fully conservative momentum formulation, alongwith proper discretization techniques for the convective operator, as already announcedin Section 4.2, should be considered in a general case. However, the formulation of thesurface tension is the most challenging term. Not being cast into a conservative form, itrelies on the accurate capturing of the interface to produce a closed, and thus conserva-tive, force field. In summary, the use of a finite grid prevents us from resolving the finestscales of the interface, represented by the marker function θ . This under-resolution of θ , either induced both by the mesh and the advection scheme, induces subsequent er-rors in the computation of both ˆ η i and κ , as stated by Magnini et al. [27]. These errorsspread into the momentum equation, which can be seen as a back-scatter of energyfrom the finest, unresolved, surface representation scales into larger kinematic ones,manifesting itself as an inappropriate momentum balance, which ultimately leads to aninaccurate kinetic energy level. This is a well-known issue in multiphase flows and theobject of ongoing research [48, 49].Nevertheless, despite the lack of linear momentum conservation, mechanical en-ergy is conserved and thus the stability of the system is guaranteed up to temporalintegration. From this perspective, the novel technique may provide extra reliability33or surface energy governed phenomena, particularly those involving surface break-upor coalesce, as it may occur in atomization processes or Plateau-Rayleigh instabilities,among others.Recompression schemes, despite producing an energetic imbalance, as has beenshown in Section 5, are common in the level-set community. They preserve a coherentmarker field at the expenses of introducing non-physical energy to the system. Evenwhen the proposed method enforces the energetic consistency between marker and mo-mentum transport equations, the inclusion of recompression prevents us from obtaininga fully energy-preserving scheme. Following the spirit described in Section 4, a firstapproach may be to modify the recompression step to produce not only a conservative,but an energetically neutral resharpening. Enforcing a null contribution to potentialenergy of equation (53), if possible, would allow an arbitrary number of recompres-sion steps, avoiding any penalty in terms of energetic balances. Although this wouldbe desirable, it requires to re-formulate a mass- and energy- conservative recompres-sion scheme which e ff ectively moves the interface irrespective of advection, which isdefinitively not obvious.Others have tried to include recompression within the advection step to yield asingle-step method. After all, recompression is included to fix the distortion producedby interface advection. This leads to phase-field-like methods [50, 51]. Interpretingthis idea as a custom-made high resolution scheme, these approaches can eventuallybe cast into a convective form like that in equations (46) and proceed to obtain theequivalent curvature interpolation as in equation (47). A variant of this model may beto approach the advection of the marker function as a regularization problem [52].Lastly, both a review of the well-known symmetry-preserving scheme and the de-velopment of the energy-preserving scheme have been approached from an algebraicpoint of view. Aside from the advantages in terms of algebraic analysis, the use ofan algebra-based discretization provides an opportunity for High Performance Com-puting (HPC) optimization, parallelization and portability [53]. By casting di ff erentialforms into algebraic ones, (i.e., matrices and vectors), it has been shown in [53] thatnearly 90% of the operations comprised in a typical FSM algorithm for the solutionof incompressible Navier-Stokes equations can be reduced to Sparse Matrix-Vector34ultiplication (SpMV), generalized vector addition (AXPY) and dot product (DOT).In this regard, the present formulation falls within a smart strategy towards portable,heterogeneous, HPC.
7. Acknowledgements
This work has been financially supported by the Ministerio de Economía y Compet-itividad, Spain (ENE2017-88697-R and ENE2015-70672-P) as well as an FI AGAUR-Generalitat de Catalunya fellowship (2017FI_B_00616) and a Ramón y Cajal postdoc-toral contract (RYC-2012-11996). The authors are grateful to Prof. Roel Verstappen,Prof. Arthur Veldman and MSc Ronald Remmerswaal for their enriching discussions.
Appendix A. Inner products
Inner products are bilinear maps from a vector space to its base field (i.e., ( · , · ) S : S × S → K ). Inner products can be defined over both continuum and discrete spacesas ( f , g ) S = (cid:90) S f gdS ∀ f , g ∈ S (A.1)This definition can be readily applied to discrete fields, yielding the definition ofinner products for discrete vectors within metric spaces as( f s , g s ) S = f s T M S g s (A.2)where M S takes over the role of integrating in space, whereas the transpose of the firstelement provides with the appropriate order to perform the subsequent products andsums. This can be seen by expressing f and g as a finite sum of piecewise defined basefunctions.Within this framework, we can define skew-symmetry as the property of operatorssatisfying ( φ, A ψ ) = − ( A φ, ψ ) ∀ φ, ψ ∈ S A : S → S (A.3)where, in the discrete setting, A must be a skew-symmetric matrix. Similarly, we candefine duality as( φ, A ψ ) = ( A ∗ φ, ψ ) ∀ φ ∈ S ψ ∈ T A : T → S A ∗ : S → T (A.4)35y using the aforementioned definitions and the well-known Gauss-Ostrogradskytheorem, it provides with (cid:0) f , ∇ · (cid:126) g (cid:1) = (cid:90) Ω f ∇ · (cid:126) g = − (cid:90) Ω ∇ f · (cid:126) g + (cid:90) ∂ Ω f (cid:126) g · ˆ n (A.5)where we can see that, assuming that there are no contributions from the boundary, theusual relation (cid:0) f , ∇ · (cid:126) g (cid:1) = − (cid:0) ∇ f , (cid:126) g (cid:1) holds as usual. However, if there is a discontinuityin either f or (cid:126) g as a consequence of, say, an interface ( Γ ) in the domain, this preventsus from using equation (A.5) directly, but rather first in both sub-domains separatelyand then sum them together. This results in an explicit expression as (cid:0) u , ∇ · (cid:126) v (cid:1) = (cid:90) Ω u ∇ · (cid:126) v = (cid:90) ∂ Ω u (cid:126) v · ˆ n + (cid:90) Γ (cid:2) u (cid:126) v (cid:3) · ˆ n − (cid:90) Ω ∇ u · (cid:126) v (A.6)where the discontinuity is now explicitly included in the system. Note then that for adiscrete system, the aforementioned gradient-divergence duality is (cid:16) u s , D (cid:126) v s (cid:17) S = − (cid:16) Gu s , (cid:126) v s (cid:17) S + (cid:90) Γ (cid:104) u s (cid:126) v s (cid:105) (A.7)where the extra rightmost term captures the corresponding jump of the variables underconsideration. Note that a proper approximation of Γ is required in order to obtainaccurate solutions. References [1] S. Tanguy, T. Ménard, A. Berlemont, A Level Set Method for vaporizing two-phase flows, J. Comput. Phys. 221 (2007) 837–853.[2] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set / ghostfluid method for simulating turbulent atomization, J. Comput. Phys. 227 (2008)8395–8416.[3] B. P. Van Poppel, O. Desjardins, J. W. Daily, A ghost fluid, level set methodologyfor simulating multiphase electrohydrodynamic flows with application to liquidfuel injection, J. Comput. Phys. 229 (2010) 7977–7996.364] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin, Benchmarks and numer-ical methods for the simulation of boiling flows, J. Comput. Phys. 264 (2014)1–22.[5] M. Lepilliez, E. R. Popescu, F. Gibou, S. Tanguy, On two-phase flow solvers inirregular domains with contact line, J. Comput. Phys. 321 (2016) 1217–1251.[6] E. Gutiérrez, N. Balcázar, E. Bartrons, J. Rigola, Numerical study of Taylorbubbles rising in a stagnant liquid using a level-set / moving-mesh method, Chem.Eng. Sci. 164 (2017) 158–177.[7] C. W. Hirt, B. D. Nichols, Volume of fluid (VOF) method for the dynamics offree boundaries, J. Comput. Phys. 39 (1981) 201–225.[8] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: Al-gorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988)12–49.[9] D. M. Anderson, G. B. McFadden, Di ff use-Interface Methods in fluid mechanics,C’Est Une Rev. (1997).[10] B. S. Mirjalili, S. S. Jain, M. S. Dodd, Interface-capturing methods for two-phaseflows : An overview and recent developments, Annu. Res. Briefs (2017) 117–135.[11] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, J.Comput. Phys. 210 (2005) 225–246.[12] E. Olsson, G. Kreiss, S. Zahedi, A conservative level set method for two phaseflow II, J. Comput. Phys. 225 (2007) 785–807.[13] N. Balcázar, L. Jofre, O. Lehmkuhl, J. Castro, J. Rigola, A finite-volume / level-set method for simulating two-phase flows on unstructured grids, Int. J. Multiph.Flow 64 (2014) 55–72.[14] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput.22 (1968) 745–762. 3715] R. Verstappen, A. Veldman, Direct numerical simulation of turbulence at lowercosts, J. Eng. Math. (1997) 143–159.[16] R. Verstappen, A. Veldman, Symmetry-preserving discretization of turbulentflow, J. Comput. Phys. 187 (2003) 343–368.[17] K. Lipnikov, G. Manzini, M. Shashkov, Mimetic finite di ff erence method, J.Comput. Phys. 257 (2014) 1163–1227.[18] N. Robidoux, S. Steinberg, A discrete vector calculus in tensor grids, Comput.Methods Appl. Math. 11 (2011) 23–66.[19] E. Tonti, Why starting from di ff erential equations for computational physics?, J.Comput. Phys. 257 (2014) 1260–1290.[20] J. B. Perot, Discrete Conservation Properties of Unstructured Mesh Schemes,Annu. Rev. Fluid Mech. 43 (2011) 299–318.[21] F. X. Trias, A. Gorobets, A. Oliva, Turbulent flow around a square cylinder atReynolds number 22,000: A DNS study, Comput. Fluids 123 (2015) 87–98.[22] L. Paniagua, O. Lehmkuhl, C. Oliet, C. D. Pérez-Segarra, Large eddy simulations(LES) on the flow and heat transfer in a wall-bounded pin matrix, Numer. HeatTransf. Part B Fundam. 65 (2014) 103–128.[23] H. Giráldez, C. D. Pérez Segarra, C. Oliet, A. Oliva, Heat and moisture insulationby means of air curtains: Application to refrigerated chambers, Int. J. Refrig. 68(2016) 1–14.[24] D. Fuster, An energy preserving formulation for the simulation of multiphaseturbulent flows, J. Comput. Phys. 235 (2013) 114–128.[25] M. Sussman, K. M. Smith, M. Y. Hussaini, M. Ohta, R. Zhi-Wei, A sharp inter-face method for incompressible two-phase flows, J. Comput. Phys. 221 (2007)469–505. 3826] B. Lalanne, L. R. Villegas, S. Tanguy, F. Risso, On the computation of viscousterms for incompressible two-phase flows with Level Set / Ghost Fluid Method, J.Comput. Phys. 301 (2015) 289–307.[27] M. Magnini, B. Pulvirenti, J. R. Thome, Characterization of the velocity fieldsgenerated by flow initialization in the CFD simulation of multiphase flows, Appl.Math. Model. 40 (2016) 6811–6830.[28] D. Jacqmin, Calculation of Two-Phase Navier–Stokes Flows Using Phase-FieldModeling, J. Comput. Phys. 155 (1999) 96–127.[29] D. Jamet, D. Torres, J. U. Brackbill, On the Theory and Computation of SurfaceTension: The Elimination of Parasitic Currents through Energy Conservation inthe Second-Gradient Method, J. Comput. Phys. 182 (2002) 262–276.[30] D. Jamet, C. Misbah, Toward a thermodynamically consistent picture of thephase-field model of vesicles: Curvature energy, Phys. Rev. E - Stat. Nonlin-ear, Soft Matter Phys. 78 (2008) 1–8.[31] G. Tryggvason, B. Bunner, a. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber,J. Han, S. Nas, Y.-J. Jan, A Front-Tracking Method for the Computations ofMultiphase Flow, J. Comput. Phys. 169 (2001) 708–759.[32] M. Sussman, P. Smereka, S. Osher, A Level Set Approach for Computing Solu-tions to Incompressible Two-Phase Flow, J. Comput. Phys. 114 (1994) 146–159.[33] S. Chen, B. Merriman, S. Osher, P. Smereka, A Simple Level Set Method forSolving Stefan Problems, J. Comput. Phys. 135 (1997) 8–29.[34] T. Wacławczyk, A consistent solution of the re-initialization equation in the con-servative level-set method, J. Comput. Phys. 299 (2015) 487–525.[35] J. Brackbill, D. Kothe, C. Zemach, A continuum method for modeling surfacetension, J. Comput. Phys. 100 (1992) 335–354.3936] T. Frankel, The Geometry of Physics, 3 ed., Cambridge UniversityPress, Cambridge, 2012. URL: http://ebooks.cambridge.org/ref/id/CBO9781139061377 . doi: .[37] J. E. Hicken, F. E. Ham, J. Militzer, M. Koksal, A shift transformation for fullyconservative methods: Turbulence simulation on complex, unstructured grids, J.Comput. Phys. 208 (2005) 704–734.[38] F. X. Trias, O. Lehmkuhl, A. Oliva, C. D. Pérez-Segarra, R. W. C. P. Verstappen,Symmetry-preserving discretization of Navier-Stokes equations on collocated un-structured grids, J. Comput. Phys. 258 (2014) 246–267.[39] F. Capuano, G. Coppola, L. Rández, L. de Luca, Explicit Runge–Kutta schemesfor incompressible flow with improved energy-conservation properties, J. Com-put. Phys. 328 (2017) 86–94.[40] J. T. Davies, E. K. Rideal, Interfacial Phenomena, Academic Press Inc., London,1961. doi: .[41] W. Rozema, J. C. Kok, R. W. Verstappen, A. E. Veldman, A symmetry-preservingdiscretisation and regularisation model for compressible flow with application toturbulent channel flow, J. Turbul. 15 (2014) 386–410.[42] N. Valle, X. Álvarez, F. X. Trias, J. Castro, A. Oliva, Algebraic implementationof a flux limiter for heterogeneous computing, in: Tenth Int. Conf. Comput. FluidDyn., Barcelona, July, 2018.[43] J. K. Patel, G. Natarajan, A generic framework for design of interface capturingschemes for multi-fluid flows, Comput. Fluids 106 (2015) 108–118.[44] O. Ubbink, R. Issa, A Method for Capturing Sharp Fluid Interfaces on ArbitraryMeshes, J. Comput. Phys. 153 (1999) 26–50.[45] S. H. Lamb, Hydrodynamics, Dover Publications, Inc., New York, NY, 1945.[46] D. Fyfe, E. Oran, M. Fritts, Surface tension and viscosity with lagrangian hydro-dynamics on a triangular mesh, J. Comput. Phys. 76 (1988) 349–384.4047] S. Popinet, Numerical Models of Surface Tension, Annu. Rev. Fluid Mech. 50(2017) 49–75.[48] J. Kim, A continuous surface tension force formulation for di ff use-interface mod-els, J. Comput. Phys. 204 (2005) 784–804.[49] M. O. Abu-Al-Saud, S. Popinet, H. A. Tchelepi, A conservative and well-balanced surface tension model, J. Comput. Phys. 371 (2018) 896–913.[50] J. L. Guermond, M. Q. de Luna, T. Thompson, An conservative anti-di ff usiontechnique for the level set method, J. Comput. Appl. Math. 321 (2017) 448–468.[51] S. Mirjalili, C. B. Ivey, A. Mani, A conservative di ffff