On a marching level-set method for extended discontinuous Galerkin methods for incompressible two-phase flows

In this work a solver for instationary two-phase flows on the basis of the extended Discontinuous Galerkin (extended DG/XDG) method is presented. The XDG method adapts the approximation space conformal to the position of the interface. This allows a sub-cell accurate representation of the incompressible Navier-Stokes equations in their sharp interface formulation. The interface is described as the zero set of a signed-distance level-set function and discretized by a standard DG method. For the interface, resp. level-set, evolution an extension velocity field is used and a two-staged algorithm is presented for its construction on a narrow-band. On the cut-cells a monolithic elliptic extension velocity method is adapted and a fast-marching procedure on the neighboring cells. The spatial discretization is based on a symmetric interior penalty method and for the temporal discretization a moving interface approach is adapted. A cell agglomeration technique is utilized for handling small cut-cells and topology changes during the interface motion. The method is validated against a wide range of typical two-phase surface tension driven flow phenomena including capillary waves, an oscillating droplet and the rising bubble benchmark.


Introduction
Within the past decades the significance of direct numerical simulations in context of multi-phase flows is increasing, not at least due to the increase of computational power.Up to date incompressible flow problems described by the Navier-Stokes equations (NSE) only allow analytical solutions for reduced or simplified settings.Considering in particular immiscible two-phase flows the complexity is even increased due to: discontinuous fluid properties across the interface, low-regularity solutions for flow properties and the presence of interfacial forces.Thus, sophisticated and specialized numerical methods are in need for the direct simulation of such problems.
Discontinuous Galerkin (DG) methods have gained quite some interest in the context of computational fluid dynamics in recent years.Reasons for that is the combination of favorable numerical properties known form other well established methods as the finite volume method (FVM) and the finite element method (FEM).DG discretizations are locally conservative due to the flux formulation and enables a high-order method even on unstructured meshes by adapting the local ansatz functions.This cell locality is especially favorable for local mesh refinement and results in local mass matrices and sparse operator matrices.However, a major drawback of DG methods is the notably increase of degrees of freedom (DOF), when comparing to FVM on the same grid and interpolation order.
Considering the discretization of the sharp interface two-phase NSE with a DG method, the convergence order degenerates for such low-regularity solutions (i.e.kinks on the velocity field and jumps in the pressure field) if the approximation space is not conformal to the interface position.However, in context of transient flow problems with moving interfaces one wants to avoid costly remeshing every time step in order to adapt the solution space.A first method that enriches the approximation space such that the solution may exhibit discontinuities inside a finite element is presented in Moës et al. [1] introducing the extended FEM (XFEM) for the simulation of crack growth in solid mechanics.The first application to incompressible two-phase flows is found in Groß and Reusken [2], where the pressure is discretized with XFEM and exhibits a jump due to the surface tension force.A discretization of both pressure and velocity with XFEM is provided by Fries [3], which additionally allows the representation of the kink in the velocity field.
The first extended method for DG is presented in Bastian and Engwer [4] for the discretization of elliptic scalar problems.In Heimann et al. [5] this approach is applied to incompressible Navier-Stokes two-phase flows.The unfitted DG (UDG) method is based on the non-symmetric interior penalty method, where the cut out meshes for both subdomains are based on a piece-wise linear approximation of the interface.The XDG method proposed by Kummer [6] for steady two-phase flows uses a high-order approximation of the interface in combination with a quadrature technique (i.e. the Hierarchical Moment Fitting [7]) for the implicitly defined domains.The discretization is based on the symmetric interior penalty (SIP) method [8] and the stabilization against small cut-cells is ensured by cell agglomeration.Another approach by Saye [9] is the use of implicitly defined meshes with curved elements that are interface-conforming.The implicit mesh DG method provides high-order accuracy for the interface jump conditions in combination with the interfacial gauge method [10], where the numerical coupling between the fluid velocity, pressure, and the interface position is reduced.
In this work we follow the XDG method according to Kummer [6].For the extension to the transient two-phase problem, we apply the XDG adapted moving interface approach introduced by Kummer et al. [11].This temporal discretization allows for high orders of accuracy in time without a full discretization on space-time elements as done e.g. in Lehrenfeld [12] and in Cheng and Fries [13] for XFEM.The moving interface approach allows to keep the spatial discretization unchanged and to obtain a conservative scheme in time.The conservation property is not given for splitting approaches such as the Strang splitting used in Heimann et al. [5].In contrast to Kummer [6] the quadrature method proposed by Saye [14] is used.Furthermore, the numerical treatment of the surface tension force is in this work handled via the Laplace-Beltrami formulation.
For the representation of the interface we employ a signed-distance levelset field.The numerical discretization is done by a standard DG formulation resulting in a high-order and sub-cell accurate approximation of the interface.For the evolution of the level-set field a extension velocity field is used in order to preserve the signed-distance property and reduce the need for reinitialization.For the construction of such an extension velocity field, one can employ direct approaches [15,16,17], fast-marching [18]/fast-sweeping methods [19] or PDEbased techniques [20,21].In this work we present a two-staged construction with an high-order elliptic PDE-based approach on the interface cells and a low-order fast-marching on the neighboring cells.
The range of level-set reinitialization methods is similar to the one for the extension velocity problem.It includes direct geometry-based methods [22,14], iterative fast-marching [23]/fast-sweeping [24] methods and PDE-based approaches [25,26,27].For our method we apply again an elliptic PDE-based reinitialization [28] on the interface cells and perform a fast-marching procedure on the neighboring cells.
The outline of this work is as follows.In Section 2 the continuous problem statement for the transient incompressible two-phase flow for immiscible fluids given.An introduction to the XDG space is provided in Section 3. The numerical representation of the level-set field and its evolution by the two-staged evolution algorithm is considered in Section 4. Here, we differentiate between the interface evolution in context of the XDG method (Section 4.1) and the fastmarching algorithm for the level-set field (Section 4.2).The XDG discretization of the considered problem is given in Section 5. Furthermore, we briefly outline the implemented cell-agglomeration technique in terms of graph theory.The over-all solver structure and the adaptive local mesh refinement are described in Section 6.The presented solver is validated against a range of typical twophase surface tension driven flows in Section 7. We conclude this work in Section 8.
2 Problem statement -Transient incompressible two-phase flow for immiscible fluids with a dividing sharp interface We consider the two-phase setting within a sharp interface formulation and define the computational domain Ω ⊂ R 2 as the disjoint partitioning of the time-dependent fluid bulk phases A(t) and B(t) and the moving interface I(t) by We restrict ourselves to the two-dimensional case and assume that the sharp interface where u = u(x, t) is the velocity vector, p = p(x, t) the pressure and g = g(x, t) describes a volume force, e.g.gravity.The physical properties density ρ = ρ(x) and dynamic viscosity µ = µ(x) are piece-wise constant defined by The corresponding jump conditions for the NSE (2) at a material interface I are given by with a given function u D = u D (x, t) on ∂Ω D .In order to close the initial value problem we set an initial condition for the velocity field by where the initial interface position I(0) is given.The material interface evolves according to the bulk velocity u(x, t) at x ∈ I(t).

The extended Discontinuous Galerkin space
In this section the extended Discontinuous Galerkin (extended DG/XDG) space according to Kummer [6] is briefly introduced.In this case the interface is represented by a sufficiently smooth level-set function ϕ(x, t) that is almost everywhere C 1 (Ω)-continuous, but at least ϕ(x, t) ∈ C 0 (Ω).Thus, the partitioning of the computational domain Ω in ( 1) is implicitly defined by Such a representation of the interface via a level-set function allows the direct computation of the interface normal n I by The numerical discretization of ϕ and its evolution in context of the XDG method are discussed in the next Section 4. In order to formulate the XDG space, first, some basic definitions and operators are introduced.We assume that the computational domain Ω is polygonal and simply connected.Thus, the following basic definitions framework hold.
Definition 1 (basic notations).For a polygonal and simply connected computational domain Ω = A(t) ∪ I(t) ∪ B(t) ∈ R 2 we define: • the numerical background mesh K h = {K 1 , ..., K J } that covers the whole domain Ω = ∪ j K j with non-overlapping cells ( Kj ∩K l 1 dV = 0, l = j), where h denotes the maximum diameter of all cells K j .
• the set of all edges in the mesh is given by Γ(t) := ∪ j ∂K j ∪ I(t).This set can be subdivided into denotes the set of all internal edges, Γ D and Γ N the set of edges imposed with Dirichlet and Neumann boundary conditions, respectively.
• a normal field n Γ on Γ, where it is n Γ = n I on I and represents an outer normal on ∂Ω with n Γ = n ∂Ω .
• the broken gradient ∇ h f defines for f ∈ C 1 (Ω \ Γ) the gradient on the domain Ω \ Γ.According to that, the broken divergence ∇ h • f is defined.
The introduction of the interface I within the two-phase setting allows the definition of cut-cells and a corresponding cut-cell mesh.
Definition 2 (cut-cell and cut-cell mesh).For a numerical background mesh K h = {K 1 , ..., K J } we define the time-dependent cut-cells as for a species s(t) ∈ {A(t), B(t)}.The set of of all cut-cells defines the timedependent cut-cell mesh Note that in a background cell K j where the interface is located we end up with two cut-cells K X j,A and K X j,B , otherwise the original background cell is recovered for the respective species s, i.e.K X j,A (t) = K j and K X j,B (t) = ∅, or K X j,A (t) = ∅ and K X j,B (t) = K j .The XDG method is essentially a DG method applied on the cut-cell mesh.Thus, we define the corresponding XDG space as follows.
Definition 3 (XDG space).The broken cut-polynomial space P X k of total degree k is defined as where P k denotes the standard broken polynomial space of total degree k.
The corresponding local approximation of a field property f j in a background cell K j occupied by both phases is given by where 1 s (x, t) defines the characteristic function for species s(t), i.e. 1 s (x, t) = 1 for x ∈ s(t) and zero everywhere else.Thus, the corresponding cut-polynomial basis functions for s are given by Φ j,l (x)1 s (x, t) = Φ X j,l,s (x, t) ∈ P X k (K h , t).Note that the time-independent local basis functions Φ j,l ∈ P k (K h ) are not altered.In such cut-cells we provide for each phase separate coefficients fj,l,A (t) and fj,l,B (t), also denoted as degrees of freedom (DOF).
In order to formulate an XDG discretization for the considered problem (Section 5.1) we introduce at last the jump and average operator on the edges Γ.
Definition 4 (inner and outer value, jump and average operator).For a field f ∈ C 0 (Ω \ Γ) we define the inner and outer values, f − and f + , at the edges Γ by: Thus, the jump and average operator are defined as: {{f }} :=

Level-set discretization and interface evolution
We define the level-set function by its signed-distance formulation leading to the property that |∇ϕ| = 1.In order to allow a high-order and sub-cell accurate interface representation, the level-set function is approximated by a standard DG space P k (K h ).However, such a discretization is in general discontinuous and not suitable for a spatial discretization in the XDG space, which requires at least that ϕ(x, t) ∈ C 0 (Ω).Thus, we introduce two level-set fields The first one is used for handling the level-set evolution (Section 4.2) and we choose the same polynomial degree as for the velocity field with k (see Eq. ( 29)).The second level-set field is used for the discretization (e.g. the interface normal n I and numerical integration on cut-cells) and is given as the L 2 -projection of ϕ evo constrained by continuity conditions at the cell boundaries (Section 4.1).The polynomial degree is chosen to be k + 1.
The following operations related to the level-set, resp.interface, evolution are only performed on a subset around the interface denoted as the narrow-band K nb .This narrow-band K nb := K cc ∪ K near includes the set of cut-cells K cc and its point-neighbors K near , see Figure 1.According to that, both level-set fields ϕ evo and ϕ exhibits the signed-distance property only on the narrow-band.On the far field cells K far the value is either set to ϕ evo = −1 for cells in A or ϕ evo = +1 in B.

Interface evolution algorithm for XDG -Enforcing continuity of the level-set field
The evolution equation for the level-set field ϕ evo in the computational domain Ω is given by where the evolution velocity for a material interface I is set to u evo = u.However, in context of the XDG method the approximation space for the bulk velocity is in P X k (K h ), whereas the level-set evolution needs a velocity field u evo ∈ P k (K h ) according to ϕ evo .In order to generate a velocity field for the advection equation (17), we introduce a density-averaged velocity field u ρAver ∈ P k (K h ), where the velocity components u d j , d = {1, 2}, in each cell K j ∈ K nb are constructed as (18) The coefficients ũd j,l,A and ũd j,l,B are the same coefficients for u ∈ P X k .Recalling the local XDG approximation in (12), we omit the identity function to recover the standard approximation in P k (K h ).A marching procedure is employed for the level-set advection and is presented in Section 4.2.
Enforcing continuity of the level-set field ϕ Since we are dealing with two level-set fields ϕ evo and ϕ, the latter one needs to be constructed from ϕ evo after the advection.Therefore, we employ a L 2 -projection of ϕ evo onto ϕ with additional continuity constraints that is described as the following quadratic optimization problem (OP) with The projection is restricted to the narrow-band K nb and the continuity conditions are defined on the set of internal edges in the narrow-band Γ int,nb .Let be ϕ = φT Φ with Φ ∈ P k+1 (K h ) and φ the sought-after coordinate vector and ϕ * evo = φT evo Φ the projection of ϕ evo ∈ P k (K h ) onto P k+1 (K h ).Thus, the OP is equivalently stated as where M = I, since Φ is a orthonormal basis, d = φT evo φevo and A describes the constraint matrix.The equality constraints (20b) are incorporated via Lagrange multipliers λ following the Karush-Kuhn-Tucker conditions of first order.The solution of the resulting OP can be described as the solution of the following system of linear equations This saddle point system is solved via the Schur complement method, The Schur complement in this case is given by AM −1 A T = AA T and describes a symmetric positive definite matrix.Thus, the saddle point system is reduced to solve AA T λ = A φevo (22) for λ and evaluate the solution vector φ via Note that the problem size of ( 22) is reduced to the number of the constraint conditions.The continuity constraints (20b) at the internal edges Γ int,nb are given as the following equality conditions at a sufficient amount of points pt at Γ i with where φ− i,n and φ+ i,n denote the coefficients of the corresponding inner and outer cell at Γ i .The condition (24) represents one row in the constraint matrix A. The amount of conditions is chosen according to the approximation space of ϕ = P k+1 (K h ).Such a formulation of suitable equality conditions for the continuity projection is also applicable on meshes with hanging nodes, which is necessary for the implementation of adaptive mesh refinement as described in Section 6.
Interface evolution algorithm for a XDG method The entire interface evolution algorithm including the update of the integration domains and the narrow-band is given in Algorithm 1.The actual advection of the level-set field ϕ evo is performed in MarchingLevelSetEvolution (line 4) and is presented in the next section.

A marching level-set advection method
Considering the advection of the level-set field ϕ evo by equation (17) with u evo = u, this would not preserve the desired signed-distance property |∇ϕ evo | = 1.Furthermore, one should note that for the advection of the interface the only valid velocity for the evolution should be the one at the interface itself.Thus, in order to preserve the desired signed distance property, at least sufficiently far from non-smooth regions in ϕ evo and for some finite time, we construct the evolution velocity u evo as an extension velocity field u ext .Therefore, we solve the extension velocity problem (EVP) for each velocity component u d ext with d = {1, 2}: In words, the extension velocity field u ext , resp.its components u d ext , is constructed by propagating the interface velocity u I from I perpendicular into the bulk domain Ω \ I ∈ K nb .
The presented evolution algorithm is divided into two stages.In the first stage the extension velocity is computed monolithic with high-order accuracy on the cut-cells K cc using the reformulation of the EVP (25) into an elliptic PDE.In this case the solution describes the steady-state limit of an anisotropic diffusion problem, see Utz and Kummer [21].In our implementation we follow the work of Utz [29] and normalize the level-set gradient in (25a) with ∇ϕ evo /|∇ϕ evo | leading to more stable results.The discretization is done by the standard SIP method.The corresponding linear problem in u ext is solved by a linear solver.
In the second stage the extension velocity on the near field K near is constructed geometrically in each cell separately using an adapted fast-marching [23,18] procedure, see Algorithm 2. For the low-order local solver in line 5 boundary values are provided on edges that are adjacent to cells in Accepted.The solution is then constructed node-wise from the accepted edges into the cell K l and projected on P k ({K l }).Note that Algorithm 2 is executed twice, once for negative near cells with ϕ evo < 0 and once for ϕ evo > 0.
The entire two-staged level-set evolution algorithm is presented in Algorithm 3. Since the signed-distance property of ϕ evo is only given on the narrow-band, a reinitialization is applied on new cells entering the narrow-band during the last advection.Therefore, a fast-marching procedure (see Algorithm 4) is executed in line 4.After constructing the extension velocity field u ext on the narrowband K nb (line 7 and 8) the new level-set field ϕ n+1 evo is computed according to the evolution equation ( 17) by Algorithm 2 Fast-marching procedure for the EVP on near cells for l = 1, ..., L do 5: solve Eq. (25a) node-wise and project solution on K l 6: end for 9: return u near ext 10: end procedure where Γ int,nb denotes the internal edges in K nb .The stabilization term f upwind describes a simple upwind-flux for the flux ϕ evo u ext .The temporal term is discretized via the total variation diminishing Runge-Kutta scheme of third order according to Gottlieb and Shu [30].At last, the advected level-set field ϕ n+1 evo is additionally stabilized by penalizing jumps at inner edges, i.e.Γ int,nb (η evo /h) ϕ n+1 evo [[ϑ]] dS with η evo = 10, via an implicit Euler scheme with ∆t pnlty = 0.001∆t.Thus, the deviation to the projected continuous level-set field ϕ used for the spatial discretization is reduced.
Algorithm 3 A two-staged marching level-set evolution algorithm on the narrow-band monolithic solve of elliptic formulation of EVP (25) 8: 26) Runge-Kutta scheme (TVD3) return ϕ n+1 * evo 13: end procedure Level-set reinitialization As mentioned before the signed-distance property needs to be initialized for new cells entering the narrow-band (line 4 in Algorithm 3).Such a reinitialization problem is described by the following Eikonal equation Therefore, a fast-marching [23,18] procedure is utilized as given in Algorithm 4.
The local solver for one cell K j in line 5 is based on the elliptic reinitialization method by Basting and Kuzmin [27] and adapted for DG in Utz [28].A low-order initial solution for ϕ evo is obtained node-wise by a direct geometric approach.Furthermore, note that at the end (line 11) of the marching algorithm 3 an optional full reinitialization of the level-set field ϕ evo on the narrow-band may be performed.In this case the elliptic reinitialization [28] is solved monolithic on the cut-cells and the fast-marching procedure 4 is executed on the near field.
Algorithm 4 Fast-marching procedure for the reinitialization problem on near cells while Close = ∅ do return ϕ near evo 11: end procedure 5 The XDG method for the instationary twophase NSE In this section the XDG discretization of the considered transient two-phase flow problem is provided (Section 5.1) and the concept of the cell agglomeration for resolving stability issues and handling topology changes during interface motion is described (Section 5.2).
The formulation of an XDG method requires the numerical integration of cut-cell integrals such as which represent integrals over cut-cell volumes and surfaces and integrals along the interface I.An essential prerequisite of the XDG method is a high-order accurate integration technique allowing for a stable and robust numerical integration on implicitly defined domains.In contrast to Kummer [6] we use in this work the quadrature method proposed by Saye [14].It is generally faster compared to the Hierarchical Moment Fitting, but restricted to hyperrectangles.

XDG discretization of the incompressible two-phase NSE
Before introducing the discretization for the transient two-phase flow problem, we define the function space of the ansatz and test functions for the velocity field u ∈ R 2 and the pressure field p ∈ R by where k = {k, k, k } = k γ with γ = 1, 2, 3 describes the degree vector.Here, the velocity components are discretized by an XDG space of order k and the pressure field of order k = k − 1.Thus, we comply with the Ladyzenskaja-Babuska-Brezzi condition [31,32] for the Stokes system with ρ A = ρ B and Following the spatial discretization in Kummer [6], we propose the following discretization for the transient two-phase incompressible Navier-Stokes equations (2) with the interface jump conditions (4) and boundary conditions (5): The bilinear forms denoted by m n (−, −) represent the mass matrix components for the corresponding time step t n with A third order backward differentiation formula (BDF) is employed, where the coefficients for each time step are given by d 0 = 6/11, d 1 = 18/11, d 2 = −9/11, d 3 = 2/11.One should note that in the context of the moving interface approach by Kummer et al. [11] a spatial approximation of degree k theoretically requires at least a time integration scheme of order 2k for linear equations.BDF schemes are not A-stable above order 2 including more eigenvalues z with (z) < 0 for higher orders.However, the unstable eigenvalues for the third order scheme are comparable small, i.e. near the imaginary axis.Since we are dealing with comparably small time steps due to the capillary time step restriction (see Eq. ( 43)), we choose the use of the third order BDF scheme to comply with the requirement for the moving interface approach.
In the remainder of this section the integration domains of the remaining terms are set to the time-level t n+1 , i.e. the interface is fixed at I = I(t n+1 ).The trilinear form c(u * , u, v) describes the discretization of the convective term, where a local Lax-Friedrichs flux is employed (32) For the choice of the stabilization parameter λ we refer to Kummer [6].Note that the interface integral is excluded from the surface discretization.This contribution is canceled out by which is introduced within the moving interface approach [11].The discretized problem (30) needs to be solved iteratively due to the non-linearity of the convective term.Thus, the whole system is linearized in each iteration with u * .The iteration process in combination with the interface evolution is presented in Section 6.
The bilinear form b(−, −) represents both, the discretization of the pressure gradient and the continuity term, i.e. velocity divergence, An extension to the standard symmetric interior penalty (SIP) method [8] is used for the discretization of the viscous terms.In this case the transposed term ∇u T within the divergence operator is included and since The corresponding penalty parameter η is chosen according to where the local penalty parameter η is computed by The safety factor η 0 is set to η 0 = 4, see Kummer [6] for details.The geometric factor |∂K X s |/|K X s | is directly provided by the construction of the cut-cell quadrature rules, see Eq. (28).
Finally, we specify the term s(v, q) on the right-hand side of the variational formulation (30), which summarizes the Dirichlet boundary conditions and force terms The first two terms in the upper line describe the Dirichlet boundary conditions for the convective part (32) and the viscous parts (35).The first term in the bottom line corresponds to the Dirichlet boundary condition for the continuity equation (34) and the second term denotes the discretization of the volume term ρg.
For the numerical treatment of the surface tension force (remaining two terms in the second line) we do not evaluate the curvature κ, see (4a), but choose the Laplace-Beltrami formulation, see e.g.Gross and Reusken [33], Cheng and Fries [34] and Sauerland and Fries [17].In this case the surface tension force is rewritten to its interface divergence form with σκn I = ∇ I • (σP I ), where P I := I − n I ⊗ n I describes the projection tensor onto the interface I and the interface gradient operator ∇ I is defined as ∇ I := P I ∇.The vector τ denotes the tangential to I and the corresponding integral denotes point-measures in R 2 .Those boundary terms ensure force conservation of the surface tension along closed interfaces.Note that we do not employ a semi-implicit discretization as introduced in Dziuk [35] or a regularization within the continuum surface force framework [36,37,38].

Cell agglomeration
The discretization by the XDG method may introduce arbitrarily small cut-cells K X j,s , where the volume fraction |K X j,s |/|K j | is small compared to the background cell K j .In context of the SIP discretization for the viscous terms (35) such cutcells give rise to undesirably high condition numbers and consequently leads to stability issues.Therefore, cell agglomeration is employed to remove small cut-cells from the discretized system [6].Furthermore, the cell agglomeration is used in the context of the moving interface time discretization [11].In this case, appearing and vanishing cut-cells are handled during the interface motion.
The cell agglomeration and resulting meshes are described in terms of graph theory according to Kummer et al. [39].Basic definitions for the agglomeration are given in the following.An example is provided is in Figure 2.
Definition 5 (logical edge and cell agglomeration).For a cut-cell mesh K X h we define: • a logical edge Edg({K X j,s , K X g,s }) for two neighboring cut-cells K X j,s and K X g,s .
• an undirected graph G(K X h ) := K X h , Edg(K X h ) , where Edg(K X h ) denotes all logical edges in the cut-cell mesh K X h .• a cluster of neighboring cells a = {K X l,s , ..., K X L,s } and a corresponding aggregated cell K X a,s := ∪ ∂ K X l,s ∈a K X l,s .The modified union is defined as All elements in a are related to the same species.
• an aggregation mesh Agg(K X h , A) including all aggregated cells given by an aggregation map A ⊂ Edg(K X h ) and including non-aggregated cells.
Definition 6 (agglomerated XDG space).For some aggregation map A ⊂ Edg(K X h ) and corresponding aggregation mesh Agg(K X h , A) the agglomerated broken cut-polynomial space of total degree k is defined as: Note that P k (Agg(K X h , A)) is a sub-space of the original XDG space P k (K X h ).A corresponding basis to the aggregation mesh Agg(K X h , A) is denoted by Φ X,A .
The interface I crosses a background cell edge during on time step resulting in a vanishing cut-cell K X j,B , i.e. |K X j,B (t 0 )| > 0 and |K X j,B (t 1 )| = 0, and an appearing cut-cell K X g,A , i.e. |K X g,A (t 0 )| = 0 and |K X g,A (t 1 )| > 0. Thus, the agglomeration maps A(t 0 ) = {Edg van } and A(t 1 ) = {Edg new } ensure equal topology on both time steps t 0 and t 1 .
Agglomeration for removing small cut-cells In order to remove small cutcells form the numerical mesh K X h , we introduce an agglomeration map A α ⊂ Edg(K X h ) with edges Edg({K X j,s , K X g,s }) matching the following two conditions: • The volume fraction of K X j,s is smaller than the agglomeration threshold α, i.e 0 < |K X j,s |/|K j | < α.
• The cell K X g,s marks the neighbor with the largest volume fraction in the same species.
In this work the agglomeration threshold α is set to α = 0.1 for all presented simulations.A condition number study with varying agglomeration thresholds is presented in Kummer [6].
Agglomeration for topology changes during interface motion Ensuring the same topology throughout all mass matrices (31) at the different time steps, appearing and vanishing cut-cells due to the interface motion are handled via cell agglomeration, see Figure 3. Therefore, the agglomeration map A α additionally exhibits logical edges Edg({K X j,s , K X g,s }) matching the conditions: • The cut-cell K X j,s appears, i.e. |K X j,s (t n )| = 0 and |K X j,s (t n+1 )| > 0, or vanishes, i.e. |K X j,s (t n )| > 0 and |K X j,s (t n+1 )| = 0.
• The cell K X g,s marks the neighbor with the largest volume fraction in the same species.

Solver structure -Coupling flow solver and interface evolution
In this section the coupling of the flow solver for the discretizied system (30) and the interface evolution is presented.Additionally, the adaptive mesh refinement (AMR) procedure is briefly discussed at the end.In Algorithm 5 the procedure for a single solver run is given.There are two options for the coupling of the interface evolution: explicit and implicit.First, we describe the explicit coupling, where the interface is only updated once at the beginning (Line 3).An initial solution on the new domain Agg(K X h , A n+1 α ) is given by extrapolation of the old solutions u n and p n .The extrapolation operator is defined as follows Note that the coefficients f n do not change.This initial solution further provides a first linearization u * (Line 7) for the convective terms (32).After that, the procedure enters the iterative solution process via Picard iterations until the prescribed convergence criterion NSE is satisfied or a maximum number of iterations i max is exceeded (Line 9).
Within each iteration the linear system of the discretizied Navier-Stokes equations (30) is updated (Line 15).The corresponding matrix formulation is written as with s n+1 = (ũ, p) denoting the sought-after solution vector including the coefficients of the velocity fields ũ and the pressure field p.The corresponding operator matrix Op has the structure of a saddle point problem given by where Op c,a denotes the convective trilinear form c(−, −, −) in ( 32) and viscous bilinear form a(−, −) in (35), and Op b the bilinear forms for the pressure gradient and velocity divergence b(−, −) in (34).After performing the cell agglomeration (Line 16) the resulting system is preconditioned further reducing the condition number (Line 17).Therefore, a block Jacobi preconditioning is used resulting in diagonal matrices containing only 0, −1 and +1 entries for the symmetric part of the block diagonal in Op c,a , see Kummer [6] for details.The preconditioned linear system is solved using the direct solver MUMPS [40,41].
The procedure described above corresponds to an explicit coupling, where the interface position is updated once and is determined by the velocity field of the previous time step, i.e. u evo = u n .Employing an implicit coupling, i.e.

24:
PostProcessing(t n ,∆t) 25: end procedure u evo = u n+1 , the level-set update is additional performed within the iteration process (Line 10).Note that in this case we do not update the level-set in every iteration, but when the Navier-Stokes solution is converged.The coupled iteration process converges when both the Navier-Stokes solution with NSE and the level-set field with LS converge.
Adaptive mesh refinement The XNSE-solver is extended to allow adaptive mesh refinement (AMR) during the simulation.Therefore, a background cell is divided into four equal sized sub-cells.For higher refinement levels we ensure that neighboring cell also refine in order to exhibit always a 2:1 cell ratio on every edge.Thus, we counteract undesired locking-effects on the refined cells by much larger cells.Furthermore, we do not coarsen on cut-cells.The adaption of the mesh and the corresponding approximated DG-fields is performed before RunSolverOneStep and requires a soft-restart of the current simulation.The indication of the current refinement level of each cell is predefined.In general we employ a constant refinement on the narrow band, but additional refinement may be used depending on the flow characteristics.

Numerical results
The following simulations are all computed with an explicit coupling of the level-set evolution and no additional reinitialization is performed.Since approximating the interface, resp.the level-set field, by a polynomial basis of degree k we set the minimal resolved capillary wave length to h k+1 .The resulting capillary time step according to Denner and van Wachem [42] is given by The numerical time step size ∆t always satisfies ∆t < ∆t σ , if not stated otherwise.All simulations in this section are performed with the XNSE-Solver within the open-source solver framework BoSSS [39] developed at the department of fluid dynamics at the TU Darmstadt.The source code is available under the Apache License at http://github.com/FDYdarmstadt/BoSSS.The presented results are available online at https://doi.org/10.25534/tudatalib-327 In this results section we often derive scalar measures over the simulation time from the computed numerical DG fields.In order to quantify the error against a reference solution we define the l 2 -error according to Hysing et al. [43] by where q N describes the scalar measure at time step t n , i.e. q n = q(t n ) with n = 1, ..., NTS and NTS denoting the total number of time steps.q n,ref describes a suitable reference solution at t n , which is either an analytical solution or the solution on the finest grid of a convergence study.If there are less time steps than provided by the reference solution, we use linear interpolation to generate missing data.We determine the experimental order of convergence (EOC) by the slope of a linear regression in the corresponding log-log plot.The regression coefficients are estimated by an ordinary least-squares method.Further, we define the rate of convergence (ROC) for a single refinement level l according to Hysing et al. [43] by In the following we present the results of capillary waves and comparing to the analytical solution for the amplitude height (Section 7.1).Simulations regarding a droplet in equilibrium and non-equilibrium state (Section 7.2).At last we compare our solver against the benchmark groups of the rising bubble test case (Section 7.3).

Capillary wave
We consider the damped oscillations of a capillary wave initiated by a sinusoidal disturbance with a wavelength of λ = 1 and an initial amplitude height of a 0 = 0.01 λ, see Figure 4.The physical properties in both phases are set equal.For such a setting, i.e. small-amplitude limit and same kinematic viscosities ν = µ/ρ, Prosperetti [44] presents an analytical solution of the amplitude height a(t) for the corresponding initial value problem.Comparing to the analytical solution we follow the work of Popinet [45] and set the computational domain . The lower and upper boundaries are imposed by a wall boundary condition and the side boundaries are periodic in x-direction.In order to verify against a representative range of physical regimes, from an overdamped oscillation (low La) to a highly oscillatory behavior (high La), we investigate a study of the Laplace number defined as La = σρλ/µ 2 .The study consists of four values with La = {3, 1.2•10 2 , 3•10 3 , 3•10 5 } and the corresponding physical properties are given in Table 1.
Note that both phases exhibit the same physical properties.All physical settings are computed on three meshes with equidistant mesh sizes of λ/h = {8, 16, 32}.The time step sizes on each mesh are set according to the capillary time step restriction (43) with a spatial discretization of order k = 2.
In Figure 5 the numerical solutions on the finest grid are depicted and compared to the analytical reference solution [44].The reference solution is computed using MATLAB.Throughout the study all numerical results show very good agreement to the analytical solution.However, one should note that the agreement of the damping rate for the most dynamic case La = 3 • 10 5 is slightly underestimated resulting in higher amplitude heights for the last peaks.The frequency is still in very good agreement.
In Table 2 the l 2 -errors and the corresponding EOC value for each study are presented.Looking at the finest solutions of each study one should note the increasing error norm for increasing Laplace number.This is additionally attended by decreasing EOC values.Two more studies for the most dynamic case (La = 3•10 5 ) are given in Table 2.The first one is a mesh study performed with a fixed time step size for all simulations that corresponds to a capillary time step size with an grid size of λ/h = 64.In this case the error norms and the EOC value are slightly improved.For the second study the polynomial degree is set to k = 3 and again the time step size is fixed.As expected the error norms for k = 3 are smaller compared to k = 2, but the study shows a poorer EOC value.This may result from the theoretically insufficient integration scheme, i.e third order BDF, compared to the spatial approximation in context of the moving interface approach [11].Furthermore, one should note that quite small time steps are needed in order to obtain the desired convergence rates.

Droplet simulations
In this section two settings for the transient simulation of a droplet are presented.The first test case considers a droplet in its equilibrium state, i.e. circular shape.For the second one the droplet is initially deformed into an ellipsoidal shape.

Droplet in equilibrium
The equilibrium state, i.e. u = 0, of a droplet in R 2 is described by an circular shape, where the pressure inside a droplet (phase A) of radius r is given by the Young-Laplace equation with However, due to numerical inaccuracies regarding the surface tension force the computed solutions do not provide a zero-velocity field.In order to quantify the discretization error Smolianksi [46] proposed the following setup.The droplet with radius r = 0.25 is set in the middle of the computational domain of Ω = [−0.5,0.5]×[−0.5, 0.5], where all boundaries describe a wall boundary condition.Both phases share the same physical properties with: density ρ = 10 4 , dynamic viscosity µ = 1 and surface tension coefficient σ = 1.This results in a Laplace number of La = σρ2r/µ 2 = 5 • 10 3 .We investigate a mesh study with the following mesh sizes 1/h = {20, 40, 60, 80}.All simulations are performed until t = 125 with a fixed time step size of ∆t = 0.01.Setting the polynomial degree to k = 2 the capillary time step restriction ( 43) is satisfied for all runs.
In Table 3 the L 2 -error norm of the spurious velocities ||u|| L 2 and the error against the exact Young-Laplace solution ||p − p exact || L 2 are given.The norm is computed in the bulk Ω \ I on the terminal time step at t = 125.
One should note that the error norms seem to converge already on the second mesh.This may result from the long simulation time and fixed time step size.However, we take a closer look on some energy related properties.In Figure 6 the kinetic energy in the bulk 1  2 ||ρu • u|| L 2 and the corresponding discrete  .One observes that the numerical solutions on the coarser meshes exhibit a larger initial deviation from the zero-velocity field and that the decay of such spurious velocities takes a longer period of time.The same behavior is shown for the bulk dissipation.One should note that the term is always negative demonstrating the energetic stability of the presented discretization.However, taking the surface divergence, i.e. change rate of the interface area, into account the coarser meshes exhibit an enlarging or diminishing interface area until the end.The finest solution seem to produce a stable interface.
Figure 7: The temporal evolution of the semi-axis in x-direction of an oscillating droplet for t = [0, 1000] (left) and t = [800, 1000] (middle).On the right the temporal evolution of the kinetic bulk energy and surface energy for the finest solution with 1/h = 80.

Oscillating droplet
In this section following the work of Hysing [38] we consider a droplet initialized in an ellipsoidal shape with where the semi-axes are set to a = 1.25r and b = 0.8r.The corresponding equilibrium shape is described by an circle with radius r = 0.25.The computational domain and the physical properties are same as for the equilibrium test case, except for the surface tension coefficient with σ = 0.1.This results in a smaller Laplace number of La ≈ 500.A mesh study on the following mesh sizes 1/h = {10, 20, 40, 60, 80} is performed and the simulations run until t = 100 with a fixed time step size of ∆t = 0.5.Setting the polynomial degree to k = 2, the capillary time step restriction is already exceeded on the second mesh.In Figure 7 on the left the temporal evolution of the semi-axis in x-direction a x is plotted for the whole mesh study.It is remarkable that all simulations show a stable oscillating behavior until the end, see Figure 7 in the middle.Note that on the finest mesh the time step size is 10 times larger than the capillary time step restriction (43).
Taking a closer look on the finest solution, the temporal evolution of the kinetic bulk energy and the surface energy E σ = σ|I| with |I| describing the interface area is plotted in Figure 7 on the right.Note that the plotted surface energy is subtracted by the minimal surface, i.e.E σ,min = σ2πr.Thus, the exchange of kinetic and surface energy is observable in one plot.The maximum values of the surface energy, i.e. largest extend along one semi-axis, corresponds directly to the minimum value of the kinetic energy.Between these states the kinetic energy reaches its maximum always slightly before the minimum value of the surface energy.Due to the small Laplace number the oscillations are strongly damped by viscous dissipation in the bulk.

Rising bubble benchmark
In this last section we compare our results against the rising bubble benchmark test case established by Hysing et al. [43].The initial setup is depicted in Figure 8, where a circular bubble (radius r = 0.25) is positioned at (0.5, 0.5) in the domain Ω = [0, 1] × [0, 2].The lower and upper boundary are imposed by the wall boundary condition.The boundaries at the sides describe a free-slip boundary condition given by where τ ∂Ω denotes the tangent vector on the boundary ∂Ω.The discretization in context of the presented variational problem (30) is given in the Appendix A. The driving force for the rise of the bubble ( ρ A < ρ B ) is the gravity force g = −ge y oriented in negative y-direction.For this initial setting two test cases are defined with different physical properties resulting in an ellipsoidal terminal shape for the first test case and a dimpled cap with filaments for the second test case.The corresponding physical parameters are given in Table 4.For both test cases the simulation is performed for 3 time units.
Benchmark groups and comparing quantities In Hysing et al. [43] an extensive data-set is provided by the following three research groups: TP2D (Transport Phenomena in 2D), FreeLIFE (Free-Surface Library of Finite Element) and MoonMD (Mathematics and object-oriented Numerics in Magdeburg).For details on the methodology of each solver the reader is referred to the benchmark paper [43].The plotted data of the benchmark groups are taken from http://www.featflow.de/en/benchmarks/cfdbenchmarking/bubble/bubble_reference.html.Furthermore, we compare in the following some results to the work of Heimann et al. [5], where the unfitted DG method with a piece-wise linear approximation of the interface is employed.
In order to compare the different numerical methods three scalar measures for the temporal evolution of the rising bubble are introduced [43].First, the center of mass is considered given as The second measure is denoted as circularity and defined by The circularity is ¡ c = 1 for a perfectly circular bubble and gets smaller for more deformed bubble shapes.The third measure is the mean bubble velocity where the component in y-direction is denoted as the rise velocity V c .The error quantification is done via the l 2 -error norm (44) and the corresponding ROC values (45).Note that the in the original benchmark paper furthermore l 1 and l ∞ -error norms are considered.For completeness the definitions and corresponding results are found in the Appendix B.Besides the three scalar quantities the terminal shape of the bubble at t = 3 is compared.

Test case 1 -Ellipsoidal shape
For the first test case two mesh studies were performed for polynomial degrees of k = {2, 3} on meshes with equidistant cell sizes of 1/h = {10, 20, 40, 60, 80}.The corresponding number of background cells, total number of DOF and the number of time steps are given in Table 5.The total number of DOF indicates the maximum number encountered during the simulation.Since the number of cut-cells is depending on the interface form, the number of double DOF is changing according to the number of cut-cells.The time step sizes are chosen for each setting according to the capillary time step restriction (43).
The temporal evolution of the bubble shape at times t = {1.2,2.1, 3.0} is shown in Figure 9. Here, the magnitude of the corresponding velocity field is plotted for the finest solution of k = 3. Comparing the terminal shape at t = 3 with the respective finest solution of benchmark groups (Figure 13 in Appendix B), there are no significant differences visible throughout all methods.The l 2 -error norms for the benchmark quantities against the finest solution and the corresponding ROC values on each refinement level are given in Table 6.The ROC values of all benchmark quantities for k = 2 show a high-order convergence up to 2.5, whereas the study for k = 3 shows poorer results with ROC below 2. This is again a result regarding the moving interface approach (see Section 7.1), since the temporal discretization with a third order BDF does not match the theoretical requirement for a time integration scheme of at least 2k.Thus, for k = 3 we cannot expect the higher-order convergence rates.However, the higher polynomial order still provides smaller error norms on the same mesh sizes.Note that a reinitialization is performed every 50th time step on the coarsest mesh for k = 3 in order to stabilize the interface evolution.This additional operation may lead to the lower error norm compared to the next finer mesh resulting in a negative ROC value.Note that the error norms l 1 and l ∞ with the corresponding ROC values are found in the Appendix B Tables 7  and 8.The above statements are also valid for both other error norms.
The results of the benchmark quantities on the finest mesh for each study are compared in the following against the finest solutions of the benchmark groups, see Figure 10.In the top row of Figure 10 the temporal evolution of the center of mass y c (t) is plotted.Overall, all numerical solutions agree excellently on the evolution, only a zoom for t = [2.75,3.0] exhibits two distinct trends.Both our numerical solutions coincide very well with FreeLIFE.Likewise, the other two benchmark groups coincide, but are constantly above ours.The terminal rise height at t = 3 for both studies (k = {2, 3}) is given in Appendix A, see Table 9.Besides our results, the reader finds the results of the benchmark groups and of Heimann et al. [5] for a direct comparison.All results are in the range of y c (t = 3) = [1.0799,1.0817].
The temporal evolution of the rise velocity V c (t) depicted in Figure 10 in the middle row shows a very agreement between all groups.A closer look at the maximum velocity at around t = 0.9 displays that our numerical solutions slightly underestimate the peak velocity compared to the benchmark groups.However, one should note that FreeLIFE overestimates both other benchmark groups in the same range, i.e. the deviation is around 4 • 10 −4 .Furthermore, one observes that higher peak values are obtained later in time, see Table 9 in the Appendix.Looking at Heimann et al. [5], their result also tends to underestimate the benchmark groups.
Considering the circularity (Figure 10 in bottom row) the overall qualitative agreement is still very good, but our solutions clearly overestimate the minimum value at around t = 1.9.Zooming in, the results of TP2D and MooNMD indistinguishably coincide.FreeLIFE reaches its minimum value slightly before the other groups.Note that our results are qualitatively closer to the solution of FreeLIFE exhibiting a more gentle slope in the beginning and getting steeper to the end.A reason for the overestimation may be the mass production during the simulation which is, however, comparably small for the finest solution, i.e. the relative mass error is around 10 −3 .The minimum circularity and the corresponding time are given in the Appendix in Table 9.

Test case 2 -Dimpled cap with filaments
For the second test case a mesh study with k = 2 is done on meshes with 1/h = {20, 40, 80}.All simulations are performed with an additional constant mesh refinement on the narrow band with level 1, since the bubble shape exhibits stronger deformations for this case.In Figure 11 the temporal evolution at times t = {1.2,2.1, 3.0} is depicted on the left for the finest solution, where the magnitude of velocity field |u| is plotted.Compared to the first test case (Figure 9) the bubble exhibits a considerably more concave deformation.This evolves into a cap-like shape, where thin filaments starts to emerge from the bubble and get longer and thinner close to the bubble.The terminal shape at t = 3 of the bubble assumes a dimpled cap.A close-up of the right filament is displayed in the middle of Figure 11, where the actual mesh of our simulation is shown.Note that there is a buffer layer of refined cells outside the narrow band, which mark previous cells in the narrow band.The thickness of the filament in gray is determined by the spatial resolution of the adapted mesh.
A comparison between the benchmark groups is given on the right in Figure 11.Since such filaments are strongly mesh dependent, one cannot expect good agreement in such regions.The result of TP2D (1/h = 320) even exhibit a break up of the filament with an additional satellite droplets.Taking a closer look on the transition between the main bubble and the filament, two different shapes may be characterized.The filament is slightly more extended to the outside of the bubble in our solution and the one from MooNMD.On the other side TP2D and FreeLIFE show a more inward curved shape at the transition region.
Comparing the benchmark quantities, the agreement for the center of mass between all methods is again very good, see Figure 14 in the Appendix B. The temporal evolution is not much changed compared to the first test case, but the 0 0.5 1 0 0.5  deviations to the end are larger, where our numerical result lies in between the benchmark groups.The terminal rise height and the other distinct values for the benchmark groups and Heimann et al. [5] are again found in the Appendix B in Table 10.
The temporal evolution of the rise velocity (Figure 12 in top row) exhibits an additional peak around t = 2, which is captured by all methods.The second peak is due to the emerging filaments form the main bubble, see Figure 11.Taking a closer look on both peaks, the deviations at the second one are clearly larger again revealing two distinct trends.Our solution agrees very well with MooNMD underpredicting the whole evolution of rise velocity compared to both other groups.
The qualitative agreement between BoSSS and MooNMD also holds for the circularity, see Figure 12 in bottom row.For the second test case the benchmark groups agree very well up to t ≈ 1.75.After that they show large deviations, especially TP2D which is due to the retractions of the filaments after the break up.All other groups do not exhibit a break up of the filaments.9: Benchmark quantities (test case 1) at distinct values in time for the finest solutions of the respective groups: BoSSS for k = {2, 3}, Heimann et al. [5] for the space pair (X k,2 , M k,1 ) and the benchmark groups TP2D, FreeL-IFE and MooNMD.The values denote: the minimum circularity ¡ c min and the corresponding point in time t| £ c= £ cmin , the maximum rise velocity V c,max and corresponding time t| Vc=Vc,max , and the terminal rise height y c (t = 3).

BoSSS
Heimann et al. [5] benchmark groups [43]   Table 10: Benchmark quantities (test case 1) at distinct values in time for the finest solutions of the respective groups: BoSSS for k = 2, Heimann et al. [5] for the space pair (X k,2 , M k,1 ) and the benchmark groups TP2D, FreeLIFE and MooNMD.The values denote: the minimum circularity ¡ c min and the corresponding point in time t| £ c= £ cmin , both maximum values for the rise velocity V c,max1 , V c,max2 and the corresponding times t| Vc=Vc,max1 , t| Vc=Vc,max2 , and the terminal rise height y c (t = 3).

[
[u]] = 0 on I. (4b) The jump operator [[•]] is defined in (14) and the interface normal n I is pointing from A to B. In the momentum jump condition (4a) on the right-hand side σ denotes the surface tension coefficient and κ the mean curvature of the interface I.At the domain boundary ∂Ω = ∂Ω D ∪ ∂Ω N , describing a disjoint decomposition into Dirichlet ∂Ω D and Neumann regions ∂Ω N , the corresponding boundary conditions read

Figure 1 :
Figure 1: The narrow-band K nb := K cc ∪ K near around the interface includes cut-cells K cc (dark gray cells) and the neighbors K near (light gray) considering the vertices of the cut-cells.The white cells are denoted as far-field cells K far .

Figure 4 :
Figure 4: Numerical setup for the capillary wave test case with an sinusoidal initial perturbation of wavelength λ = 1 and amplitude height a 0 λ.

Figure 8 :
Figure 8: Numerical setup for the rising bubble benchmark test case.

Figure 9 :
Figure 9: Temporal evolution of the rising bubble benchmark test case 1 for times t = {1.2,2.1, 3.0}.The shown results are computed with a polynomial degree of k = 3 on a mesh with 1/h = 60.The plotted field describes the magnitude of the velocity vector field |u|.

Figure 10 :
Figure 10: Comparison to the benchmark groups for the temporal evolution of the y-component of the center of mass y c (top row), of the rise velocity V c (middle row) and the circularity ¡ c (bottom row).In the left column the overall solutions are depicted and the right column provides a zoom for specific time intervals.The plotted solutions show the finest solution of each group for test case 1.

Figure 11 :
Figure 11: Left: Temporal evolution of the rising bubble benchmark test case 2 for times t = {1.2,2.1, 3}.The shown results are computed with a polynomial degree of k = 2 on a mesh with 1/h = 80 and AMR level 1 on the narrow band.The plotted field describes the magnitude of the velocity vector field |u|.Middle: Close-up on the right filament showing the actual numerical mesh (1/h = 80) with AMR level 1. Right: Comparison between the benchmark groups.BoSSS (crosses) with 1/h min = 160, TP2D (squares) with 1/h = 640, FreeLIFE (circles) with 1/h = 160, MooNMD (triangles) with NDOF int = 900.

Figure 12 :
Figure 12: Comparison to the benchmark groups for the temporal evolution of the rise velocity V c (top row) and the circularity ¡ c (bottom row).In the left column the overall solutions are depicted and the right column provides a zoom for specific time intervals.The plotted solutions show the finest solution of each group for test case 2.

Figure 14 :
Figure 14: Comparison to benchmark groups for the temporal evolution of the ycomponent of the centre of mass (test case 2).On the right the overall evolution is depicted and on the left a zoom for t = [2, 3] is provided.The plotted solutions show the finest solution of each group.
Algorithm 5 One time step of the XNSE-solver 1: procedure RunSolverOneStep(t n ,∆t) solve ϕ, u and p for t n+1 = t n + ∆t

Table 1 :
Physical parameters of the Laplace number study for capillary waves

Table 4 :
Physical parameters for both test cases of the rising bubble benchmark.

Table 5 :
Number of background cells (NEL), total number of DOF (NDOF) and number of time steps (NTS) for the rising bubble benchmark test case 1 for polynomial degrees k = {2, 3}.

Table 7 :
l 1 -error norms and ROC for the mesh study of the rising bubble benchmark test case 1.

Table 8 :
l ∞ -error norms and ROC for the mesh study of the rising bubble benchmark test case 1.