An adaptive fully discontinuous Galerkin level set method for incompressible multiphase flows

Ali Karakus (Department of Mathematics, Virginia Tech, Blacksburg, Virginia, USA)
Tim Warburton (Department of Mathematics, Virginia Tech, Blacksburg, Virginia, USA)
Mehmet Haluk Aksel (Department of Mechanical Engineering, Middle East Technical University, Ankara, Turkey)
Cuneyt Sert (Department of Mechanical Engineering, Middle East Technical University, Ankara, Turkey)

International Journal of Numerical Methods for Heat & Fluid Flow

ISSN: 0961-5539

Publication date: 4 June 2018



This study aims to focus on the development of a high-order discontinuous Galerkin method for the solution of unsteady, incompressible, multiphase flows with level set interface formulation.


Nodal discontinuous Galerkin discretization is used for incompressible Navier–Stokes, level set advection and reinitialization equations on adaptive unstructured elements. Implicit systems arising from the semi-explicit time discretization of the flow equations are solved with a p-multigrid preconditioned conjugate gradient method, which minimizes the memory requirements and increases overall run-time performance. Computations are localized mostly near the interface location to reduce computational cost without sacrificing the accuracy.


The proposed method allows to capture interface topology accurately in simulating wide range of flow regimes with high density/viscosity ratios and offers good mass conservation even in relatively coarse grids, while keeping the simplicity of the level set interface modeling. Efficiency, local high-order accuracy and mass conservation of the method are confirmed through distinct numerical test cases of sloshing, dam break and Rayleigh–Taylor instability.


A fully discontinuous Galerkin, high-order, adaptive method on unstructured grids is introduced where flow and interface equations are solved in discontinuous space.



Karakus, A., Warburton, T., Aksel, M. and Sert, C. (2018), "An adaptive fully discontinuous Galerkin level set method for incompressible multiphase flows", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 28 No. 6, pp. 1256-1278.

Download as .RIS



Emerald Publishing Limited

Copyright © 2018, Emerald Publishing Limited

1. Introduction

Multiphase flows occur in many areas of practical importance such as flow-structure interaction, air–water interfacial dynamics, phase change problems and reacting flows. It is hard to define multiphase flows precisely but it generally denotes the systems in which more than one phase is simultaneously present. Phases refer to fluids of the same substance such as a liquid and its vapor or fluids of different substances such as liquid and gas or fluid and solid, all coexisting in the same domain and at the same time (Andrea and Grétar, 2009). In this work, we investigate a specific multiphase flow, i.e. incompressible and immiscible flow where the phases are separated with a well definite, sharp interface. Numerical prediction of the deformable interfaces in these complex flows are challenging because of discontinuity in the material properties, varying topology and wide range of scales.

Numerous numerical methods were proposed to represent the phase dynamics in incompressible flows. These methods can be classified as interface tracking and interface capturing. The former group is Lagrangian or semi-Lagrangian, where the mesh or massless particles explicitly represent the interface (Unverdi and Tryggvason, 1992; Harlow and Welch, 1965). Interface tracking approaches are generally accurate and robust but difficult to use when the interface encounters topological changes. Interface capturing methods, Volume of Fluid (VOF) (Hirt and Nichols, 1981) and Level Set (LS) (Osher and Sethian, 1988) are Eulerian, i.e. an implicit function defined on a fixed grid represents the interface. Specifically, VOF methods are widely used in multiphase flow simulations (Kleefsman et al., 2005; Caron et al., 2015) because of their natural conservation properties and efficiency, but reconstruction of the interface from volume fraction data and obtaining geometry dependent properties, such as interface normal and curvature, are difficult.

LS method addresses the problems of interface tracking and VOF methods. In the LS method, an interface between the phases is represented by a specified contour level of a continuous function. This implicit interface representation offers many advantages such as straightforward extension from 2D to 3D, simple handling of topological changes and easy calculation of geometric properties. Because of these advantages, LS method is also widely used to model different multiphase flows (Osher and Fedkiw, 2002; Marchandise and Remacle, 2006; Wang et al., 2011; Stapór, 2016). The main difficulty with use of the LS method is the need to control mass loss present in the method. Various techniques have been proposed to improve the conservation properties of the original method. Popular ones are combining LS with VOF (Sussman and Puckett, 2000) or with a semi-Lagrangian particle method (Enright et al., 2002) and applying additional mass correction procedures (Sussman and Fatemi, 1999; Wang et al., 2012). Common to all these approaches is the fact that simplicity and computational efficiency of the original LS method are partly lost. In fact, mass loss problem of the LS method is more inherent in the discretization scheme used than the formulation itself (Marchandise et al., 2006; Karakus et al., 2016a).

The discontinuous Galerkin (DG) methods ((Hesthaven and Warburton, 2008) and references (therein) are a class of finite element methods that make use of completely discontinuous spatial discretization. In DG method, trial and test spaces are restricted to each element, and all communications occur weakly by numerical fluxes, which offer compact stencil size with a substantial flexibility. The DG methods have also excellent dissipation properties achieved by local high-order polynomial approximations to overcome the mass loss problem of the LS formulation. For this reason, using DG for LS advection gives more accurate results compared with the widely used essentially non-oscillatory type schemes (Sussman and Hussaini, 2003; Karakus et al., 2016a). DG LS modeling is also applied in incompressible multiphase flow simulations, i.e. it is combined with a stabilized finite element (Marchandise and Remacle, 2006) and a conservative finite difference (Owkes and Desjardins, 2013) flow solvers. These hybrid methods require interpolation of the velocity field obtained from the lower-order flow solver to higher-order discontinuous polynomial space in each time step. Fine grids are needed for flow solver, which may become computationally inefficient when high-order interpolation orders are used for LS evolution.

Multiphase flows are generally highly dynamic in the sense that location and topology of the interface change considerably during the simulation. In a fixed grid approach, resolution should be increased in the whole domain to capture the full physics of the interface, which obviously creates unnecessary burden in computational time. Adaptive mesh refinement (AMR) techniques decrease the computational effort by increasing the resolution in the vicinity of the interface by keeping coarse grid at less dynamic regions. There are some adaptive LS methods for multiphase flows based on the low-order finite difference/volume schemes on structured Cartesian grids and tree type data structures (Sussman et al., 1999; Sussman, 2005). Unstructured anisotropic mesh refinement is also coupled with the low-order finite element flow solution with higher-order DG interface modeling (Marchandise and Remacle, 2006). In this conformal discretizations, each face can only be shared by two elements. Then, AMR algorithm should handle the complicated interpolation between grids and mesh transition to ensure the conformity. Global re-meshing or iteratively optimized local modifications may also become inefficient when frequent adaptation is needed as in the multiphase flow simulations. Because of relaxed strong elemental connectivity, DG methods can easily be used on unstructured elements with hanging nodes (Kopera and Giraldo, 2014) that enable the fast adaptation with accurate interpolation between the grids.

In this study, we present a high-order, fully DG method for immiscible, incompressible multiphase flows with sharp interfaces on adaptive unstructured grids. The proposed method allows to capture interface topology accurately in simulating wide range of flow regimes with high density/viscosity ratios and offers good mass conservation even in relatively coarse grids while keeping the simplicity of the LS interface modeling. Outline of this paper is as follows: the governing equations of the problem are given in Section 2. Section 3 represents the notation used in the numerical scheme, time and high-order DG spatial discretizations and matrix-free solution techniques. Local interface model with reinitialization of the level set function and AMR strategy are described in Section 4. Finally, Section 5 gives distinct numerical test cases to show accuracy and mass conservation of the method.

2. Governing equations

Incompressible, laminar flow of two immiscible, non-reacting fluids are governed by the Navier–Stokes equations defined on the domain Ω = Ω 1 ∪ Ω2. Here, subscripts 1 and 2 denote the sub-domains of the first and the second fluids, respectively. Domain boundary is represented by Ω while the Γ = Ω 1 ∩ Ω 2 denotes the interface separating the fluid phases. Problem domain, Ω, is fixed but two fluid domains, Ω1 and Ω2, and the interface, Γ, change in time.

The non-dimensional, incompressible Navier–Stokes and the continuity equations on Ω can be written as:

(1) ut+(u·)u=pρ+1Reρ·(2μS)+egFr2+κnΓWe,·u=0,xΩ
where, u denotes the velocity field, p is the static pressure, ρ is the density, μ is the dynamic viscosity, S=12(u+uT) is the deformation rate tensor, eg is the direction of the gravitational acceleration, nΓ is the unit normal of Γ and κ = ∇ · nΓ is the local curvature of the interface. The non-dimensional Reynolds, Froude and Weber numbers are defined using the reference variables denoted with subscript R as Re=ρRURLR/μR,Fr=UR/gLR and We = ρRURLR/σR. Hereafter, F=egFr2 is used as the non-dimensional body force term. In this study, surface tension forces are neglected but the method still governs many practical, large scale problems where curvature is too small to contribute to the governing equations.

LS method is used to model phase dynamics where the interface is represented by at least a Lipschitz continuous function, ϕ, which is positive in one fluid domain and negative in the other. The zero contour of the implicit LS function, i.e. ϕ(x,t) = 0, defines the current location of interface. Using this definition, ρ and μ can be globally defined as:

(2) ρ(ϕ)=ρ1H(ϕ)+ρ2(1H(ϕ))μ(ϕ)=μ1H(ϕ)+μ2(1H(ϕ))
where, H(ϕ) is the Heaviside step function. Time evolution of the interface can be obtained by simple advection of the LS function with the known velocity field, as given below:
(3) ϕt+uϕ=0

In the interface model, LS is kept as signed distance function such that |∇ϕ| = 1 is fulfilled. Evolution of the interface destroys regularity of LS function and creates very large or small gradients near the interface. Reinitialization is required to replace distorted LS with a more desired signed distance function to enhance accuracy of the interface representation and given in section 4.

3. Numerical method

3.1 Preliminaries and notations

We assume that d dimensional domain Ω ⊂ Rd is well approximated by the computational domain, Ωh, having boundaries, Ωh, which can be Dirichlet, ΩD, or Neumann type, ΩN. Ωh is partitioned into non-overlapping, possibly nonconforming triangular/tetrahedral elements, Ωh=k=1KDk. Two element domains, Dk and Dk+, of the triangulation Ωh has a common face if DkDk+, where ∂Dk denotes the element interface. n = –n+ be the unit outward normal vector of ∂Dk. φk and φk+ denote the traces of a scalar function, φ when evaluated at Dk and Dk+, respectively. According to this definition, average and jump operators for a scalar-valued function can be defined as:

(4) {φ}=φ+φ+2,[[φ]]=φφ+

When φ is a vector-valued function, the above operators act on it component wise. On Dirichlet or Neumann boundaries, normal vectors point outward to Ωh, and average and jump operators are equal to the trace of function φ unless otherwise stated.

The discontinuous approximate spaces are VN={vL2(Ω)|v|DkPN(Dk),DkΩh} and its vector version, VNd. PN(Dk) is the space of discontinuous piecewise polynomial functions of degree N ≥ 1 on each elemental domain. Lagrange polynomials are constructed as the basis for this polynomial space using Warp and Blend nodes (Warburton, 2006) and orthonormal Proriol–Koornwinder–Dubiner (Proriol, 1957; Koornwinder, 1975; Dubiner, 1991) polynomials. Let (·,·)Dk represent the inner product computed over the volume of the element k. Similarly, let (·,·)∂Dk denote the inner product taken along the element boundaries.

3.2 Temporal discretization

A high-order splitting (Karniadakis et al., 1991) is used for the temporal discretization of the flow equations. The scheme is semi-implicit, in which non-linear terms are integrated explicitly and linear terms are treated implicitly. Second-order backward differentiation is used for the unsteady term and second-order extrapolation is used for the non-linear terms. With this implementation, equation (1) can be advanced from time levels, n to n + 1 by solving the following semi-discrete equations decomposed into three fractional steps for u:

(5a) u^α0unα1un1Δt=β0((un·)un+Fn)+β1((un1·)un1+Fn1)
(5b) u^^u^Δt=1ρ(ϕ)pn+1
(5c) γ0un+1u^^Δt=1ρ(ϕ)Re·μ(ϕ)un+1

The coefficients, α, β and γ correspond to a second-order stiffly stable scheme, and their values can be found in Karniadakis et al.’s (1991) and Karniadakis and Sherwin’s (2005) studies. In the splitting scheme, pressure is decoupled from the velocity, which avoids spurious pressure modes and enables to use equal order approximations for both flow fields.

In the first step, nonlinear terms and body forces are advanced explicitly by a second-order stiffly stable extrapolation scheme. Then, incompressibility constraint is enforced in equation (5b) by forcing the second intermediate velocity field, u^^ to be divergence free. Finally, a modified Helmholtz equation is solved for the viscous terms to obtain next time-level velocities. In Guermond et al.’s (2006) study, it was demonstrated that the splitting scheme preserves the optimal second-order accuracy.

3.3 Nonlinear treatment

The first step of splitting scheme requires approximation of the nonlinear term, which is written in divergence form as (u · ∇)u = u · ∇u + u∇ · u = ∇ (uu), where uv = uivj,i,j = 1,…,d. Let u be approximated by uhVNd. Multiplying the nonlinear term with test functions, vhVNd, integrating over element domain and performing integration by parts, we obtain the following local statement.

(6) (·(uhuh),vh)Dk=(uhuh,·vh)Dk+(n·[[uhuh]],vh)Dk

Because of discontinuous approximation space, flux function, n · [[uhuh]], in the boundary contribution is not uniquely defined and hence replaced with the local Lax–Friedrich numerical flux, F͂c

(7) F˜C=n·{uhuh}+12ΛDk[[uh]]
here, ΛDk=maxu|n·(uhuh)uh| is the absolute value of maximum eigenvalue of the flux Jacobian. The first term of the flux function denotes central part and the second one is the dissipative contribution. ΛDk determines the required artificial diffusion to stabilize the system. This flux choice leads to compact stencil size such that the degrees of freedom of an element couples only with its immediate neighbors. Dirichlet boundary conditions are applied weakly by defining an exterior ghost state where average jump operators become {uhuh}DkΩD=0.5(uhuh+uDuD) and [[uh]]DkΩD=(uhuD).

Explicit treatment of nonlinear terms introduces a CFL type restriction on the time step size. CFL estimate for high-order methods gives Δt ≈ O (h/UN2) for advection model problems, with U and h referring to the characteristic velocity and mesh size, respectively.

3.4 Implicit treatment

The time splitting scheme requires solution of implicit pressure Poisson equation for velocity projection and modified Helmholtz equation for velocity correction. Especially, Poisson equation needs special attention because of poor conditioning of the system for efficient solution and obtaining divergence free velocity field for long-term stability.

In the splitting scheme, incompressibility constraint is enforced by taking the divergence of equation (5b), which leads to the following variable density pressure Poisson equation:

(8) ·(u^Δt)=·(1ρ(ϕ)pn+1)

Boundary conditions of this Neumann problem can be derived from the original equation. To preserve the temporal accuracy and fulfill the compatibility condition, Neumann boundary conditions have to be approximated to the same order as the time discretization scheme, which requires extrapolation of the Neumann boundary data:

(9) n·(ρ1p)n+1=i=01βi(ut+(u·)u(ρRe)1·(2μS)F)ni·n

There are different DG methods developed for the discretization of the equations with second-order (or higher-order) operators, most of which are based on the mixed formulation (Arnold et al., 2002). In this study, symmetric interior penalty (SIP) (Arnold, 1982) method is preferred because of its simplicity, computational efficiency and compact stencil size. Let pressure p be approximated by phVN, and vhVN be the test functions. SIP discretization of equation (8) for an elemental domain reads as follow:

(10) (ρ1ph,vh)Dk(n·{ρ1ph},vh)DkΩN([[ph]],n·ρ1vh)+(γ[[ph]],vh)DkΩN=(·u^Δt,vh)Dk(pD,n·ρ1vh)DkΩD+(γpD,vh)DkΩD+(pN,vh)DkΩN
where, γ is the penalty parameter and pD and pN are Dirichlet and derived Neumann boundary conditions. In the equation, time dependence of the pressure is dropped for clarity.

Penalty parameter plays an important role for the stability of scheme and should be selected sufficiently large to enforce coercivity. On the other hand, selecting arbitrarily large γ values increases condition number of the system and degrades performance of the linear solvers. Penalty parameter is defined as:

(11) γDk=12(N+1)(N+2)max(hDkρDkhDk,hDk+ρDk+hDk+)
where, h∂Dk and hDk denote surface area and volume of the element, respectively. On the non-conformal faces with hanging nodes, where more than two elements are connected, γ is computed using the collection of all adjacent elements, which guarantees coercivity of the bilinear form (Shahbazi, 2005). Solution of the projection step is challenging and time consuming because of construction of global system matrices in each time step. Also, discretization matrix of the Poisson problem is poorly conditioned in high density ratio flows even with optimum penalty parameter selection.

The pressure system arising from the SIP discretization can be written in the generic form, Ax = b, where A is the symmetric and positive definite coefficient matrix that enables the solution of the system efficiently using the preconditioned conjugate gradient method. To avoid costly set-up time and minimize the memory requirements, no global matrices are constructed in the conjugate gradient (CG) solver. Instead, matrix vector multiplications are performed explicitly by doing the calculations per-element as given in equation (10). Because the system is poorly conditioned in high-density ratio flows, efficient preconditioners are required to converge the system within small iteration numbers.

Because of its known advantages, algebraic multigrid (AMG) is used as the CG preconditioner in this study. AMG has set-up and solve phases. Set-up phase includes construction of hierarchical grids and defining restriction and interpolation operators to transfer data between fine to coarse and coarse to fine grids, respectively. Hierarchy of the different grids can be computed at the beginning of the solution and stored for fixed grid and constant coefficient problems, where sparsity pattern and spectral properties of the coefficient matrix do not change in time. But this is not the case for adaptively refined meshes, where computational burden of the set-up phase of the preconditioner may suppress the actual solution time.

A p-multigrid preconditioner is designed for the pressure Poisson system to reduce the computational time and memory requirements. Construction of the full matrix and related hierarchical levels are avoided. Required residuals are computed with a matrix-free approach starting from the approximation order of N until reaching first-order coarse system. Restriction operator, Ipq that projects the errors from the polynomial space of order p to a space of order q for p > q is easily constructed because of discontinuous interpolating functions:

(12) Ipq=(Mq)1Mqp
where corresponding matrices can be computed in index notation as:
(13) Mi,jq=(viq,vjq)RandMi,jqp=(viq,vjp)R

Note that restriction operator is defined on the reference element and can be applied element-wise through the whole domain. Prolongation operator that transforms state vectors from the low-order polynomial space to the higher one can be defined similarly, leading Iqp=(Ipq)T. With the use of Lagrange basis space, these operators are full but small matrices in the size of polynomial space dimension, Nq × Np, implying required matrix-vector multiplication can be performed efficiently for an element. Because residuals are computed in matrix-free form, we do not need to define residual restriction/prolongation operators where they can be different from the state vector operators. Instead, residual functions handling different approximation orders are constructed. Geometric mappings and metric identities are constant for straight sided elements, which enables the computation of residuals without storing extra information.

After reducing approximation order to the first order, q = 1, the system matrix is assembled, and related hierarchical grids are constructed at this level based on the recently introduced AMG method (Gandham et al., 2014). One pre- and one post-smoothing step, Sl are used to remove high-frequency errors at the corresponding level, l. Damped Jacobi smoother given below is used because of its comparatively low set-up cost:

(14) x=x+ωD1(bAx),ω=431ϱ(D1A)
where, D is the diagonal of A and ϱ(D−1A) is the spectral radius of the matrix D– 1A estimated from the Arnoldi iterations. K-multigrid cycles are used as the preconditioner for CG iterations as illustrated in Algorithm 1. Number of coarse level iterations increases with K-cycles but it is required to achieve grid-independent convergence and improve overall run-time performance. Note that combining matrix-free p-multigrid with AMG reduces the storage and computational effort significantly. If the system is solved with matrix form, NN × NN × (1 + Nf) non-zero entries have to be stored for an element only at the highest level, where Nf stands for the total face connection pairs. This requirement is reduced to N1×N1×(1 + Nf) so that matrix free approach becomes more effective in 3D problems and high orders.

      Algorithm 1 K-Cycle Multigrid Preconditioner      

1: xl,p ← K-Cycle ((l,p), bl,p, xl,p)      

2: Input: initial guess x0, order n, level l      

3: Output: Updated Solution xlN      

4: if q ≥ 1 then {Operations in Matrix-free Form}      

5: xpSp(bp, Ap, xp) {pre-smoothing}      

6: rqbqAqxq {compute residual for 1 < nN}      

7: rqIpqrq {restrict residual to coarse-grid}      

8: xq ← K-Cycle (q, rq, xq) {inner CG iteration}      

9: else {Operations in Matrix Form}      

10: xlSp(lp, Al, xl) {pre-smoothing}      

11: rlblAlxl {compute residual in matrix form}      

12: rl+1Ill+1rl {restrict residual to coarse-grid}      

13: if l + 1 = L then {coarsest grid exact solution}      

14: xl+1Al+11rl+1      

15: else      

16: cl+1 ← K-Cycle(l + 1, rl+1, xl+1) {inner CG iteration}      

17: vl+1Al+1cl+1, a1cl+1Tvl+1      

18: α1cl+1Trl+1,r˜l+1rl+1α1a1vl+1      

19: ifr͂l + 1‖ ≤ TOL‖rl + 1then      

20: xl+1α1a1cl+1      

21: else      

22: dl+1K-Cycle(l+1,r˜l+1,xl+1) {inner second CG iteration}      

23: wl+1Al+1dl+1,γdl+1Tvl+1      

24: βdl+1Twl+1,α2dl+1Tr˜l+1,a2βγ2a1      

25: xl+1(α1a1γα2a1a2)cl+1+α2a2dl+1      

26: if l ≠ 0 then {Operations in Matrix Form}      

27: rlIl+1lrl+1 {prolongate}      

28: rlrlAlxl {compute residual}      

29: xlSp(rl, Al, xl) {post-smoothing}      

30: else {Operations in Matrix-Free Form}      

31: rpIqprp {prolongate}      

32: rprpApxp {compute residual}      

33: xpSp(rp,Ap,xp) {post-smoothing}

After obtaining the next time-level pressure, phn+1, second intermediate velocity field is computed with the use of equation (5b). Required pressure gradient, G=phn+1, is evaluated in the weak form:

(15) (G,vh)Dk=(phn+1,vh)Dk+(n·{phn+1},vh)Dk
where the flux function in the surface integral contribution is replaced with the numerical central flux.

Finally, time discretization of the flow equations is completed by applying viscous correction through the solution of a modified Helmholtz equation given in equation (5c) which can be recast into following form:

(16) γ0ρReΔtun+1·μun+1=ρReΔtu^^¯

Actually, the equation is composed of scalar Helmholtz equations for each velocity component closed with the appropriate velocity boundary conditions at time, tn+1. For the scalar form, this system is very similar to pressure Poisson equation with additional scaled un+1 term on the left-hand side. Let u be a component of the velocity vector approximated by uhVN, and vhVN be the test functions. Similar to the pressure equation, SIP discretization of an elemental domain can be written as:

(17) (γ0ρReΔtuh,vh)Dk+(μuh,vh)Dk(n·{μuh},vh)DkΩN([[uh]],n·μvh)+(γ[[uh]],vh)DkΩN=(ρReΔtu^^¯,vh)Dk(uD,n·μvh)DkΩD+(γuD,vh)DkΩD+(uN,vh)DkΩN
where penalty parameters and boundary conditions are as defined before. The system is symmetric, positive definite and solved with CG method in matrix-free form. Because solution of velocity system is considerably easier than the pressure equation because of scaled mass matrix, block-Jacobi preconditioner is preferred.

4. Interface modeling

4.1 Discontinuous fluid properties

Discontinuous fluid properties are avoided by smoothing their variations in the vicinity of interface with thickness of size, ϵ. There are many definitions of smoothing functions, such as piecewise continuous (Grooss and Hesthaven, 2006) or continuous trigonometric (Sussman et al., 1999). Here, hyperbolic tangent function is used as the regularized Heaviside function:

(18) Hϵ(ϕ)=tanh(πϕϵ)
which smooths variations over the distance 2ϵ. Heaviside function is infinitely differentiable and more appropriate for high-order discretizations.

Interface thickness, ϵ, should be selected as small as possible to get sharp profile but large enough to capture variations and to stabilize the system. Classical low-order methods require to choose ϵ at least in the order of characteristic mesh size, ϵ = O(h). Because of high-order interpolation of the DG method, variations can be resolved and integrated stably in the smaller distance of ϵ = O(h/N). Our numerical tests showed that ϵ = 2h/N offers sharp but smooth profiles with good integrability of the equations.

After defining the smoothed Heaviside function, discontinuous material properties can be replaced with globally defined continuous counterparts. For the isothermal, incompressible flows, only density and viscosity are to be smoothed as follows:

(19) ρ(ϕ)=12ρ1(1+Hϵ(ϕ))+12ρ2(1Hϵ(ϕ))μ(ϕ)=12μ1(1+Hϵ(ϕ))+12μ2(1Hϵ(ϕ))
which gives, ρ(ϕ) = ρ1 for ϕ < 0 and ρ(ϕ) = ρ2 for ϕ > 0 with smooth transition around ϕ = 0. LS function has to be symmetric with respect to the zero-level set to accurately represent the material properties, which requires constant gradient at least around the ϵ-neighborhood of the interface which is achieved by the use of signed distance function.

4.2 Level set advection

For incompressible velocity field, using the identity, ∇ · (uϕ) = ϕ∇ · u + u · ∇ϕ, LS equation can be written in the following conservative form:

(20) ϕt+(uϕ)=0

Similar to the previous section, let ϕ be approximated by ϕhVN, and vhVN are smooth test functions. Weak form of equation (20) can be obtained by integration by parts:

(21) (ϕht,vh)Dk(uhϕh,vh)Dk+(F˜L,vh)Dk=0
where F͂Lh,uh) is the numerical flux to approximate the normal trace of non-uniquely defined flux function, n · [[uhϕh]]. For F͂L, upwinding is used to approximate ϕh on ∂Dk, such that:
(22) F˜L(ϕh,uh)=(nuh){ϕh}+12|nuh|[[ϕh]]

LS equation, hence the interface location, is evolved in time with an explicit TVD Runge–Kutta (RK) method. Stability region of the RK method is larger than the Adams–Bashforth method used in the explicit integration of the flow equations. This guarantees the stability of LS advection under the time step restriction of the semi-explicit flow solution.

DG method has excellent mass conservation properties because of the low numerical dissipation achieved by the use of local high-order polynomial approximations (Karakus et al., 2016b). So that, the formulation enables to capture the motion of interface accurately without using any special treatments or modifications of the original LS formulation.

4.3 Reinitialization

While advancing the LS equation in time, level sets adjacent to interface may move with velocities different from the zero LS which distorts the scalar LS field. Losing the regularity and symmetry of the LS function around the interface may cause numerical instabilities, errors in the computation of fluid properties and precisely locating interface location. LS function should be reinitialized to signed distance at least around the interface. But, it is not required to apply reinitialization at the end of each time step to improve the mass loss unlike the standard LS methods (Sussman et al., 1994; Smolianski, 2005) because of the high-order, conservative interface modeling. We reinitialized the LS function only when it produces too steep or flat profile near the interface based on the indicator relying |∇ϕ|.

The flow-based reinitialization (Sussman et al., 1994) given below is used here because of its flexibility, high-order accuracy and efficiency:

(23) ϕτ+sgnϵ(ϕ0)(|ϕ|1)=0,ϕ(x,0)=ϕ0
where, τ is the pseudo time which is not related with the physical time, t and sgnα is the regularized signum function. Equation (23) defines an artificial flow to obtain signed distance function in steady state. We choose to regularize sign function similar to Heaviside function:
(24) sgnϵ(ϕ0)=tanh(πϕ0|ϕ0|ϵ)
where, band thickness, ϵ, is selected as in the smoothing material properties. Because LS function is not reinitialized in each time step, it is not exactly signed distance function in the whole computation time. The scaling factor, |∇ϕ0|, is added to avoid very small coefficients resulting in small characteristic speeds, when the LS function becomes flat or to improve the accuracy when the LS function is too steep around the interface.

After regularizing the signum term, reinitialization equation can be written as the following Hamilton–Jacobi form with smooth coefficient:

(25) ϕτ+H(ϕ,x)=0,ϕ(x,0)=ϕ0
where, the Hamiltonian, H(∇ϕ,x) denotes sgnϵ(ϕ0)(|∇ϕ| – 1) and x = xi and i = 1, … , d are the Cartesian coordinates. In Karakus et al. (2016b), we introduced a direct, adaptive, high-order DG method for the LS reinitialization with artificial diffusion stabilization which preserves the optimal accuracy. Our technique is updated here for accounting the local interface modeling for efficient multiphase flow simulations.

Solution of HJ equation requires accurate approximation of derivatives. We followed the local DG scheme (Yan and Osher, 2011) to approximate ∇ϕ by defining two new variables for each derivative component and applied the standard upwind DG method. Let pxl and pxr be auxiliary variables used to approximate ϕxi when fluxes are chosen from the left and right upwinding sides such that, pxil,rϕxi=0. Multiplying the relation with smooth test functions, vh, and integrating by parts gives the following upwind DG scheme in weak form:

(26) (pxil,r,vh)Dk+(ϕh,vhx)Dk(FRil,rnxi,vh)Dk=0
the left, FRil, and right, FRir, upwind fluxes associated with the ith component are defined as:
(27) FR,il={ϕ+fornxi0ϕfornxi<0andFR,ir={ϕfornxi0ϕ+fornxi<0

Although two auxiliary variables are defined for each component, it is obvious that volume terms are the same and only the flux functions differ between directional derivative approximations. Finally, equation (25) is integrated over the cell Dk to get:

(28) (ϕht,vh)Dk+(H¯(pxl,pxr),vh)Dk=0

H¯(pxl,pxr) is the Godunov numerical Hamiltonian approximating H(∇ϕ,x). For a piecewise constant approximation, this scheme is monotone and converges to the entropy solution. However, stabilization is needed for higher order approximations. An artificial diffusion approach is used to damp out the high frequencies in the solution (Karakus et al., 2016b).

Characteristic velocities of the reinitialization equation, as mentioned before, point outwards from the interface in the direction of normals. In other words, LS function is reinitialized to signed distance starting from interface. Multiphase flows require LS function to be signed distance around interface only. This means that we do not need to solve the reinitialization equation to steady state over the whole domain instead over small neighborhood of interface. We kept signed distance function on the thickness of 2ϵ which requires approximately 2ϵ/Δτ time steps. Because characteristic speed is slightly lower than 1 because of regularization of the signum function, we use more time steps than the estimated value. Numerical tests show that 10 per cent extra time step is satisfactory to obtain the signed distance function in the band thickness. Constant LS function is used outside of interface band. In practice, a few explicit RK step is satisfactory to get the signed distance function around the interval so that reinitialization does not create any computational burden in total simulation time.

4.4 Mesh adaptivity

A non-conformal discretization on unstructured triangular/tetrahedral elements is used in this study. The refinement strategy enables to obtain fast and local adaptation and accurate interpolation of field variables between the grids. In the adaptive scheme, the computational mesh consists of elements in a range of predefined levels with l0 denoting initial coarse level and lM being the maximum level. Refinement and coarsening are performed dynamically during the solution. The level of refinement and the elemental dependencies are stored in a hierarchical tree for efficient h type adaptivity.

Refinement is carried out in an isotropic way, i.e. a parent element is divided into siblings by connecting the mid-edges resulting with four and eight children for a triangle and a tetrahedron, respectively. A threshold value, γ, is selected to mark the elements for refinement. If min|ϕk| ≤ γ holds and refinement level of the element is smaller than the predefined maximum refinement level (lk < lM), the element is marked for refinement and the approximation on the parent element is projected onto its siblings. If min|ϕk| > γ holds and refinement level of all four siblings are larger than the initial coarse level, these elements are marked for coarsening and combined by removing all siblings and their newly created vertices. Solution of the siblings are then projected to their parent element and level of refinement information is updated. Projection from and to siblings are local operations, which can be performed efficiently by local matrix–vector multiplications.

The DG method supports arbitrary number of hanging nodes per face, but it is restricted to decrease computational complexity in the flux evaluations (Karakus et al., 2016a). Adaptive mesh is 3:1 and 2:1 balanced for triangular and tetrahedral elements, respectively. In other words, a face of a triangular and a tetrahedral element can connect to four and three elements at most, respectively. Any 1:1 connection (conformal pair) is a member of 2:1 balanced grid and so on.

5. Numerical tests

Verification problems in two and three spatial dimensions are solved to investigate convergence and mass conservation characteristics of the developed numerical framework for multiphase flows with various density and viscosity ratios and interface topologies.

5.1 Poiseuille two-phase flow

We first consider the flow of two immiscible fluids between parallel plates where exact solution is known at steady state. This test case does not include the interface deformation but helps to characterize errors in the flow field (Marchandise and Remacle, 2006). Density of fluids is the same, whereas kinematic viscosities are ν1 = 0.1 and ν2 = 0.02. The problem is solved on 2D computational domain of [0,4] × [–1,1]. Interface between the fluids is located at the half of the channel height. To ensure the constant gradient of –2, pressure is imposed on the left and the right boundaries of the domain, while standard wall conditions are used on the upper and the lower boundaries.

We performed convergence studies for both successive mesh refinements and approximation order enrichment. Figure 1(a)-(b) illustrate the exact and numerical computed velocity field along x = 0 line on three different semi-structured meshes with the characteristic sizes of h, h∕2 and h∕4 for h = 0.1 and polynomial order, N = 3. It can be seen that the numerical scheme captures the analytical solution well in all cases. The zoomed view shows that errors are larger around the interface as expected. Increasing resolution improves the quality of approximation, and numerical result matches with the exact solution. Finally, Figure 1(c) shows L error on velocity field for the same meshes and approximation orders of N = 3, N = 4 and N = 5. Numerical results show spectral convergence which is around the optimal N + 1 rate for this undeformed interface problem.

5.2 Sloshing in rectangular tank

This test case has a liquid column oscillating in a rectangular tank, with an incompressible gas phase over it. The driving forces are gravitational acceleration and viscous dissipation, causing the liquid column eventually reach its lowest potential energy level from the initially perturbed shape, which is given by the following sinusoidal function:

(29) Γ=d+a0sin(π(0.5dx))
where d is the mean liquid depth and a0 is the amplitude of the initial wave. Reference velocity, UR and length, LR for this problem are gd and d, respectively.

The problem is first solved in the computational domain of [0, d] × [0, 2d] for the small amplitude, a0 = 0.01, where an analytical solution based on the linearized Navier–Stokes equations is reported by Wu et al. (2001). The initial grid has elements with a characteristic element length of h = 1∕3 (K = 50), which is refined adaptively near y = 1 with lM = 3 as shown in Figure 2(a)-(b). Slip boundary condition is assigned to the bottom and side walls while zero pressure is imposed at the top wall. Zero velocities and hydrostatic pressure distribution are used as the initial condition. Density and viscosity ratios of the liquid and gas phases are both 100. The non-dimensional Froude and Reynolds numbers are taken as 1 and 100, respectively. Figure 2(c) shows variation of the normalized wave elevation, ηa0 in time for computed and analytical solutions. Numerical result matches well with the analytical one.

Figure 3 shows the lM = 3 locally adapted mesh structures and interfaces at different simulation times for the high initial amplitude case with a0 = 0.2. Polynomial order is set to five and all other parameters are kept the same as the previous low amplitude solution. Time integration is carried out until t = 40, where the liquid column comes to rest. Mesh adaptivity used here always keeps the interface in the highest refinement level elements.

Table I shows the liquid phase mass fluctuations introduced by the scheme for long time integration. Different refinement levels and order of approximations are used in the numerical experiments. Mass fluctuations are computed very accurately by the adaptive contouring algorithm (Remacle et al., 2007) using the formula, Af = 100 × (max (At) – min (At))/Aexact, where At is the time history of the area of the liquid phase. Mass is well conserved for this long time integration solution and increasing resolution in the vicinity of interface significantly improves it.

5.3 Rayleigh–Taylor instability

In this commonly used test problem, a heavy fluid overlies a lighter one. Instability arises when the initial interface is perturbed. Because of the vertical gravitational field, fluids start penetrating into each other with increasing amplitude in time. Problem is solved in the 2D rectangular domain of [0, d/2] × [0, 4d], with the characteristic length d being 1. Boundary conditions are slip at the side walls, no-slip at the bottom wall and prescribed zero pressure at the top wall. Initial conditions are taken as zero velocity field and hydrostatic pressure distribution. Interface between the fluids is initially perturbed with a cosine function of amplitude 0.1, η = 2.0d + 0.1cos (2πx). Dynamic viscosities of the fluids are equal and their density difference is given by the Atwood ratio, At = (ρlρg)/(ρl + ρg). To use the same notation as used by Tryggvason (1988), reference time is chosen as tr=d/(gAt).

Non-dimensional tip positions of the rising and dropping fluid columns, xt, are given in Figure 4 for Re = ρld3/2g1/2/μl and At = 0.5. For comparison, results obtained by the Lagrangian–Eulerian vortex method (Tryggvason, 1988) and the variable density finite element projection method (Guermond and Quartapelle, 2000) is also included. Numerical results are obtained for a fixed grid with characteristic length h = 1∕3 (K = 212) and locally adapted grid with lM = 2. As seen, present results compare well with the previous works and increasing resolution near the interface improves the accuracy.

Figure 5 shows the interface shapes for lM = 2 locally adapted grid at different simulation times. Results compare well with those reported by Puckett et al. (1997), Popinet and Zaleski (1999) and Xie et al. (2014). Mesh adaptation details for Re = 500 solution are shown in Figure 6. Effective local adaptation of the initial coarse mesh around the complex interface is clearly seen. At the shown instant, element numbers for fixed and adapted grids with lM = 1 and lM = 2 are 212, 305 and 425, respectively.

Percentage mass fluctuations of the liquid phase of the Rayleigh–Taylor instability problem are given in Table II. Values are computed as described previously for the sloshing problem. Starting with the same initial coarse grid, simulations are conducted for different approximation orders and refinement levels. Results obtained for the same problem with a front tracking method (Popinet and Zaleski, 1999) and stabilized finite element discontinuous LS formulation (FEM-DGLS) (Marchandise and Remacle, 2006) are also included in the table. In adaptive simulations, time average of the element numbers are used for the comparison with the reference fixed-grid solutions. Use of adaptive grid improved the mass loss problem significantly, and very accurate results are obtained with increasing refinement level. It is worthwhile to mention that the degree of freedom for the finest solution (lM = 2, N = 5) is very close to those of the reference studies, but lower mass loss is achieved with the current approach.

5.4 3D dam break with obstacle

To test the developed solver on 3D, complex interface problems, a dam break problem with a rectangular obstacle is selected. For direct comparison with the experimental study of Maritime Research Institute Netherlands (MARIN) (SPHERIC, 2016), computational domain is taken as a box of size 3.22 × 1 × 1 m. Rectangular obstacle with dimensions 0.161 × 0.403 × 0.161 m is placed 2.476 m downstream of the water column. Rectangular water column with height and width of 0.55 m and 1.228 m is released at t = 0 and collapses under the action of gravity, which creates a highly unsteady flow field and complex interface deformations.

Slip boundary conditions are applied to the bottom and side walls. Top wall is modeled as open boundary with zero atmospheric pressure and normal velocity gradients. Water and air properties are assigned to the liquid and gas phases as ρl = 1,000 kg m−3, ρg = 1 kg m−3 and μl = 10−3Pa s, μg = 10−5Pa s. The problem is first solved on a fixed tetrahedral grid with characteristic length of 0.2 (K ≈ 20000).

Figure 7 shows snapshots of the interface at different solution times on lM = 1 locally adapted grid. The figure demonstrates how the developed solver is capable of tracking highly deformable interfaces with strong topological changes.

In the experimental work of MARIN (SPHERIC, 2016), pressure histories on the rectangular obstacle are measured by eight sensors, four located on the front wall and four on the top wall. Figure 8 presents computed pressure history on the first and third sensor location (SPHERIC, 2016) for fixed and lM = 1 grids. The instant that the water wave hits the obstacle is measured as t = 0.4 s, which is well approximated by the current solution. Simulation and the experiment show good agreement except the magnitude of the impact pressure which is under-predicted by the numerical solution, especially with the fixed grid. This difference was also observed in other reference solutions using different techniques such as VOF (Kleefsman et al., 2005) and smoothed particle hydrodynamics (Lee et al., 2010).

6. Conclusion

A high-order, fully discontinuous Galerkin method for the solution of incompressible multiphase flows on unstructured adaptive meshes is presented. With the proposed numerical framework, the mass is well-conserved even in the coarse grid solutions without introducing any special treatment or modification of the level set formulation. Velocity, pressure and interface modeling uses the same high-order polynomial space preventing interpolation of the field variables, increasing the efficiency. All implicit systems are solved with a matrix-free, p-multigrid approach, reducing the memory requirement. As a future work, parallelization of the method on all GPU–CPU platforms will be considered using recently created unified multi-threading language, OCCA (Medina et al., 2014).


Comparison between analytical and numerical solutions computed on successively refined meshes for two-phase Poiseuille flow

Figure 1.

Comparison between analytical and numerical solutions computed on successively refined meshes for two-phase Poiseuille flow

Sloshing problem for the small amplitude case

Figure 2.

Sloshing problem for the small amplitude case

lM = 3 Locally adapted grid and the interface at different simulation times for the high amplitude sloshing

Figure 3.

lM = 3 Locally adapted grid and the interface at different simulation times for the high amplitude sloshing

Tip positions of the rising and dropping fluids for Rayleigh–Taylor instability

Figure 4.

Tip positions of the rising and dropping fluids for Rayleigh–Taylor instability

Interface evolution of the Rayleigh–Taylor instability for lM = 2 adaptive solution at (a) t = 0; (b) t = 1.28; (c) t = 1.71; (d) t = 2.18 and (e) t = 2.70

Figure 5.

Interface evolution of the Rayleigh–Taylor instability for lM = 2 adaptive solution at (a) t = 0; (b) t = 1.28; (c) t = 1.71; (d) t = 2.18 and (e) t = 2.70

Interface shapes of the Rayleigh–Taylor instability for (a) fixed grid (b) lM = 1 and (c) lM = 2 locally adapted grids

Figure 6.

Interface shapes of the Rayleigh–Taylor instability for (a) fixed grid (b) lM = 1 and (c) lM = 2 locally adapted grids

Snapshots of interface topology for 3D broken dam problem for lM = 1, N = 3 and h = 0.2

Figure 7.

Snapshots of interface topology for 3D broken dam problem for lM = 1, N = 3 and h = 0.2

Comparison of computed pressure histories with experimental measurements (SPHERIC, 2016)

Figure 8.

Comparison of computed pressure histories with experimental measurements (SPHERIC, 2016)

Percentage mass fluctuations of the liquid phase for the sloshing problem obtained using different refinement levels and approximation orders

N lM
0 1 2
3 4.52 × 10– 1 1.11 × 10– 1 2.12 × 10– 2
4 1.51 × 10– 1 8.27 × 10– 2 6.60 × 10– 3
5 6.33 × 10– 2 1.55 × 10– 2 2.01 × 10– 3

Percentage mass fluctuations of the liquid phase of the Rayleigh–Taylor instability problem for different refinement levels and approximation orders (Re = 500, At = 0.5)

Method lM N Element No. % Mass loss
Present method 0 3 212 0.259
5 212 0.148
3 836 0.093
5 836 0.032
1 3 296 0.096
5 294 0.037
2 3 490 0.025
5 482 0.0091
Front Tracking [45] 32 × 265 0.14
FEM-DGLS [9] 1 32 × 265 0.17


Andrea, P. and Grétar, T. (Eds) (2009), Computational Methods for Multiphase Flow, 1st ed., Cambridge University Press, Cambridge; New York.

Arnold, D. (1982), “An interior penalty finite element method with discontinuous elements”, SIAM Journal on Numerical Analysis, Vol. 19 No. 4, pp. 742-760.

Arnold, D., Brezzi, F., Cockburn, B. and Marini, L.D. (2002), “Unified analysis of discontinuous Galerkin methods for elliptic problems”, SIAM Journal on Numerical Analysis, Vol. 39 No. 5, pp. 1749-1779.

Caron, P.A., Cruchaga, M.A. and Larreteguy, A.E. (2015), “Sensitivity analysis of finite volume simulations of a breaking dam problem”, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 25 No. 7, pp. 1718-1745.

Dubiner, M. (1991), “Spectral methods on triangles and other domains”, Journal of Scientific Computing, Vol. 6 No. 4, pp. 345-390.

Enright, D., Fedkiw, R., Ferziger, J. and Mitchell, I. (2002), “A hybrid particle level set method for improved interface capturing”, Journal of Computational Physics, Vol. 183 No. 1, pp. 83-116.

Gandham, R., Esler, K. and Zhang, Y. (2014), “A GPU accelerated aggregation algebraic multigrid method”, Computers & Mathematics with Applications, Vol. 68 No. 10, pp. 1151-1160.

Grooss, J. and Hesthaven, J.S. (2006), “A level set discontinuous Galerkin method for free surface flows”, Computer Methods in Applied Mechanics and Engineering, Vol. 195 Nos 25/28, pp. 3406-3429.

Guermond, J.L. and Quartapelle, L. (2000), “A projection FEM for variable density incompressible flows”, Journal of Computational Physics, Vol. 165 No. 1, pp. 167-188.

Guermond, J.L., Minev, P. and Shen, J. (2006), “An overview of projection methods for incompressible flows”, Computer Methods in Applied Mechanics and Engineering, Vol. 195 Nos 44/47, pp. 6011-6045.

Harlow, F.H. and Welch, F.E. (1965), “Numerical calculations of time dependent viscous incompressible flow of fluid with a free surface”, Physics of Fluids, Vol. 8 No. 12, pp. 21-82.

Hesthaven, J.S. and Warburton, T. (2008), Nodal Discontinuous Galerkin Methods: algorithms, Analysis, and Applications, Springer, Berlin.

Hirt, C.W. and Nichols, B.D. (1981), “Volume of fluid (VOF) method for the dynamics of free boundaries”, Journal of Computational Physics, Vol. 39 No. 1, pp. 201-225.

Karakus, A., Warburton, T., Aksel, M.H. and Sert, C. (2016a), “A GPU-accelerated adaptive discontinuous Galerkin method for level set equation”, International Journal of Computational Fluid Dynamics, Vol. 30 No. 1, pp. 56-68.

Karakus, A., Warburton, T., Aksel, M.H. and Sert, C. (2016b), “A GPU accelerated level set reinitialization for an adaptive discontinuous Galerkin method”, Computers & Mathematics with Applications, Vol. 72 No. 3, pp. 755-767.

Karniadakis, G. and Sherwin, S. (2005), Spectral/Hp Element Methods for CFD, Oxford University Press.

Karniadakis, G., Israeli, M. and Orszag, S.A. (1991), “High-order splitting methods for the incompressible Navier-Stokes equations”, Journal of Computational Physics, Vol. 97 No. 2, pp. 414-443.

Kleefsman, K.M.T., Fekken, G., Veldman, A.E.P., Iwanowski, B. and Buchner, B. (2005), “A volume-of-fluid based simulation method for wave impact problems”, Journal of Computational Physics, Vol. 206 No. 1, pp. 363-393.

Koornwinder, T. (1975), “Two-variable analogues of the classical orthogonal polynomials”, Theory and Application of Special Functions, pp. 435-495.

Kopera, M.A. and Giraldo, F.X. (2014), “Analysis of adaptive mesh refinement for IMEX discontinuous Galerkin solutions of the compressible Euler equations with application to atmospheric simulations”, Journal of Computational Physics, Vol. 275, pp. 92-117.

Lee, E.S., Violeau, D., Issa, R. and Ploix, S. (2010), “Application of weakly compressible and truly incompressible SPH to 3-D water collapse in waterworks”, Journal of Hydraulic Research, Vol. 48 No. 1, pp. 50-60.

Marchandise, E. and Remacle, J.F. (2006), “A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows”, Journal of Computational Physics, Vol. 219 No. 2, pp. 780-800.

Marchandise, E., Remacle, J.F. and Chevaugeon, N. (2006), “A quadrature-free discontinuous Galerkin method for the level set equation”, Journal of Computational Physics, Vol. 212 No. 1, pp. 338-357.

Medina, D., St-Cyr, S.A. and Warburton, T. (2014), OCCA: A Unified Approach to Multi-Threading Languages, arXiv,1403.0968.

Owkes, M. and Desjardins, O. (2013), “A discontinuous Galerkin conservative level set scheme for interface capturing in multiphase flows”, Journal of Computational Physics, Vol. 249, pp. 275-302.

Osher, S. and Fedkiw, R.P. (2002), Level Set Methods and Dynamic Implicit Surfaces, 1st ed., Springer, Berlin.

Osher, S. and Sethian, J.A. (1988), “Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-jacobi formulations”, Journal of Computational Physics, Vol. 79 No. 1, pp. 12-49.

Popinet, S. and Zaleski, S. (1999), “A front-tracking algorithm for accurate representation of surface tension”, International Journal for Numerical Methods in Fluids, Vol. 30 No. 6, pp. 775-793.

Proriol, J. (1957), “Sur une famille de polynomes ‘a deux variables orthogonaux dans un triangle”, Comptes rendus de l’Académie des Sciences, Vol. 257 No. 1, pp. 2459-2461.

Puckett, E.G., Almgren, A.S., Bell, J.B., Marcus, D.L. and Rider, W.J. (1997), “A high-order projection method for tracking fluid interfaces in variable density incompressible flows”, Journal of Computational Physics, Vol. 130 No. 2, pp. 269-282.

Remacle, J.F., Chevaugeon, N., Marchandise, E. and Geuzaine, C. (2007), “Efficient visualization of high-order finite elements”, International Journal for Numerical Methods in Engineering, Vol. 69 No. 4, pp. 750-771.

Shahbazi, K. (2005), “An explicit expression for the penalty parameter of the interior penalty method”, Journal of Computational Physics, Vol. 205 No. 2, pp. 401-407.

Smolianski, A. (2005), “Finite-element/level-set/operator-splitting (FELSOS) approach for computing two-fluid unsteady flows with free moving interfaces”, International Journal for Numerical Methods in Fluids, Vol. 48 No. 3, pp. 231-269.

SPHERIC (2016) “SPHERIC: 3D schematic dam break and evolution of the free surface”, available at: (accessed 5 July 2016).

Stapór, P. (2016), “A two-dimensional simulation of solidification processes in materials with thermo-dependent properties using xfem”, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 26 No. 6, pp. 1661-1683.

Sussman, M. (2005), “A parallelized, adaptive algorithm for multiphase flows in general geometries”, Computers & Structures, Vol. 83 Nos 6/7, pp. 435-444.

Sussman, M. and Fatemi, E. (1999), “Efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow”, SIAM Journal on Scientific Computing, Vol. 20 No. 4, pp. 1165-1191.

Sussman, M. and Hussaini, M.Y. (2003), “A discontinuous spectral element method for the level set equation”, Journal of Scientific Computing, Vol. 19 Nos 1/3, pp. 479-500.

Sussman, M. and Puckett, E.G. (2000), “A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows”, Journal of Computational Physics, Vol. 162 No. 2, pp. 301-337.

Sussman, M., Smereka, P. and Osher, S. (1994), “A level set approach for computing solutions to incompressible two-phase flow”, Journal of Computational Physics, Vol. 114 No. 1, pp. 146-159.

Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H. and Welcome, M.L. (1999), “An adaptive level set approach for incompressible two-phase flows”, Journal of Computational Physics, Vol. 148 No. 1, pp. 81-124.

Tryggvason, G. (1988), “Numerical simulations of the Rayleigh-Taylor instability”, Journal of Computational Physics, Vol. 75 No. 2, pp. 253-282.

Unverdi, S.O. and Tryggvason, G. (1992), “A front-tracking method for viscous, incompressible, multi-fluid flows”, Journal of Computational Physics, Vol. 100 No. 1, pp. 25-37.

Wang, C.-Y., Teng, J.-T. and Huang, G.P.G. (2011), “Numerical simulation of sloshing motion inside a two dimensional rectangular tank by level set method”, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 21 No. 1, pp. 5-31.

Warburton, T. (2006), “An explicit construction of interpolation nodes on the simplex”, Journal of Engineering Mathematics, Vol. 56 No. 3, pp. 247-262.

Wang, Y., Simakhina, S. and Sussman, M. (2012), “A hybrid level set-volume constraint method for incompressible two-phase flow”, Journal of Computational Physics, Vol. 231 No. 19, pp. 6438-6471.

Wu, G.X., Taylor, R.E. and Greaves, D.M. (2001), “The effect of viscosity on the transient free-surface waves in a two-dimensional tank”, Journal of Engineering Mathematics, Vol. 40 No. 1, pp. 77-90.

Xie, Z., Pavlidis, D., Percival, J.R., Gomes, J.L.M.A., Pain, C.C. and Matar, O.K. (2014), “Adaptive unstructured mesh modelling of multiphase flows”, International Journal of Multiphase Flow, Vol. 67, pp. 104-110.

Yan, J. and Osher, S. (2011), “A local discontinuous galerkin method for directly solving Hamilton-Jacobi equations”, Journal of Computational Physics, Vol. 230 No. 1, pp. 232-244.


The first author thanks Dr Rajesh Gandham, Dr David Medina and Prof Jesse Chan for their helpful discussions and valuable comments.

Corresponding author

Ali Karakus can be contacted at: