The extended discontinuous Galerkin method adapted for moving contact line problems via the generalized Navier boundary condition

In this work, an extended discontinuous Galerkin (extended DG/XDG also called unfitted DG) solver for two‐dimensional flow problems exhibiting moving contact lines is presented. The generalized Navier boundary condition is employed within the XDG discretization for the handling of the moving contact lines. The spatial discretization is based on a symmetric interior penalty method and the numerical treatment of the surface tension force is done via the Laplace–Beltrami formulation. The XDG method adapts the approximation space conformal to the position of the interface and allows a sub‐cell accurate representation within the 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. No adaption of the level‐set evolution algorithm is needed for the extension to moving contact line problems. The developed solver is validated against typical two‐dimensional contact line driven flow phenomena including droplet simulations on a wall and the two‐phase Couette flow.


INTRODUCTION
Considering a droplet on a surface, the interface between the liquid phase and the surrounding phase forms a three-phase contact line L on the surface. Restricted to the two-dimensional case in this work, we deal with single contact points on the surface. Within a continuum mechanical approach the no-slip boundary condition would be the standard choice for the surface. However, describing a contact line that is moving along the boundary, for example, a droplet on a tilted surface, the standard approach is no longer well-posed. Huh and Scriven 1 show that the no-slip boundary condition leads to a diverging stress field, which is non-integrable when approaching the moving contact line. This results further into a non-physical divergence in the energy dissipation rate. In order to regularize the singularity at the moving contact line the introduction of a Navier-slip boundary condition 2 may be used, for example, see Huh and Scriven 1 and Dussan. 3 However, this relaxation of the no-slip boundary condition still features a weak, that is, an integrable singularity of the pressure field. 4 Another feature of slip models is that they exhibit low-velocity regions near the contact line, see Shikhmurzaev. 5 However, such non-physical stagnation regions prevent the development of a rolling motion, which are first observed in Dussan and Davis 6 and further investigated in Dussan 3 and Chen et al. 7 An approach aiming to remove both artificial features is introduced with the interface formation model by Shikhmurzaev. 8,9 This model extends the continuum mechanical description via additional surface mass densities at the interface between all three phases. For each pair, individual mass and momentum balance equations are incorporated and the surface tension also depends on the surface mass distribution. A reduced model that relaxes the impermeability condition and allows mass exchange between the interfaces and the bulk area is presented in Lukyanov and Pryer. 10 Nevertheless, the use of a Navier-slip boundary condition in context of macroscopic flow simulations is widespread, for example, see Bao et al., 11 Gao and Wang, 12 and Shen et al. 13 In this work, we also impose the Navier-slip boundary condition, more precisely the generalized Navier boundary condition, for flow problems with moving contact lines.
Discontinuous Galerkin (DG) methods have gained quite some interest in the context of computational fluid dynamics in the recent years due to their favorable numerical properties. The discretization of a DG method is based on a flux formulation allowing to be locally conservative and the local ansatz functions enable a high-order method even on unstructured meshes. This cell locality is especially favorable for local mesh refinement and results in local mass matrices and sparse operator matrices.
In the context of multi-phase flow simulations with a moving sharp interface, extended methods were developed that adapt the approximation space conformal to the position of the interface in order to adapt costly remeshing. 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. 14 introducing the extended FEM (XFEM) for the simulation of crack growth in solid mechanics.
The first extended method for DG is presented in Bastian and Engwer 15 for the discretization of elliptic scalar problems. In Heimann et al., 16 this approach is applied to incompressible Navier-Stokes two-phase flows. The unfitted DG (UDG) method applies the 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 17 for steady two-phase flows uses a high-order approximation of the interface in combination with a high-order quadrature technique for the implicitly defined domains. The discretization is based on the symmetric interior penalty (SIP) method 18 and the stabilization against small cut-cells is ensured by cell agglomeration. In Smuda and Kummer,19 this method is extended to the instationary case, where the XDG adapted moving interface approach by Kummer et al. 20 is used for the temporal discretization. A signed-distance level-set field is used for the representation of the interface and discretized by a standard DG formulation. The evolution is described by an extension velocity field.
Considering the extension to flow problems including moving contact lines, Reusken et al. 21 present a variational formulation of the GNBC for the stabilized XFEM. The interface is tracked implicitly by the level-set method which does not need any special boundary conditions. An extension of the UDG method is given in Heimann 22 in the context of a conservative level-set method. In this case, the boundary condition at the contact line needs to be adapted for the recompression of the phase-field. Further works in the context of a finite element discretization with moving contact lines are found in Gerbeau and Lelièvre, 23 Manservisi and Scardovelli, 24 Ganesan and Tobiska, 25 and Ganesan. 26 In this work, we adapt the XDG method presented in Smuda and Kummer 19 to allow the direct numerical simulation of moving contact line problems for the two-dimensional case. For this purpose, we employ the generalized Navier boundary condition for the XDG discretization. Therefore, we follow the formulation by Reusken et al. 21 Thus, no additional adaption of the original level-set evolution algorithm is needed. Note that this method was already used for a comparative study of the transient capillary rise, see Gründing et al. 27 The outline of this work is as follows. In Section 2, the continuous problem statement for the transient incompressible two-phase flow with moving contact lines is given. An introduction to the XDG method and the extended discretization for the considered problem is provided in Section 3. The numerical representation of the level-set field and its evolution by the two-staged evolution algorithm is briefly outlined in Section 4. Furthermore, the solver structure is given describing the coupling of the interface evolution with the flow solver. A mesh convergence study is presented for the static case and a validation of the instationary solver is performed by droplet simulation under gravity and the two-phase Couette flow in Section 5. We conclude this work in Section 6.

PROBLEM STATEMENT-TRANSIENT INCOMPRESSIBLE TWO-PHASE FLOW WITH DYNAMIC CONTACT LINES VIA SLIP BOUNDARY CONDITION
We consider the two-phase flow setting within the sharp interface formulation and define the computational domain Ω ⊂ R 2 as the following disjoint partitioning where (t) and (t) denote the time-dependent fluid bulk phases and ℑ(t) the moving interface. Furthermore, we allow especially that the interface is not a closed surface, that is, ℑ ≠ ∅, and we define the contact line at the domain boundary Ω by L(t) ∶= ℑ(t) ∩ Ω = ℑ(t). Since we restrict ourselves to the two-dimensional case, we assume that the sharp interface ℑ(t) ∶= (t) ∩ (t) is a one-dimensional manifold and L(t) accordingly a zero-dimensional manifold. The considered flow problems are described by the transient incompressible Navier-Stokes equations (NSE) denoting the momentum balance equation (2a) and the continuity equation (2b) by where u = u(x, t) is the velocity vector, p = p(x, t) the pressure and g = g(x, t) describes a volume force, for example, gravity. The physical properties density = (x) and dynamic viscosity = (x) are defined piece-wise constant in each phase by The corresponding jump conditions for the two-phase NSE (2) at a material interface ℑ are given by The jump operator [ [⋅] ] is defined in (16) and the interface normal vector n ℑ is pointing from to . The surface tension force on the right-hand side of the momentum jump condition (4a) is given by the surface tension coefficient and the mean curvature of the interface ℑ. The domain boundary Ω is given by the following disjoint decomposition with where Ω D and Ω N describes Dirichlet and Neumann regions. The corresponding boundary conditions read with a given function u D = u D (x, t) on Ω D . On Ω S , we impose the standard Navier-slip boundary condition with u ⋅ n S = 0 and P S (∇u + ∇u T )n S = − S P S u on Ω S .
Here, n S denotes the outer normal on Ω S and P S ∶= I − n S ⊗ n S the orthogonal projection onto Ω S , see Figure 1. The term on the right-hand side of the second condition describes a dissipative friction force linear dependent on u with the phase coefficient of friction S ≥ 0. Note that S introduces a slip length l S with S = ∕l S . For S = 0, the Navier-slip boundary condition reduces to the free-slip boundary condition. At the contact line L, we enforce the following condition F I G U R E 1 Contact line setting on a slip wall Ω S where stat denotes the given static contact angle and the current contact angle. The term on the right-hand side describes a dissipative line force with the contact line velocity u ⋅ n L , where n L is the normal to L and tangential to Ω S . Condition (8) in combination with the Navier-slip boundary condition yields the generalized Navier-slip boundary condition (GNBC). For L = 0, the GNBC describes a quasi-static contact angle model implicitly imposing the static contact angle at all times. In order to close the initial value problem we set an initial condition for the velocity field by where the initial interface position ℑ(0) is given. The material interface evolves according to the bulk velocity u(x, t) at x ∈ ℑ(t).

THE XDG METHOD FOR MOVING CONTACT LINE PROBLEMS
In this section, the introduction of the generalized Navier boundary condition into the extended discontinuous Galerkin (extended DG/XDG) method is presented. First, the XDG space according to Kummer 17 is briefly introduced in Section 3.1 and the corresponding XDG discretization is given in Section 3.2.

The extended discontinuous Galerkin space
Within the XDG method the interface is defined by a sufficiently smooth level-set function (x, t) ∈ C 0 (Ω) that is almost everywhere C 1 (Ω)-continuous. The partitioning of the computational domain Ω is thus given by In this case, the interface normal n ℑ is computed via Assuming that the computational domain Ω is polygonal and simply connected the following basic definitions and operators hold.
Definition 1 (Basic notations). For a polygonal and simply connected computational domain Ω = (t)∪ ℑ(t)∪ (t) ∈ R 2 we define: • the numerical background mesh h = {K 1 , … , K J } that covers the whole domain Ω = ∪ j K j with non-overlapping cells (∫ K j ∩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 ∪ ℑ(t). This set can be subdivided into Γ(t) = Γ int (t) ∪ Γ D ∪ Γ N ∪ Γ S , where Γ int = Γ(t) ⧵ Ω denotes the set of all internal edges, Γ D and Γ N the set of edges imposed with Dirichlet and Neumann boundary conditions, respectively, and Γ S corresponds to the Navier-Slip boundary condition.
• a normal field n Γ on Γ, where it is n Γ = n ℑ on ℑ and represents an outer normal on Ω with n Γ = n Ω .
• for a fixed time t, the broken gradient ∇ h f defines for f ∈ C 1 (Ω ⧵ Γ(t)) the gradient on the domain Ω ⧵ Γ(t). According to that, the broken divergence ∇ h ⋅ f is defined.
Definition 2 (Cut-cell and cut-cell mesh). For a numerical background mesh h = {K 1 , … , K J } we define the time-dependent cut-cells as for a species (t) = { (t), (t)}. The set of all cut-cells defines the time-dependent 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, and K X j, , otherwise the original background cell is recovered for the respective species , that is, K X j, (t) = K j and K X j, (t) = ∅, or K X j, (t) = ∅ and K X j, (t) = K j . 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. Note that the XDG method is essentially a DG method applied on the cut-cell mesh X h (t). Definition 4 (Inner and outer value, jump and average operator). For a field f ∈  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 }} ∶=

Numerical integration on cut-cells
An essential prerequisite of the XDG method is a high-order accurate numerical integration of cut-cell integrals such as which represent integrals over cut-cell volumes and surfaces and integrals along the interface ℑ. In this work, we use the quadrature method proposed by Saye. 28

Cell agglomeration
Obviously, the discretization by the XDG method may introduce arbitrarily small cut-cells, where the volume of the cut-cell K X j, is small compared to the background cell K j , since the interface position with respect to the background mesh is also completely arbitrary. In context of the SIP discretization for the viscous terms (A5) such cut-cells give rise to large penalty parameters (A7) and undesirably high condition numbers which consequently leads to stability issues. This can be handled by agglomerating small cut-cells to their largest neighbor. 17 Furthermore, the cell agglomeration is used in the context of the temporal discretization to ensure that all involved cut-cell meshes X h (t n−2 ), … , X h (t n+1 ) share the same mesh topology. 20 Precisely, some cut-cell K X j, is marked for agglomeration if one or more of the following conditions are met: • The volume fraction of the cut-cell is below a certain threshold, that is, |K X j, (t n+1 )|∕|K j | ≤ . In this work, the agglomeration threshold is set to = 0.1 for all presented simulations.
As a target for agglomeration, the largest neighboring cell within the same species, which is not itself marked for agglomeration, is selected. If no such cell is available, that is, if all neighboring cells should also be agglomerated, a chain of agglomerations must be formed. Note that agglomeration chains or pairs always stay within the same species, that is, there is no agglomeration across the interface. For the results presented in this work, the occurrence of agglomeration chains is a rare event and the chains are usually short. For more complex cases, for example, droplet collisions in three dimensions, one would expect much larger chains. The suitability of the agglomeration for these cases will be the issue of future works.
Formally, the agglomeration is a projection from the space P X k ( h , t) onto the agglomerated sub-space: be the agglomerated mesh, which is obtained from the cut-cell mesh X h (t) by following the agglomeration procedure described above. Then, the agglomerated XDG space is defined as From an implementation perspective, we do not modify the mesh in the agglomeration procedure. Instead, the assembly of operator matrices and resp. operator evaluations are first performed on the non-agglomerated cut-cell mesh. Then, matrices and right-hand-side vectors are projected onto the agglomerated mesh, in a fashion that is very similar to a geometric multigrid method. For details on the mesh agglomeration algebra, we refer to the work of Kummer et al. 29

XDG discretization of the incompressible two-phase NSE with GNBC
The function space of the ansatz and test functions for the velocity field u ∈ R 2 and the pressure field p ∈ R is defined by where k describes a degree vector with k = {k, k, k ′ } = k and = 1, 2, 3. In this work, the velocity components are discretized by an agglomerated XDG space of order k and the pressure field of order k ′ = k − 1.
Following the XDG discretization in Smuda and Kummer, 19 we propose the following discretization for the moving contact line problem (2), (4) and boundary conditions (6) with the GNBC (7) and (8): The first line including the discretizations of the temporal term m(−, … ), the convective term c(−, −, −), and the pressure gradient and the continuity term, both denoted by b(−, −), are unchanged from the work of Smuda and Kummer. 19 Furthermore, the terms a(−, −) and s(−, −) describing the discretizations of the viscous terms and the summarized source terms are unchanged. Since this work focuses on the introduction of the GNBC, those terms are given in Appendix A. In the remainder of this section, the integration domains of the remaining terms are set to the time-level t n+1 , that is, the interface is fixed at ℑ = ℑ(t n+1 ).
Introducing the GNBC into the XDG discretization, we follow the variational formulation of Reusken et al. 21 Considering the Navier slip boundary condition (7) at Ω S we write v = (v ⋅ n S )n S + P S v. This way the normal components enforce the no-penetration condition, that is, u ⋅ n S , and for the tangential parts we introduce the friction force − S P S u. Thus, the resulting boundary form a S (−, −) in context of the extended SIP formulation (see (A5)) reads The corresponding penalty parameter is given in Appendix A. Enforcing the localized condition (8) at the contact line L, we use the same approach by using v = (v ⋅ n S )n S + P S v for the boundary term ∫ L L ⋅ v dl, which reads Noting that P S L = cos( )n L , we introduce the static contact angle stat and the dissipative line force L (u ⋅ n L )n L via (8). Since the line force is linear dependent on u, this term is added to the left-hand side of (21) with The prescribed static contact angle stat and the remaining normal component are added to the right-hand side with

INTERFACE EVOLUTION AND COUPLING WITH NSE
In this section, the level-set representation is introduced and the evolution described (Section 4.1). Furthermore, the solver structure for coupling the level-set evolution with the NSE is presented (Section 4.2). In context of the XDG method for transient two-phase flow problems, 19 two level-set fields are introduced The first one is used for handling the level-set evolution and is discretized by a standard DG space with the same polynomial degree k as the velocity components in u. The second one is used for the discretization, since the XDG method requires at least that (x, t) ∈ C 0 (Ω). Therefore, it is given as the L 2 -projection of evo constrained by continuity conditions at the cell boundaries, see Equation (31). The polynomial degree is chosen to be k + 1.

Level-set evolution with contact lines
The level-set function itself is defined by its signed-distance formulation leading to the property that |∇ | = 1. The evolution equation for the level-set field evo is given by where the evolution velocity is given by an extension velocity u ext field constructed from the bulk velocity u. In Smuda and Kummer, 19 a two-staged marching algorithm is presented, where a high-order elliptic PDE-based approach on the F I G U R E 2 The extension problem (28) is not defined in regions where ∇ ⋅ n Ω > 0 (gray area). Expanding the velocity vector at the contact line leads to a vanishing-viscosity solution interface cells and a low-order fast-marching on the neighboring cells is performed. Note that the level-set fields only have valid signed-distance values on a subset around the interface denoted as the narrow-band nb . This narrow-band nb ∶= cc ∪ near includes the set of cut-cells cc and its neighbors near that share at least one point. On the far field cells far , the value is either set to evo = −1 for cells in or evo = +1 in .
The entire interface evolution algorithm within the XDG method is essentially subdivided into four steps: • Step 1: Construct the extension velocity field u ext by (28) via a two-staged marching algorithm on the narrow-band. • Step 2: Advect evo via (30) using a Runge-Kutta scheme. • Step 3: Project evo onto the continuous level-set field by the constrained L 2 -projection (31).

•
Step 4: Update narrow band nb and integration regions.
In the following, a few aspects of the above evolution algorithm are presented in more detail.

4.1.1
Step 1: Construction of the extension velocity field u ext Considering the advection of the interface, the only valid velocity for the evolution should be the one at the interface itself. Thus, we construct the evolution velocity in (27) 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}: Note that in context of the XDG method the approximation space for the bulk velocity is in P X k ( h ), whereas the level-set evolution needs a velocity field u ext ∈ P k ( h ) according to evo . Therefore, we introduce a density-averaged velocity directly at the interface (28a). The rationale behind this choice is that, for example, in the (limit) case of free-surface flows, the interface motion should be mainly determined by the heavier fluid. Obviously, one could think of other choices, which we do not investigate here.
In the context of moving contact line problems, the interface cuts the computational domain Ω at L. Thus, there is a region where the extension problem is not defined, that is, the gray area with ∇ ⋅ n Ω > 0 in Figure 2. In our implementation, this case is handled via expanding the velocity vector at the contact line like an expansion fan, see Figure 2. Thus, the extension problem (28a) in the gray area is not equals zero, similar to a vanishing viscosity solution.
The evolution algorithm presented in Smuda and Kummer 19 is divided into two stages. In the first stage, the extension velocity is computed monolithic with high-order accuracy on the cut-cells cc using the reformulation of the EVP (28) into an elliptic PDE. 30 In the second stage, the extension velocity on the near field near is constructed geometrically in each cell separately using an adapted fast-marching 31,32 procedure. The solution is then constructed node-wise from the accepted edges into the cell K l and projected on P k ({K l }).

Level-set reinitialization
The reinitialization problem is described by the following Eikonal equation During the level-set evolution new cells entering the narrow-band need to be reinitialized. Therefore, a fast-marching 31,32 procedure is utilized, where the local solver for one cell is based on the elliptic reinitialization method by Basting and Kuzmin 33 and adapted for DG in Utz et al. 34 A low-order initial solution for evo is obtained node-wise by a direct geometric approach. Furthermore, an optional full reinitialization of the level-set field evo on the narrow-band may be performed at predefined time periods. In this case, the elliptic reinitialization 34 is solved monolithic on the cut-cells and the fast-marching procedure is executed on the near field.

4.1.3
Step 2: Advection of evo The new level-set field n+1 evo is computed according to the evolution equation (27) by where Γ int,nb denotes the internal edges in 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. 35 At last, the advected level-set field n+1 evo is additionally stabilized by penalizing jumps at inner edges, that is, [ [ ] ] 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.

4.1.4
Step 3: Enforcing the interface continuity For the construction of the continuous level-set field , 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 nb and the continuity conditions are defined on the set of internal edges in the narrow-band Γ int,nb .

Solver structure-Coupling flow solver and interface evolution
In this work, an explicit coupling between the interface evolution and the NSE is employed. In this case, the interface is only updated once at the beginning. Note that in this case the surface tension force is also treated explicitly. An initial solution on the new domain Agg( X h (t n+1 )) is given by extrapolation of the old solutions u n and p n . The extrapolation operator is defined as follows: where X, (t) is a basis of the space P X, k ( h , t). Since the spaces for t n and t n+1 have the same mesh topology, due to the agglomeration of vanishing and appearing cut-cells, there is no ambiguity: the interface always stays within the same agglomerated cells, therefore any extrapolation during movement of the interface is uniquely defined by the local polynomial representation. In (32), it is assumed that X, (t) is inferred from a basis on the background mesh which does not depend on time. Then the DG coefficientsf n do not change, and thus the implementation is trivial. 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. Within each iteration the linearized system (linearized via u * in the convective terms, see (A3)) of the discretized NSE (21) is updated. In order to reduce the condition number of the system, the cell agglomeration is performed and a block Jacobi preconditioning is used, see Kummer 17 for details. In this work, the preconditioned linear systems are solved using the direct solver MUMPS, 36,37 since the number of degrees of freedom throughout all presented simulations (Section 5) is small enough to be solved most effectively by a direct solver.

4.2.1
Adaptive mesh refinement The presented method allows the use of local adaptive mesh refinement (AMR) during the simulation. In this case, a background cell is divided into four equal sized sub-cells, see Figure 3. For adaptions on higher levels, it is always ensured that neighboring cells exhibit a 2:1 cell ratio on the edges. Thus, we counteract undesired locking-effects on the refined cells by much larger cells. Furthermore, we do not coarsen on cut-cells. 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. In the case of moving contact line problems, we additionally refine on boundaries that are imposed by the GNBC, see Figure 3 on the right.

NUMERICAL RESULTS
All

Droplet on slip wall
First, static simulations of a droplet set on a slip wall are investigated. In this case, the droplet assumes a fixed non-equilibrium state. Therefore, we neglect the temporal term and switch off the level-set evolution. The fixed non-equilibrium state of the droplet is given by

Condition number study
The condition number of the operator matrix cond(Op), that is, in this case the Stokes-operator, is a measure for the numerical solvability of a linear system with Op(x) = b. For small condition numbers, relative small errors in the right-hand side b only lead to relative small errors in the solution x. This property is essential for the stability and accuracy of a numerical method. In context of iterative solver methods, a small condition number ensures that the reduction of the residual in each iteration also results in a convergent solution. However, absolute condition number values are not sufficient to evaluate the solvability and stability of the discretized linear system. Thus, a mesh study is conducted for the static droplet test case described above. The study with polynomial degree k = 2 is performed on meshes with 8 × 16, 16 × 32, 32 × 64, 64 × 128, and 96 × 192 cells. The condition numbers are computed with MATLAB via condest(A) which computes a lower bound for the 1-norm condition number of a square matrix A. The results for the four cases with stat = {90 • , 120 • } and S = {5, 0} are shown in Figure 4. The corresponding lines exhibit slopes of 1.99 for S = 0 and of 1.90 for S = 5 w.r.t. the grid size h. This means that the total condition number does not increase more than the size of the corresponding discretized operator matrix. For the two-dimensional case, we expect a maximum slope of 2, since the increase of the grid resolution in each dimension by a factor of 2 leads to a 4 times larger operator matrix. Exhibiting a slightly smaller value, the presented discretization does not introduce grid-dependent numerical artifacts and behaves as expected under grid refinement. One should observe that the results are independent of the prescribed static contact stat = {90 • , 120 • }. The absolute values of the total condition numbers for S = 0 describing a free-slip boundary condition are larger compared to S = 5.

5.1.2
Convergence study Looking at the EOC of both cases with stat = 120 • one can see that for all tested polynomial degrees k = {2, 3, 4} the convergence order degenerates to EOC ≈ 1 for the velocity norm and, respectively, to EOC ≈ 0.1 for the pressure norm. This is a result of the weak pressure singularity with p ∼ ln(r) near the contact line, where r denotes the distance to the contact line L (see Sprittles and Shikhmurzaev 4 ). Such a solution cannot be approximated accurately by polynomial ansatz functions as used in our discretization. The representation of the pressure singularity on the finest mesh is shown in Figure 5 top row for both cases S = {5, 0} with stat = 120 • . The singularity is positive in the droplet phase and negative in the surrounding phase. Note that the pressure singularity inside the droplet is broader for the free-slip case.
Considering the cases for stat = 90 • , where the current contact angle assumes the static equilibrium contact angle, one observes no singularity ( Figure 5, bottom row). However, for S = 5 one notices a drop of the pressure level at the contact line. The results of the pressure fields are directly reflected in the corresponding EOC (Table 1). For S = 5, the  convergence order for both the velocity and pressure norm increases by an order of 1, but still the expected order cannot be recovered. However, looking at the free-slip case we almost regain the expected orders of h k+1 . This is in accordance with the findings in Fricke et al., 38 where it is stated that a regular, non-trivial, quasi-stationary solution only exists in the free-slip case.

Spurious velocities-Transient simulations
In this part, we consider a droplet in its equilibrium state, that is, a circular droplet with a contact angle of = stat within a zero-velocity field. According to the Young-Laplace equation, the corresponding pressure inside the droplet (phase ) with radius r is given by where 1∕r describes the curvature in the two-dimensional case. In order to quantify the discretization error due to numerical inaccuracies regarding the surface tension force, for example, interface representation, we present the following transient setup similar for the static simulations. The droplet describes a circular shape with radius r = 0. 25 Table 2, 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 Ω ⧵ ℑ on the terminal time step at t = 125. The error norms seem to converge already on the second mesh, what may result from the long simulation time and fixed time step size for all simulations. In Figure 6, on the left the temporal evolution of the spurious velocities is depicted. All numerical solutions exhibit a large deviation from the zero-velocity field in the beginning, which gets smaller for finer meshes. However, one observes that there is a slight growing tendency to the end, which may be due to the fixed time step size. Further, taking the change rate of the interface area, that is, interface length in R 2 , into account ( Figure 6 on the right) it shows that for the coarser meshes the interface is slightly increasing, whereas it decreases for the finest mesh.

Droplet spreading under gravity
For a second test case, we consider a semi-circular droplet on the wall with an initial contact angle of 0 = 90 • which does not describe an equilibrium state, since the imposed static contact angle is stat < 90 • . Thus, the droplet starts to spread in order to regain its equilibrium shape. However, the terminal shape varies under the influence of an additional gravity force g = −ge y . For a study with different gravity forces g, we follow the setting proposed in Reference 39. The computational domain is given by Ω = [−3r 0 , 3r 0 ] × [0, 2r 0 ], where the initial droplet radius is set to r 0 = 0.01 m. The physical parameters are given in Table 3. Imposing a static contact angle of stat = 50 • , the shape of the droplet is only characterized by the Eötvös number Eo which is defined in this test case as where L ≪ G denotes the density of the liquid droplet and G the density of the surrounding gaseous phase. For Eo ≪ 1, the terminal shape is dominated by surface tension effects and the droplet assumes a circular cap with the given static contact angle stat . Thus, the droplet thickness, that is, highest point of the droplet shape, is given by For Eo ≫ 1, gravity forces dominates the shape and the droplet spreads to a puddle with a thickness of which is directly proportional to the capillary length l cap = √ ∕ L g. In our study, we consider the following range of Eötvös numbers Eo = {0.01, 0.1, 1, 2, 5}. The corresponding gravity values are given in Table 4.
The study was performed on a 30 × 10-mesh with AMR level 1 at the interface and the slip boundaries. The induced slip length is set to l S = h∕100. The simulations are done with a polynomial degree of k = 2 and run until t = 1 with Δt = 10 −4 . In Figure 7, the normalized thickness e * = e∕e 0 of our results are displayed and compared to the results of Dupont and Legendre 39 and Gründing. 40 The latter presents an exact solution for interface shape of the considered inverse problem. In addition to the exact solution, both asymptotic solutions for Eo ≪ 1, that is, e * = 1, and Eo ≫ 1 (Equation 37) are plotted. All results show very good agreement to the exact solution for both asymptotic regimes and in the transition region. The results of BoSSS and Gründing match excellently, only for Eo = 5 the result of Dupont and Legendre slightly underestimates (relative deviation ≈ 4%) the other results.
The terminal shapes for Eo = {0.01, 2, 5} of our simulations are shown in Figure 8. One can clearly see the circular shape for the low Eötvös region (Eo = 0.01) and the deformation to a puddle for higher Eo. In the transition, one notices the receding surface tension forces and increasing effect of gravity. All shapes exhibit the imposed static contact angle of stat = 50 • . The maximum mass loss for the terminal time step in the presented study is around 0.3% for Eo = 5.

Temporal evolution
Considering the cases for Eo = {0.01, 2, 5}, we take a closer look on the temporal evolution of the spreading process. Therefore, we performed the simulations on a finer mesh with AMR level 2, which refines additionally inside the droplet. We note that the refined simulations exhibit a more dynamic behavior, so the simulation for Eo = 5 is performed on a larger domain with Ω = [−4r 0 , 4r 0 ] × [0, 2r 0 ] in order to prevent the contact line to touch the outer boundaries. However, the terminal thickness e * and shapes remain unchanged (compare mesh study for the capillary rise in Gründing et al. 27 ). In Figure 9, the temporal evolution of the normalized droplet height at x = 0 with h * = h(x = 0)∕e 0 is depicted on the left. On the right, the normalized wetting length l * = l∕l 0 is given, where l 0 denotes the wetting length for a circular cap with One observes for both properties the stronger deforming effect due to the increasing gravitational force, that is, increasing Eötvös number. However, one should note that for smaller Eo the dynamic behavior exhibits some higher oscillations modes, especially for Eo = 0.01. Looking at the droplet height h * (Figure 9, left) an additional peak around t = 0.045 is visible. This results from a surface wave starting at the moving contact line due to the adjustment of the imposed static contact line stat = 50 • , see Figure 10. Since the gravitational force and its acceleration on the droplet is small compared to the travelling wave, the interface at the top moves upwards around t = 0.045. In Figure 11, contour plots of the magnitude velocity and pressure field for Eo = 2 are presented at times t = {0.05, 0.1, 0.4}. The contour plots for Eo = 0.01 (t = {0.05, 0.1, 0.3}) and Eo = 5 (t = {0.05, 0.1, 0.5}) are found in Appendix C, see Figure C1 and Figure C2. Looking at the velocity fields, one sees a low velocity region at the droplet center in all cases. Furthermore, there is a thin shear layer at the slip boundary, which gets thinner for higher Eötvös numbers. The high velocity regions are found above those shear layers and directly at the moving contact line. Looking at the corresponding pressure fields, one observes some irregularities in the contour lines. Reason for that is the pressure singularity at the moving contact line due to the GNBC, which affects the pressure field especially in the beginning of the simulations. Note that the irregularities reduce to a small region around the contact line at the last depicted pressure plot for each case, that is, t = 0.3 for Eo = 0.01, t = 0.4 for Eo = 2, and t = 0.5 for Eo = 5. The pressure singularity vanishes for the static equilibrium case. Note that the maximum mass loss for the terminal time step of the refined study reduced to around 0.04% (Eo = 5).

Two-phase Couette flow
For this test case, we follow the MD setup of Qian et al. 41 for a two-phase Couette flow, which is also presented in Gerbeau and Lelièvre. 23 The initial setup is depicted in Figure 12 Table 5. The upper wall is imposed with a velocity u = u wall e x in positive x-direction and in negative direction with u = −u wall e x on the lower wall. The test cases are characterized by its quasi-stationary flow pattern, where the first test case describes a symmetric (Figure 13, bottom) and the second an asymmetric one (Figure 15, bottom). Both test cases are computed up to t = 160. In the following, we present for each test case a simulation study of comparable number of degrees of freedom (NDOF) with different pairs of polynomial orders and meshes. The setup descriptions are given in Table 6. For all setups a time step size of Δt = 0.02 is chosen.

F I G U R E 13
In Figure 13, the temporal evolution of the magnitude of the velocity field is displayed. One can clearly observe the symmetric evolution of the flow field at both interfaces. Furthermore, the color scale indicates a quasi-static state, where the interfaces assume a static position due to the slip boundary condition. Between the interfaces two circulating flows emerge within each phase. One should note that the velocity at the moving walls do not resume the prescribed wall velocity, that is, u wall = 0.25 for the symmetric case. Looking in Figure 14 Figure 14 at the top left plot. Here, the results of the first setup assume a slightly larger extend of movement than the other two.
Considering the asymmetric test case, the temporal evolution is depicted in Figure 15. Here, the static contact angle is changed to stat = arccos(0.38) ≈ 67.67 • and one can clearly observe the enforcement of the quasi-static contact line model ( L = 0) in the beginning. Thus, the upper contact line of the left interface moves against the imposed direction of the upper wall and the lower one is moving faster than the wall (u wall = 0.2). The corresponding evolution of the contact line position and the velocity are given in Figure 14 in the bottom row. Note that the contact lines more quickly reach their static equilibrium state. Furthermore, the circulating flow in the middle area (phase ) exhibit a slower velocity on the moving wall due to the larger slip length induced by < . In Figure 16, the terminal shapes at t = 160 are compared. In addition, the results of Gerbeau and Lelièvre 23 (circles) and Qian et al. 41

CONCLUSION
Within this work, the XDG method for two-dimensional instationary two-phase flows is extended to flow problems exhibiting moving contact lines at the domain boundaries. For this purpose, the generalized Navier boundary condition is adapted for the XDG discretization. Here, the Laplace-Beltrami formulation of the surface tension force allows the enforcement of the prescribed static contact angle. The existing level-set evolution algorithm does not need to be adapted. The condition number study for a fixed droplet on a slip wall ensures a correct and stable implementation of the GNBC in the solver framework. This also holds for the free-slip boundary condition. The corresponding convergence study showed two characteristics. First, the high-order spatial convergence rate degenerates due to the inherent weak pressure singularity near the contact line. Second, the high-order convergence is recovered for the quasi-stationary case if a free-slip boundary condition is imposed. Regarding spurious velocities transient simulations of a droplet in equilibrium state were performed. The simulations showed stable results under mesh refinement, but evaluating the change rate of the interface length indicates some stability issues for longer simulation times.
The results for a spreading droplet under the influence of gravity forces showed excellent agreement to the exact solution for the terminal droplet height. Both regimes from a circular shape dominated by surface tension effects to a puddle dominated by gravity forces were presented. Furthermore, the temporal evolution of the droplet height and wetting length were investigated and contour plots of the velocity and pressure fields were presented.
At last a study with comparable NDOF was performed for the two-phase Couette flow test case. The results for the different pairs of polynomial order and meshes showed very good agreement among each other. This holds for both the symmetric and the asymmetric case.

DATA AVAILABILITY STATEMENT
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-330.

ORCID
Martin Smuda https://orcid.org/0000-0001-7931-5176 For the convective term c(u * , u, v), a local Lax-Friedrichs flux with the stabilization parameter , see Kummer,17 is employed where u * describes a linearization velocity for the non-linear convective term. Note that the interface integral is excluded due to the moving interface approach. The discretization of the pressure gradient and the continuity term, both denoted by b(−, −), is given by An extension to the standard SIP method 18 including the transposed term ∇u T , see Kummer,17 is used for the discretization of the viscous terms with The penalty parameter in the viscous terms (A5) is chosen according to where the local penalty parameter̃is computed bỹ= for an agglomerated cell K X , that is, K X ∈ Agg( X h (t)) The safety factor 0 is set to 0 = 4, see Kummer 17 for details. The geometric factor | K X |∕|K X | is directly provided by the construction of the cut-cell quadrature rules, see Equation (18).
We specify the term s(v, q) on the right-hand side of the variational formulation (A1) as follows The two terms in the top line describe the Dirichlet boundary conditions for the convective part (A3) and the viscous parts (A5). The first term in the bottom line corresponds to the Dirichlet boundary condition for the continuity equation (A4) and the second term denotes the discretization of the volume term g. The last two terms correspond to the Laplace-Beltrami formulation for the discretization of the surface tension force. Here, P ℑ ∶= I − n ℑ ⊗ n ℑ describes the projection tensor onto the interface ℑ and the interface gradient operator ∇ ℑ is defined as ∇ ℑ ∶= P ℑ ∇. The vector denotes the tangential to ℑ and the corresponding integral term denotes point-measures in R 2 . Note that the surface tension force may also be discretized directly in form of the momentum jump condition (4a). Thus, the curvature-based discretization is given by ∮ ℑ n ℑ ⋅ v dS. Employing the Stokes integral theorem, one can infer for a single cell K the equality This can be applied to convert the curvature-based discretization of the surface tension into the form found in Equation (A8). In our simulations, we found this formulation to be less critical with respect to the stability of the integrated time evolution of level-set and flow solver, although we cannot present systematic studies on this issue yet. It seems that curvature-based discretizations require some filtering of the curvature field. A patch-recovery filtering technique for this purpose was presented by Kummer and Warburton. 42 Unfortunately, it is not known how much filtering is sufficient. On the other hand, numerical tests indicate also that too much filtering may lead to instabilities and/or non-physical results.