Discontinuous Galerkin method for incompressible two-phase flows
Janick Thomas Gerstenberger Samuel Burbulla Dietmar Kröner
math.NA cs.NA physics.flu-dyn 76T10, 76D45, 76D05, 65M60


Discontinuous Galerkin method for incompressible two-phase flows


  Janick Gerstenberger AAM, Albert-Ludwigs-Universität Freiburg,
  Hermann-Herder-Str. 10, 79104 Freiburg, Germany,

  janick.gerstenberger@mathematik.uni-freiburg.de
  Samuel Burbulla IANS, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany,

  samuel.burbulla@mathematik.uni-stuttgart.de
  Dietmar Kröner AAM, Albert-Ludwigs-Universität Freiburg,
  Hermann-Herder-Str. 10, 79104 Freiburg, Germany,

  dietmar@mathematik.uni-freiburg.de

*
    Janick Gerstenberger, Samuel Burbulla and Dietmar Kröner
    Final Version, December 30, 2023
============================================================




  In this contribution we present a local discontinuous Galerkin (LDG) pressure-correction scheme for the incompressible Navier-Stokes equations.
  The scheme does not need penalty parameters and satisfies the discrete continuity equation exactly.
  The scheme is especially suitable for two-phase flow when used with a piecewise-linear interface construction (PLIC) volume-of-fluid (VoF) method and cut-cell quadratures.
  




§ INTRODUCTION


Sharp interface models for incompressible two-phase flows have gained in popularity in recent years.
These models combine incompressible flows in the bulk domains with jump conditions along the interfaces that separate the fluids,
which model fluid interactions and surface effects, like surface tension.

Discontinuous Galerkin methods are a popular choice for solving incompressible flow problems due to their local mass conservation property and potentially high order of convergence.

Designing schemes for solving such two-phase flows presents several challenges.
The choice of method for the phase transport has implications for the interface representation and conservation properties.
Volume-of-fluid <cit.> methods are conservative but their interface representations is in general non-continuous,
while level-set methods have continuous interface representations, but are not conservative by themselves.
There are several further methods that combine features from or generalize volume-of-fluid and level-set methods, but these are comparatively more complex.
Incompressible flow problems have a saddle-point structure, which can make them computationally difficult/expensive to solve.
Splitting methods like the various projection methods where introduced to decouple the problems into simpler ones.
Projections methods replace the saddle-point problem with a advection-diffusion equation for the velocity and Poisson problem for the pressure,
which are computationally simpler to solve.
But for most discontinuous Galerkin discretizations the discrete continuity equation is not satisfied without some postprocessing techniques like the H(div) reconstruction presented in <cit.>.
Stability in regards to high coefficient ratios and strong surface effects are further issues that need to be addressed.

We present a discontinuous Galerkin pressure-correction method for incompressible two-phase flow that is robust in regards to coefficient ratios.
The scheme is simple to implement and has shown good results for benchmarks problems and some numerical experiments with realistic data.

This contribution is structured as follows: We first briefly present the model for incompressible two-phase flow without phase transitions.
Next we give some notation, present the modified LDG method without penalization introduced in <cit.>
and present the new LDG pressure-correction scheme for incompressible two-phase flows based on it.
Then we present numerical experiments for benchmark problems from the literature <cit.>.
Lastly we give some concluding remarks and an outlook on further work.




§ MODEL

The model we consider is a simplification of the sharp interface model presented in <cit.> without phase transitions.
Let Ω⊂ℝ^n, n=2, 3 be a bounded domain, which is divided into two disjunct phases: at time t, phase i = 1, 2 occupies subdomain Ω_i(t) ⊂Ω.
The boundary between the phase is the (phase) interface Γ(t) ∂Ω_1∩∂Ω_2.

Let u⃗ denote the velocity field, p the pressure field, ν the outer normal on Ω_1, κ the mean curvature of Γ(t) and q q_2 -q_1 the jump of the variable q across Γ(t).
In addition the constants ϱ_1, ϱ_2 > 0 denote the densities, μ_1, μ_2 > 0 the viscosities of the phases, g⃗ the gravitational acceleration and σ the coefficient of surface tension.
In the following we drop the indices on density and viscosity and keep in mind that these coefficients depend on the phase.

Then the model is given by the incompressible Navier-Stokes equations in each phase,

    .
    ∂_t (ϱu⃗) + ∇·(ϱu⃗⊗u⃗ + T)   = ϱg⃗
    ∇·u⃗   = 0
      }  in Ω_1, Ω_2 ,

with the stress tensor T = pI - 2 𝒟u⃗, 𝒟u⃗ = S(∇u⃗)  (∇u⃗ + ∇u⃗^⊤) / 2,
augmented with the following jump conditions
  
    .
    Tν⃗   = -σκν⃗
    u⃗   = 0
      }  on Γ.

Additionally we impose either no-slip or free-slip conditions on the boundary of Ω.




§ DISCRETIZATION


To discretize the model above we employ a primal LDG method <cit.>.
It can be formulated in terms of a discrete gradient operator that is composed of the elementwise gradient and a lifting of the jumps into the piecewise discrete space.
By constructing the liftings one order higher than the used discrete space the method is rendered stable without penalty parameters <cit.>.
From these building blocks we construct an incremental pressure-projection scheme that satisfies the discrete continuity equation.

Since we require strict mass-conservation we choose to use a PLIC-VoF method <cit.>.
These methods have the disadvantage that the reconstructed interface is in general not continuous and its curvature needs further approximations <cit.>, but is mass conservative.

We chose this LDG method because with the discrete gradients and no penalty terms  no (explicit) evaluations of integrals over element boundaries are needed.
In our experience the non-continuous discrete interface reconstructions on the inter-element boundaries lead to issues around the interface.




 §.§ Notation and Liftings

  To derive the discretization we first introduce some notation.
  Let 𝒯_h be a triangulation of Ω into elements E.
  By Σ_h^I we denote the set of all interior intersections e of elements E^-, E^+ ∈𝒯_h with e = E^- ∩ E^+ ≠∅,
  by Σ^D_h the set of all intersections with Dirichlet boundary values,
  by Σ^N_h the set of intersections with Neumann boundary values
  and by Σ_h = Σ_h^I ∪Σ_h^D ∪Σ_h^N the set of all intersections.

  For e ∈Σ_h^I, with e = E^- ∩ E^+,  we introduce the jump and the average
  
    
    u   = u_ |_E^- - u_ |_E^+ ,
    u   = 12(u_ |_E^- + u_ |_E^+)

  

  We define the discrete spaces of piecewise polynomials of degree ≤ k
  
    𝒱_k^d = {v⃗∈[L^2(Ω)]^d|v⃗_ |_E∈[𝒫_k(E)]^d, E∈𝒯_h}

  and 𝒱_k = 𝒱_k^1 in the scalar case.
  The piecewise L^2 inner product over 𝒯_h is given by
  
    v⃗w⃗ = ∑_E∈𝒯_hv⃗w⃗_E = ∑_E∈𝒯_h∫_Ev⃗·w⃗


  The gradient lifting operators R: 𝒱_k ↦𝒱_k+1^d and R_a: 𝒱_k ↦𝒱_k+1^d
   are defined by
  
    
    R(v)w⃗
          = ∑_e ∈Σ_h^I∫_evw⃗·n⃗_e ,

    
    R_a(v)w⃗
          = R(v)w⃗
          + ∑_e ∈Σ_h^D∫_e(v-a) w⃗·n⃗_e ,

  
  for all w⃗∈𝒱_k+1^d.
  Similarly, the divergence lifting operators M: 𝒱_k+1^d ↦𝒱_k, M_b⃗: 𝒱_k+1^d ↦𝒱_k
  are defined by
  
    
    M(v⃗)w
          = ∑_e ∈Σ_h^I∫_ev⃗·n⃗_e w ,

    
    M_b⃗(v⃗)w
          = M(v⃗)w
          + ∑_e ∈Σ_h^D∫_e(v⃗-b⃗)·n⃗_e w ,

  
  for all w∈𝒱_k. Here a, b⃗ are the respective Dirichlet boundary conditions.

  With the liftings we can now define the lifted DG gradient and divergence
  
    
    ∇v    = ∇_h v - R(v),    ∇_ gv    = ∇_h v - R_g(v),     v    ∈𝒱_k,
    ∇·v⃗   = ∇_ h·v⃗ - M(v⃗),    ∇_g·v⃗   = ∇_ h·v⃗ - M_g(v⃗),    v⃗   ∈𝒱_k^d.

  
  The lifted derivatives with homogeneous Dirichlet boundary (a =0, b⃗ = 0) satisfy the following discrete integration-by-parts identities, as shown in <cit.>:
  
    
    ∇_0·v⃗w   = -v⃗∇w
    ∇·v⃗w   = -v⃗∇_0w

  
  for  all v⃗∈𝒱_k+1^d, w∈𝒱_k.
  This means the lifted derivatives are adjoint to each other.
  These identies are useful for defining of projection methods with respect to the lifted derivatives.



 §.§ Unpenalized LDG Scheme


  As a simple example we now consider the discretization of the Poisson equation -Δu = f with homogeneous Dirichlet boundary data.
  The modified LDG method <cit.> reads: Find u∈𝒱_k such that
  
    ∇_0u∇_0v = fv    ∀ v∈𝒱_k.

  This scheme with order k+1 liftings is stable without adding penalty terms and can also be written in a “strong”-form by using the integration-by-parts identity (<ref>)
  
    -Δ_1.5mu 0 uv -∇·∇_0uv = fv    ∀ v∈𝒱_k.





 §.§ Two-Phase LDG Scheme


We now present our primal LDG pressure-correction scheme for two-phase flow.

The scheme first split into an explicit (linearized) advection step for the momentum/velocity and the phase transport and an implicit Stokes step.
Because the momentum and the phase interface are transported at the same rate, we can formulate the explicit step in the velocity,
which also means this part of the scheme does not depend on the phase interface.
The Stokes step is then further split into an implicit momentum step, a pressure Poisson step and an update step.

To sharply resolve the phases the terms containing a phase dependent coefficient are integrated using cut-cell quadratures.
The cut-cell quadratures are constructed by cutting the elements containing an interface reconstruction with its interface reconstruction, subtriangulating the part of each phase and using standard simplex quadratures in the the subtriangulations.
Surface tension effect are included by integration of the jump condition over segments of the interface reconstruction.
Using the lifted derivatives we eliminate (explicit) evaluations of integrals containing phase dependent coefficients along element boundaries, which in our experience can cause instabilities.

To simplify the presentation we restrict ourself to first order time stepping (IMEX Euler) and assume the phase transport/interface reconstruction is given:
Let Γ_h^ n+1 denote the set of all interface reconstruction at time t^n+1 and κ^n+1 the approximated interface mean curvature at time t^n+1.

The complete scheme for the velocity and pressure at time t^n+1 then reads as:
Find u⃗^n+1∈𝒱_k+1^d, p^n+1∈𝒱_k such that

  
    1/Δ tû⃗̂^*-u⃗^nv⃗
        - u⃗^n⊗w⃗^n∇v⃗
        + ∑_e∈Σ_h∫_eF_e(u⃗^n;w⃗^n)·v⃗
        = g⃗v⃗ ,

  
    1/Δ tϱ(u⃗^*-û⃗̂^*)v⃗
        - 2  S(∇_0u⃗^*)∇_0v⃗   = -∇p^nv⃗
        + ∑_γ∈Γ_h^n+1∫_γσκ^n+1v⃗·n⃗_γ ,
    1/ϱ∇p^*∇q   = 1/Δ t∇_0·u⃗^*q ,
    u⃗^n+1 = u⃗^* + Δ t/ϱ∇p^*    ,      p^n+1 = p^n + p^* ,


for all v⃗∈𝒱_k+1^d, q∈𝒱_k.
Here F_e is a suitable numerical flux (e.g. local Lax-Friedrichs flux), w⃗^n = u⃗^n is the transport velocity and all terms in eqs. (<ref>), (<ref>) containing ϱ,  are integrated with cut-cell quadratures with respect to Γ_h^n+1.

We note that here the lifted gradient/divergence maps into the dg space with increased/decreased polynomial order of its argument function space (e.g. ∇_0u⃗∈𝒱_k+2^d× d, ∇_0·u⃗∈𝒱_k,).

The resulting discrete velocity field is exactly divergence-free with respect to the discrete lifted divergence operator for arbitrary density ratios (up to the accurracy of the solver used).
Applying the discrete DG divergence to the velocity field update and using eq. (<ref>)  recovers eq. (<ref>).



§ NUMERICAL EXPERIMENTS


The scheme has been implemented in DUNE-FEM <cit.>, which is part of the DUNE (Distributed and Unified Numerics Environment) framework <cit.>.
The discretization used here are orthonormal 𝒫_2 (resp. 𝒫_1) dg elements for the velocity (resp. pressure) on a cartesian grid for space and a BDF2 scheme for time stepping.
The implicit substeps are solved with a GMRES preconditioned with either an incomplete LU- or LDL^𝖳-decomposition.
The phase transport is solved using a VoF method with a geometric flux, also implemented in DUNE.

We consider the benchmark problems presented in <cit.>.
The general setup is the following:
The computational domain is given by Ω = [0,1]×[0,2] ⊂ℝ^2, Ω_2 is a circular subdomain with diameter d = 0.5 centered at (0.5, 0.5)^⊤ and Ω_1 Ω∖Ω_2.
The no-slip boundary condition is prescribed at the top and bottom boundaries, the free-slip condition is prescribed on the vertical boundaries and the system is initially at reset.

The fluid coefficients for the test problems are given in table <ref>.
Test case 1 results in a ellipsoidal bubble, while in test case 2 the bubble gets deformed significantly and eventually break ups occur.
Our results are in good agreement for the shape (see Fig. <ref>), the rise velocity and the center of mass of the bubble with the results presented in <cit.>.
The `circularity' can not be compared directly since our interface representation makes measuring the perimeter of the bubble difficult, nonetheless it is in general agreement.







Figure <ref> shows an experiment of the effect of droplet merging on a super-hydrophobic surface in a water-steam system in free-fall.
The initial configuaration consist of two equal-sized droplets, two unequal-sized droplets and a single droplet as reference.
When droplets of similar size merge the force of surface tension is strong enough that the droplet “jumps”.




§ CONCLUSIONS AND OUTLOOK


In this contribution we presented a primal local discontinuous Galerkin pressure correction scheme suitable for two-phase flows with high density and viscosity ratios.
The interface is sharply resolved by using cut-cell quadratures.
Numerical experiments show that the scheme agrees with other results presented in the literature.

Further work includes extending the model and the scheme to include the effects of phase transition and investigating if other, more compact, DG methods can recast in the same form as the presented LDG method.


  The first author would like to thank Tobias Malkmus for providing the initial code base for the fluid solver and explaining its subtleties, Martin Nolte for his ideas and help in modifying said code.


spmpsci


