# A high-order diffused-interface approach for two-phase compressible flow
simulations using a Discontinuous Galerkin framework

Niccolò Tonicello<EMAIL_ADDRESS>Matthias Ihme  Department of
Mechanical Engineering, Stanford University, Stanford, CA 94305, United States

###### Abstract

A diffused-interface approach based on the Allen-Cahn phase field equation is
developed within a high-order Discontinuous Galerkin framework. The interface
capturing technique is based on the balance between explicit diffusion and
sharpening terms in the phase field equation, where the former term involves
the computation of the local interface normal vectors. Due to the well-known
Gibbs phenomenon encountered in high-order discretisations of steep profiles
such as shocks and/or interfaces, the accurate evaluation of the normal vector
requires special consideration. To this end, a non-linear preconditioning
strategy is proposed in this work where an additional smooth level-set
function advected by the velocity field is used for the evaluation of the
normal vectors. It is shown that for appropriate choices of numerical fluxes
and parameters of the model, the phase field remains bounded without any need
for explicit regularisation. The proposed diffused-interface technique is
implemented within a five equation model for fully compressible two-phase
flows. In order to preserve isolated interfaces, a quasi-conservative
discretisation of the five equation model is employed. A series of numerical
experiments of increasing complexity are performed in order to assess the
accuracy and robustness of the developed methodology, including two-phase
flows involving viscous effects, gravitational forces, and surface tension.

###### keywords:

High-order methods, Discontinuous Galerkin, phase field, two-phase flows

††journal: Journal of Computational Physics

##  1 Introduction

The constant increase in computational power has made the use of Computational
Fluid Dynamics (CFD) an increasingly valuable tool in the design process of
complex industrial applications [1]. Commonly available commercial CFD
softwares are mostly based on low-order numerical schemes such as finite
volumes or finite differences methods. While these numerical methods are
particularly robust and reliable, they usually lack in accuracy for complex
applications. Along these lines, the development of innovative high-order
schemes, such as Discontinuous Galerkin [2, 3, 4, 5], Flux Reconstruction [6,
7, 8] and Spectral Difference [9, 10] methods, has experienced a significant
growth in the CFD community, representing a promising alternative for the next
generation of CFD commercial codes.

Whereas the use of high-order spectral element methods is becoming more common
in the simulation of compressible aerodynamics problems [11, 12, 13, 14, 15,
16, 17, 18, 19, 20], multiphase applications still largely rely on low-order
discretisations, in particular for fully compressible flows.

High-order numerical schemes, due to the low numerical dissipation and
dispersion errors [21, 22, 23, 24, 25, 26, 27, 28] can be beneficial in the
resolution of small-scale structures which are often encounter in the
simulation of multiphase flows. In contrast, because of the delicate nature of
such schemes, careful attention is needed in order to retain stability
properties.

One common approach in the simulation of two-phase flows relies on the concept
of volume fraction and its most popular formulation is commonly known as
Volume of Fluid (VOF) method [29, 30]. The volume fraction is an auxiliary
function that is advected by the velocity field which varies between zero and
one and represents the ratio of primary to secondary fluid at a given
computational grid point. The interface between two immiscible fluids is then
represented by a sharp variation of the volume fraction.

The phase field approach [31, 32, 33] follows a similar concept with respect
to VOF methodologies where a function bounded between zero and one is used to
distinguish the two phases. However, instead of relying on the explicit
reconstruction of the interface as it is custom in VOF methods, a balance
between a diffusion and an anti-diffusion term in the transport equation of
the phase field is used in order to maintain a sharp, but at the same time
sufficiently smooth, diffused profile of the interface.

Once a clear separation between the two phases is identified, it is possible
to introduce an additional set of equations that can be used to model each
phase, potentially, with considerably different transport properties and
equations of state. At the same time, localised forces or fluxes can be easily
imposed at the interface as it often happens in the modelling of surface
tension effects [34].

The objective of the present work is to develop a simple but, at the same
time, accurate and robust approach to deal with compressible two-phase flows
within the framework of the Discontinuous Galerkin method. The proposed
methodology is able to preserve important properties in the simulation of two-
phase flows such as boundedness of the volume fraction, accurate evaluation of
the interface normal vectors, low mass conservation errors and exact
resolution of contact discontinuities. From an implementational point of view,
the proposed approach does not need explicit limiting, interface
reconstruction or mass redistribution, but it is simply based on appropriate
definitions of numerical fluxes and additional equations. Finally, the high-
order spatial discretisation provides significant advantages from many
different points of view, from considerably reducing mass conservation errors,
to avoiding spurious deformations of the interface for long time integration.
Along the same lines, the generality of this method allows easy implementation
of different physical phenomena such as viscous, gravitational and surface
tension forces. In other words, the novelty of the present method resides in
its generality and flexibility: with only minor, targeted modifications on the
central kernels of the numerical scheme, it is possible to preserve a large
amount of desirable features in two-phase flows simulations, adding, at the
same time, the benefits given by the high-order spatial discretisation.

The paper is structured as follows. In section 2, the five equation model for
two immiscible compressible fluids is introduced, including all the additional
terms and equations associated to the specific interface capturing approach
employed in the present work. In section 3, the Discontinuous Galerkin method
is first outlined in its general formulation for conservation laws. Its
implementation for the specific case of the five equation model is then
presented, including additional details on the interface capturing technique,
the specific choices of numerical fluxes and their related properties.
Subsequently, section 4 is dedicated to the numerical results in which both
kinematics tests and fully coupled two-phase flows for increasing levels of
complexity are considered in order to validate the present implementation.
Within the benchmark cases herein considered, extensive studies on different
types of elements, different orders of approximation and grid-convergence
investigations are carried out in order to assess the robustness and accuracy
of the proposed methodology for a wide range of different problems in two-
phase flows. Finally, in section 5 the key conclusions of this work are
discussed.

##  2 Five equation model

The system of conservation laws considered in this work aims at modelling two-
phase compressible flows including viscous effects, gravitational forces and
surface tension. The five equation model [35] was chosen as framework that can
be used to simulate many different conditions in multiphase flows. The present
formulation slightly differs from the original model by considering an
additional equation for an advected function $\psi$, which is used to compute
the interface normal vectors in a similar fashion with respect to the work by
Al-Salami et al. [36]. It is worthwhile mentioning that the method proposed in
[36] considered a weakly compressible formulation. In this work, the same
approach was generalised to the fully compressible five equation model.
Furthermore, their choice of parameters in the conservative phase field
equation did not satisfy boundedness of the phase field variable, leading to
non-negligible mass conservation errors. As it will be shown in the numerical
results section, these errors are significantly smaller with the present
formulation.

The full system, including contributions for considering viscous,
gravitational and surface tension effects reads:

$\displaystyle\frac{\partial\phi_{1}}{\partial
t}+\textbf{u}\cdot\nabla\phi_{1}$ $\displaystyle=\nabla\cdot\textbf{a}_{1},$
(1) $\displaystyle\frac{\partial\psi}{\partial t}+\textbf{u}\cdot\nabla\psi$
$\displaystyle=0,$ (2) $\displaystyle\frac{\partial\rho_{1}\phi_{1}}{\partial
t}+\nabla\cdot(\rho_{1}\phi_{1}\textbf{u})$
$\displaystyle=\nabla\cdot\textbf{R}_{1},$ (3)
$\displaystyle\frac{\partial\rho_{2}\phi_{2}}{\partial
t}+\nabla\cdot(\rho_{2}\phi_{2}\textbf{u})$
$\displaystyle=\nabla\cdot\textbf{R}_{2},$ (4)
$\displaystyle\frac{\partial\rho\textbf{u}}{\partial
t}+\nabla\cdot(\rho\textbf{u}\otimes\textbf{u}+P\mathbb{I})$
$\displaystyle=\nabla\cdot(\textbf{f}\otimes\textbf{u})+\nabla\cdot\boldsymbol{\tau}+\sigma\kappa\widehat{\textbf{n}}\delta_{\Gamma}+\rho\textbf{g},$
(5) $\displaystyle\frac{\partial\rho E}{\partial t}+\nabla\cdot((\rho
E+P)\textbf{u})$
$\displaystyle=\nabla\cdot(\textbf{f}k)+\sum_{l=1}^{2}\nabla\cdot(\rho_{l}H_{l}\textbf{a}_{l})+\nabla\cdot(\boldsymbol{\tau}\cdot\textbf{u})$
(6)
$\displaystyle+\sigma\kappa\textbf{u}\cdot\widehat{\textbf{n}}\delta_{\Gamma}+\rho\textbf{g}\cdot\textbf{u},$
(7)

where $\phi_{l}$ is the phase field associated to the $l$-th phase, $\rho_{l}$
is the density of each phase, $\psi$ is the additional level-set function,
$\rho\textbf{u}$ is the total momentum, $P$ is the pressure and $\rho E=\rho
e+\frac{1}{2}\rho||\textbf{u}||^{2}$ is the total energy.

In addition,

$\textbf{a}_{l}=\Gamma(\epsilon\nabla\phi_{l}-\phi_{l}(1-\phi_{l})\widehat{\textbf{n}}_{l}),\quad\textbf{R}_{l}=\rho_{l}^{(0)}\textbf{a}_{l},\quad\textbf{f}=\sum_{l=1}^{2}\textbf{R}_{l},\quad
k=\frac{1}{2}||\textbf{u}||^{2},$ (8)

and $H_{l}$ is the specific enthalpy of the $l$-th phase. In terms of modeling
of buoyancy, viscous stresses and surface tension:
$\boldsymbol{\tau}=2\mu(\mathbb{S}-1/3(\nabla\cdot\textbf{u})\mathbb{I})$ is
the viscous stress tensor, with $\mu$ the dynamic viscosity of the mixture
evaluated as $\mu=\mu_{1}\phi_{1}+\mu_{2}\phi_{2}$,
$\mathbb{S}=(\nabla\textbf{u}+\nabla\textbf{u}^{\intercal})/2$ is the strain-
rate tensor, g is the gravitational acceleration, $\sigma$ is the surface
tension coefficient, $\kappa=-\nabla\cdot\widehat{\textbf{n}}$ is the
curvature of the interface and $\delta_{\Gamma}=||\nabla\phi_{1}||$ is an
approximate delta function around the interface.

The system is then closed by relating internal energy with the pressure field
using an Equation of State (EOS). A classical choice is the Stiffened-Gas EOS
[37]:

$P=\frac{\rho
e-\bigg{(}\frac{\gamma_{1}P^{\infty}_{1}}{\gamma_{1}-1}\phi_{1}+\frac{\gamma_{2}P^{\infty}_{2}}{\gamma_{2}-1}\phi_{2}\bigg{)}}{\bigg{(}\frac{\phi_{1}}{\gamma_{1}-1}+\frac{\phi_{2}}{\gamma_{2}-1}\bigg{)}},$
(9)

where $\gamma_{l}$ and $P^{\infty}_{l}$ are the parameters of the EOS. From
the stiffened-gas equation of state it is possible to write the speed of sound
and specific enthalpy of each phases as

$c_{l}=\sqrt{\gamma_{l}\bigg{(}\frac{P+P^{\infty}_{l}}{\rho_{l}}\bigg{)}}\quad\mathrm{and}\quad
H_{l}=\frac{(P+P^{\infty}_{l})\gamma_{l}}{\rho_{l}(\gamma_{l}-1)}\quad\mathrm{for}\quad
l=1,2.$ (10)

Finally, for completeness, the following mixture relations apply:

$\displaystyle\phi_{2}=$ $\displaystyle 1-\phi_{1},$ (11) $\displaystyle\rho=$
$\displaystyle\rho_{1}\phi_{1}+\rho_{2}\phi_{2},$ (12)
$\displaystyle\frac{1}{\gamma-1}=$
$\displaystyle\phi_{1}\frac{1}{\gamma_{1}-1}+\phi_{2}\frac{1}{\gamma_{2}-1},$
(13) $\displaystyle P^{\infty}\frac{\gamma}{\gamma-1}=$
$\displaystyle\phi_{1}\frac{\gamma_{1}P_{1}^{\infty}}{\gamma_{1}-1}+\phi_{2}\frac{\gamma_{2}P_{2}^{\infty}}{\gamma_{2}-1}.$
(14)

In order to preserve the _Interface Equilibrium Condition_ (IEC), the quasi-
conservative formulation proposed by Cheng et al. [38] is herein employed. In
[38], a DG discretisation was coupled with an explicit limiting technique to
avoid oscillations of the phase field in proximity of the interface. In the
present work, instead, the sharpening/diffusion balance first proposed by Chiu
& Lin [33] was used. Even though similar results can be achieved with both
approaches, the simplicity of implementation is surely one of the main
advantages of the present strategy. The same technique, in fact, can be
adapted to any numerical scheme for appropriate choices of parameters and
numerical fluxes.

The present formulation was implemented within the opensource code _Quail_
[39].

##  3 Discontinuous Galerkin discretisation

In this section the Discontinuous Galerkin (DG) method will be briefly
introduced by considering a general system of conservation laws. The same
strategy is then applied to the discretisation of the phase field equation and
subsequently to the compressible five equation model for two-phase flows.

A general set of conservation laws can be written in a compact form as

$\frac{\partial{\textbf{w}}}{\partial{t}}+\nabla\cdot\textbf{F}(\textbf{w},\nabla\textbf{w})=\textbf{S}(\textbf{w},\nabla\textbf{w}).$
(15)

Before applying the DG approach to equation 15, we first introduce the
computational domain $\Omega$ as the partition of $N_{e}$ non-overlapping
discrete elements such that $\Omega=\cup_{n=1}^{N_{e}}\Omega_{n}$. Let us also
denote the boundary of the $n$-th element as $\partial\Omega_{n}$.

Secondly, we introduce the space of test functions in which the numerical
solution will be seeked into. In particular, a classical choice for DG scheme
consist in the functional space:

$\mathcal{V}:=\\{\varphi\in
L^{2}(\Omega):\varphi\arrowvert_{\Omega_{n}}\in\mathbb{P}_{\mathrm{p}}(\Omega_{n}),\forall\Omega_{n}\\},$
(16)

where $\mathbb{P}_{\mathrm{p}}$ is the space of piecewise continuous
polynomials of order not greater than p on $\Omega_{n}$ and

$L^{2}(\Omega)=\bigg{\\{}\varphi:\Omega_{n}\rightarrow\mathbb{R}\bigg{\arrowvert}\int_{\Omega_{n}}|\varphi(\textbf{x})|^{2}d\Omega\leq\infty\bigg{\\}}$
(17)

with $\mathbb{R}$ being the space of the real numbers. Different bases can be
used to define the functional space $\mathcal{V}$, both modal or nodal [2]. In
this work, classical Lagrange polynomials are used within a modal DG
framework.

The approximation of the global solution $\textbf{w}^{\delta}$ can be defined
as

$\textbf{w}^{\delta}=\oplus_{n=1}^{N_{e}}\textbf{w}^{\delta}_{n},$ (18)

where $\textbf{w}^{\delta}_{n}$ is the local discrete solution:

$\textbf{w}^{\delta}_{n}=\sum_{j=0}^{\mathrm{p}}\textbf{w}_{j}(t)\varphi_{j}(\textbf{x}).$
(19)

The local formulation of the DG method requires $\textbf{w}^{\delta}_{n}$ to
satisfy

$\int_{\Omega_{n}}\varphi_{i}\frac{\partial{\textbf{w}^{\delta}_{n}}}{\partial{t}}d\Omega+\int_{\Omega_{n}}\varphi_{i}\nabla\cdot\textbf{F}(\textbf{w}^{\delta}_{n},\nabla\textbf{w}^{\delta}_{n})d\Omega=\int_{\Omega_{n}}\varphi_{i}\textbf{S}(\textbf{w}_{n}^{\delta},\nabla\textbf{w}_{n}^{\delta})d\Omega\quad\forall\varphi_{i}\in\mathcal{V}.$
(20)

The second term on the left-hand side of equation 20 is the flux term. Upon
performing integration by parts, this term can be expressed as

$\int_{\Omega_{n}}\varphi_{i}\nabla\cdot\textbf{F}(\textbf{w}^{\delta}_{n},\nabla\textbf{w}^{\delta}_{n})d\Omega=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{F}(\textbf{w}^{\delta}_{n},\nabla\textbf{w}^{\delta}_{n})d\Omega+\oint_{\partial\Omega_{n}}\varphi_{i}\widehat{\textbf{F}}(\textbf{w}^{\delta,+}_{n},\textbf{w}^{\delta,-}_{n},\nabla\textbf{w}^{\delta,+}_{n},\nabla\textbf{w}^{\delta,-}_{n},\widehat{\textbf{m}})dS,$
(21)

where $\widehat{\textbf{m}}$ is the outward-pointing unit normal vector,
$(\cdot)^{+}$ and $(\cdot)^{-}$ denote the right and left state with respect
to the element’s interface $\partial\Omega_{n}$ and $\widehat{\textbf{F}}$ is
the numerical flux.

After exploiting the form of $\textbf{w}^{\delta}_{n}$ (equation 19), the
local discrete weak formulation reads:

$\sum_{j=0}^{\mathrm{p}}\frac{\mathrm{d}{\textbf{w}_{j}}}{\mathrm{d}{\mathrm{t}}}\mathrm{M}_{ij}=\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{F}(\textbf{w}^{\delta}_{n},\nabla\textbf{w}^{\delta}_{n})d\Omega-\oint_{\partial\Omega_{n}}\varphi_{i}\widehat{\textbf{F}}(\textbf{w}^{\delta,+}_{n},\textbf{w}^{\delta,-}_{n},\nabla\textbf{w}^{\delta,+}_{n},\nabla\textbf{w}^{\delta,-}_{n},\widehat{\textbf{m}})dS+\int_{\Omega_{n}}\varphi_{i}\textbf{S}(\textbf{w}_{n}^{\delta},\nabla\textbf{w}_{n}^{\delta},\textbf{x})d\Omega,$
(22)

where $\mathrm{M}_{ij}=\int_{\Omega_{n}}\varphi_{i}\varphi_{j}d\Omega$
represents the $(i,j)$-th entry of the element-local mass matrix. Both volume
and surface integrals appearing in equation 22 can be either computed
analytically or using appropriate quadrature rules. Equation 22 can then be
discretised in time with explicit or implicit schemes.

The only component of the spatial discretisation which is intrinsically
dependent on the specific conservation law is the definition of the numerical
fluxes, which needs to take into account the eigen-structure of the hyperbolic
system. Lax-Friedrichs and Symmetric Interior Penalty [40] fluxes will be
considered as they offer good flexibility for a wide range of different
conservation laws.

Now that the general framework of the DG method is introduced, we will proceed
in presenting the specific numerical strategy used in the discretisation of
the five equation model. First, the numerical treatment of the phase field and
level-set equations will be considered in order to highlight the key aspects
of the interface capturing technique only. The subsequent coupling with the
five equation model will naturally follow with appropriate choices of
numerical fluxes.

The first two equations of the five equation model introduced in the previous
section read:

$\displaystyle\frac{\partial\phi_{1}}{\partial
t}+\textbf{u}\cdot\nabla\phi_{1}$ $\displaystyle=\nabla\cdot\textbf{a}_{1},$
(23) $\displaystyle\frac{\partial\psi}{\partial t}+\textbf{u}\cdot\nabla\psi$
$\displaystyle=0,$ (24)

with
$\textbf{a}_{1}=\Gamma(\epsilon\nabla\phi_{1}-\phi_{1}(1-\phi_{1})\widehat{\textbf{n}}_{1})$
and $\widehat{\textbf{n}}_{1}=\nabla\psi/||\nabla\psi||$.

For a prescribed velocity field, this system is well-posed and it can be
numerically resolved as it is common practice in kinematic tests for the
validation of interface capturing techniques. In the following discussion, in
the phase field equation, the subscript of $\phi_{1}$ will be dropped since
only one phase field needs to be resolved (the phase field for the secondary
fluid will be directly evaluated as $\phi_{2}=1-\phi_{1}$).

The local weak formulation of the classic Discontinuous Galerkin method for
this set of equations reads:

$\displaystyle\int_{\Omega_{n}}\frac{\partial\phi_{n}}{\partial
t}\varphi_{i}d\Omega+\int_{\Omega_{n}}(\textbf{u}\cdot\nabla\phi_{n})\varphi_{i}d\Omega$
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{a}_{l}d\Omega+\oint_{\partial\Omega_{n}}\widehat{\textbf{a}}_{l}\cdot\widehat{\textbf{m}}\varphi_{i}dS,$
(25) $\displaystyle\int_{\Omega_{n}}\frac{\partial\psi_{n}}{\partial
t}\varphi_{i}d\Omega+\int_{\Omega_{n}}(\textbf{u}\cdot\nabla\psi_{n})\varphi_{i}d\Omega$
$\displaystyle=0,$ (26)

where the vector $\widehat{\textbf{a}}_{l}$ is the numerical flux at the
interface between neighbouring elements which depends on the left and right
state of the conservative variables and their gradients (namely, $\phi^{\pm}$,
$\psi^{\pm}$, $\nabla\phi^{\pm}$, $\nabla\psi^{\pm}$). A symmetric interior
penalty flux was chosen in the present work to evaluate this term. The
transport term, instead, is treated as a source term. The system written as it
is can be numerically resolved for a prescribed, analytical velocity field.
Instead, when considering the fully coupled system where the velocity field is
itself an unknown of the problem and modelling the advection term as a source
term is needed in order to fulfil the interface equilibrium condition (_i.e._
exact resolution of contact discontinuities). This property is of fundamental
importance for compressible flows in order to avoid pressure and velocity
oscillations in proximity of the interface. These oscillations can ultimately
lead to global instability of the numerical scheme and thus the IEC represents
a key aspect in the success of the simulation.

Integrating by parts twice, the previous system can be written as:

$\int_{\Omega_{n}}\frac{\partial\phi_{n}}{\partial
t}\varphi_{i}d\Omega+\int_{\Omega_{n}}(\textbf{u}\cdot\nabla\phi_{n})\varphi_{i}d\Omega-\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widetilde{(\textbf{u}\phi_{n})}\varphi_{i}dS+\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widehat{(\textbf{u}\phi_{n})}\varphi_{i}dS=$
$=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{a}_{l}d\Omega+\oint_{\partial\Omega_{n}}\widehat{\textbf{a}}_{l}\cdot\widehat{\textbf{m}}\varphi_{i}dS$
$\int_{\Omega_{n}}\frac{\partial\psi_{n}}{\partial
t}\varphi_{i}d\Omega+\int_{\Omega_{n}}(\textbf{u}\cdot\nabla\psi_{n})\varphi_{i}d\Omega-\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widetilde{(\textbf{u}\psi_{n})}\varphi_{i}dS+\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widehat{(\textbf{u}\psi_{n})}\varphi_{i}dS=0$

where two additional numerical fluxes for the transport term need to be
defined (namely, $\widehat{(\cdot)}$ and $\widetilde{(\cdot)}$). The same
expressions proposed by [38] were considered in this work. In particular,

$\widehat{\textbf{m}}\cdot\widetilde{(\textbf{u}\phi)}=\begin{cases}\widehat{\textbf{m}}\cdot\textbf{u}_{n}^{-}\phi^{-}_{n}\enspace\mathrm{in}\enspace\Omega^{-}_{n}\\\
\widehat{\textbf{m}}\cdot\textbf{u}_{n}^{+}\phi^{+}_{n}\enspace\mathrm{in}\enspace\Omega^{+}_{n}\\\
\end{cases}\enspace\mathrm{and}\enspace\widehat{\textbf{m}}\cdot\widehat{(\textbf{u}\phi)}=\begin{cases}\widehat{\textbf{m}}\cdot\textbf{u}^{-}\\{\\!\\!\\{\phi\\}\\!\\!\\}-c[\\![\phi]\\!]\enspace\mathrm{in}\enspace\Omega^{-}_{n}\\\
\widehat{\textbf{m}}\cdot\textbf{u}^{+}\\{\\!\\!\\{\phi\\}\\!\\!\\}-c[\\![\phi]\\!]\enspace\mathrm{in}\enspace\Omega^{+}_{n}\\\
\end{cases}$ (27)

where $c$ represents the maximum characteristic speed of the system (_i.e._
the speed of sound for the five equation model),

$\\{\\!\\!\\{\phi\\}\\!\\!\\}=\frac{\phi^{+}+\phi^{-}}{2}\quad\mathrm{and}\quad[\\![\phi]\\!]=\frac{\phi^{+}-\phi^{-}}{2}.$

It can be proven that this choice, in addition to the quasi-conservative form
of the advection term, leads to an exact conservation of contact
discontinuities if coupled with standard Lax-Friedrichs numerical flux for the
remaining conservation laws in the five equation model (_i.e._ conservation of
mass of the $l$-th phase, total momentum and total energy). For consistency,
the same fluxes are employed also in the level-set equation even if this
choice is not fundamental in the fulfilment of the interface equilibrium
condition.

The local weak formulation for each equation of the full system can be written
for

the Allen-Cahn equation as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial\phi_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle+\int_{\Omega_{n}}(\textbf{u}_{n}\cdot\nabla\phi_{n})\varphi_{i}d\Omega-\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widetilde{(\textbf{u}_{n}\phi_{n})}\varphi_{i}dS+\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widehat{(\textbf{u}_{n}\phi_{n})}\varphi_{i}dS=$
(28)
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{a}_{1}d\Omega+\oint_{\partial\Omega_{n}}\widehat{\textbf{a}}_{1}\cdot\widehat{\textbf{m}}\varphi_{i}dS;$
(29)

the level-set equation as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial\psi_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle+\int_{\Omega_{n}}(\textbf{u}_{n}\cdot\nabla\psi_{n})\varphi_{i}d\Omega-\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widetilde{(\textbf{u}_{n}\psi_{n})}\varphi_{i}dS+\int_{\partial\Omega_{n}}\widehat{\textbf{m}}\cdot\widehat{(\textbf{u}_{n}\psi_{n})}\varphi_{i}dS=0;$
(30)

the mass conservation of phase $1$ as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial(\rho_{1}\phi_{1})_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle-\int_{\Omega_{n}}(\rho_{1}\phi_{1})_{n}\textbf{u}_{n}\cdot\nabla\varphi_{i}d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\rho_{1}\phi_{1})_{n}\textbf{u}_{n}}\cdot\widehat{\textbf{m}}\varphi_{i}dS=$
(31)
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{R}_{1}d\Omega+\oint_{\partial\Omega_{n}}\widehat{\textbf{R}}_{1}\cdot\widehat{\textbf{m}}\varphi_{i}dS;$
(32)

the mass conservation of phase $2$ as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial(\rho_{2}\phi_{2})_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle-\int_{\Omega_{n}}(\rho_{2}\phi_{2})_{n}\textbf{u}_{n}\cdot\nabla\varphi_{i}d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\rho_{2}\phi_{2})_{n}\textbf{u}_{n}}\cdot\widehat{\textbf{m}}\varphi_{i}dS=$
(33)
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot\textbf{R}_{2}d\Omega+\oint_{\partial\Omega_{n}}\widehat{\textbf{R}}_{2}\cdot\widehat{\textbf{m}}\varphi_{i}dS;$
(34)

the total momentum conservation as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial(\rho\textbf{u})_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle-\int_{\Omega_{n}}(\rho_{n}\textbf{u}_{n}\otimes\textbf{u}_{n}+P_{n}\mathbb{I})\cdot\nabla\varphi_{i}d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\rho_{n}\textbf{u}_{n}\otimes\textbf{u}_{n}+P_{n}\mathbb{I})}\cdot\widehat{\textbf{m}}\varphi_{i}dS=$
(35)
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot(\textbf{f}\otimes\textbf{u}_{n})d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\textbf{f}\otimes\textbf{u}_{n})}\cdot\widehat{\textbf{m}}\varphi_{i}dS;$
(36)

and the total energy conservation as:

$\displaystyle\int_{\Omega_{n}}\frac{\partial(\rho E)_{n}}{\partial
t}\varphi_{i}d\Omega$
$\displaystyle-\int_{\Omega_{n}}(\rho_{n}E_{n}+P_{n})\textbf{u}_{n}\cdot\nabla\varphi_{i}d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\rho_{n}E_{n}+P_{n})\textbf{u}_{n}}\cdot\widehat{\textbf{m}}\varphi_{i}dS=$
(37)
$\displaystyle=-\int_{\Omega_{n}}\nabla\varphi_{i}\cdot(\textbf{f}k+\sum_{l=1}^{2}\rho_{l}H_{l}\textbf{a}_{l})d\Omega+\oint_{\partial\Omega_{n}}\widehat{(\textbf{f}k+\sum_{l=1}^{2}\rho_{l}H_{l}\textbf{a}_{l})}\cdot\widehat{\textbf{m}}\varphi_{i}dS.$
(38)

The numerical fluxes related to the interface capturing technique (_i.e._
terms involving $\textbf{a}_{l}$) are discretised using the symmetric interior
penalty method whereas the remaining numerical convective fluxes are
discretised using the Lax-Friedrichs flux.

Finally, it is relevant to present one additional detail related to the link
between the phase field and the additional level-set function $\psi$ and their
respective dynamics. The function $\psi$ can be interpreted as a smoother
version of the phase field which is used to compute the interface normal
vectors.

However, it is easy to see that in the previous equations, for a given
velocity field, the coupling is only one-way: the level-set equation will
affect the phase field evolution due to the computation of the normal vectors
but there is no information flowing from the phase field equation back to the
dynamics of the level-set. Consequently, after sufficiently long time, the
interface defined by the two functions will start to drift apart, leading to
completely erroneous results. Consequently, the function $\psi$ needs to be
periodically re-initialised. In this work we followed the strategy proposed by
[36] where every $2000$ time steps, the following re-initialisation equation

$\frac{\partial\psi}{\partial\tau}+sgn(\psi_{0})(1-||\nabla\psi||)=\nabla\cdot(\nu_{h}\nabla\psi),$
(39)

is solved in pseudo-time from the initial condition
$\psi_{0}=\phi-\frac{1}{2}$. The quantity

$sgn(\psi_{0})=\tanh\bigg{(}\frac{\psi_{0}}{2\epsilon||\nabla\psi_{0}||}\bigg{)}$

is a smeared sign function and $\nu_{h}$ is a vanishing viscosity proportional
to the grid size. Since a smooth profile for the level-set function is only
needed inside a narrow bandwidth localised in proximity of the interface, the
re-initialisation equation is integrated only for a limited number of pseudo-
time steps. It was found that applying approximately $50$ iterations every
$2000$ physical time-steps was at the same time able to  (i) provide a
sufficiently smooth profile of the interface,  (ii) avoid decoupling between
phase field and level-set and  (iii) keep the computational cost of the re-
initialisation process relatively low.  In similarity with the phase field
equation, the numerical diffusive flux in the re-initialisation equation is
discretised using the symmetric interior penalty approach whereas the smeared
sign hyperbolic term is treated as a simple source term. It is worthwhile
mentioning that even if a re-initialisation step is herein considered, the
mass conservation errors associated to this process are identically equal to
zero since the level-set function is only used to compute the interface normal
vectors and not to track the interface (as it is normally done in classical
level-set methods).

It is important to underline that the present methodology improves the
previously cited works from several point of views, increasing either the
accuracy, the simplicity of implementation or the generality of each model.
The quasi-conservative approach presented in [38] is used to preserve contact
discontinuities in the five equation model and it is herein further improved
by avoiding explicit limiting of the volume fraction employing, instead,
explicit sharpening/diffusive terms of the Allen-Cahn equation. In this way
the proposed model is characterised by a significant simplicity of
implementation and flexibility. For example, considering complex unstructured
grids, the proposed methodology does not need any ad-hoc limiting which might
become considerably expensive and complex to implement. The present approach
also allows the use of much higher orders of approximation with respect to the
ones used in [38]. The advantages of using high orders of approximation will
be further discussed in the results section. In the same section, the proposed
methodology will also be tested on triangular elements, highlighting the
robustness of this approach for unstructured grids. Finally, the model herein
presented extends the interface capturing technique proposed by [36] to a
fully compressible formulation and it was also possible, by appropriate
choices of the parameters, to significantly reduce mass conservation errors of
the original technique.

##  4 Results

###  4.1 Kinematic tests

Regarding the validation of the present methodology, the proposed numerical
test cases will be sorted in kinematic tests and two-phase flow tests. For the
former case, the velocity field is represented by a given analytical function
which will be used in the phase field and level-set equations. In this way, it
is possible to assess the performance of the interface capturing technique in
a segregated fashion. The latter set of tests considers the full five equation
system which is used to describe the flow of compressible two-phase flows.

####  4.1.1 Linear-advection of a droplet

The first test case herein considered is a simple advection of a circular
droplet with a prescribed interface thickness. The goal of this specific test
is to verify that the additional diffusive/sharpening terms, for a given
hyperbolic tangent initial condition, do not deteriorate the order of accuracy
of the underlying numerical scheme. The hyperbolic tangent profile, in fact,
is supposed to maintain its shape and simply be transported by the constant
advection velocity.

Consequently, a spherical droplet of radius $R=0.15$, represented by the phase
field

$\phi(\textbf{x},0)=\frac{1}{2}\bigg{(}1+\tanh\bigg{(}\frac{r-R}{2\epsilon}\bigg{)}\bigg{)}\quad\mathrm{with}\quad
r=||\textbf{x}||,$ (40)

is advected by the constant advection velocity $\textbf{u}=(1,0)$ in a
$[0,1]^{2}$ periodic square for a full period. The parameter $\epsilon$ is set
to be a constant equal to the coarsest mesh resolution (_i.e._
$\epsilon=0.05$).

Similarly, the level-set, signed-distance function is initialised as:

$\psi(\textbf{x},0)=-r.$ (41)

The $L^{1}$ error of the phase field as a function of the number of elements
and different orders of approximation is shown in figure 1. It can be seen
that the optimal order of convergence of the spatial discretisation is
correctly recovered under mesh refinement.

Figure 1: $L^{1}$-error against number of elements along each direction for
different orders of approximation. Optimal orders of convergence are shown by
the dashed grey lines.

####  4.1.2 Rider-Kothe vortex

Another standard test case for interface capturing techniques consists in the
deformation of a circular bubble by a shear flow [41]. In particular, a
droplet of radius $R=0.15$ centered at $(0.0,0.25)$ in a $[-0.5,0.5]^{2}$
periodic domain is advected by the following divergence-free velocity field:

$\displaystyle u_{1}=$ $\displaystyle\sin^{2}(\pi x_{1})\sin(2\pi
x_{2})\cos\bigg{(}\frac{\pi t}{T}\bigg{)}$ (42) $\displaystyle u_{2}=$
$\displaystyle-\sin(2\pi x_{1})\sin^{2}(\pi x_{2})\cos\bigg{(}\frac{\pi
t}{T}\bigg{)},$ (43)

where $T$ denotes the characteristic period of the shear flow. The classical
value of $T=4$ is chosen for this particular test. Under the action of such a
velocity field, in the first half of the period, the initial circular bubble
is strongly deformed into a thin filament. After velocity reversal, the
filament is stretched back to the initial condition.

For this problem, we considered a polynomial discretisation of order $1$ and
$3$ within the DG framework on a series of increasingly refined meshes. In
particular, simulations involving $64^{2}$, $128^{2}$ and $256^{2}$ degrees of
freedom were considered. Notice that in spectral element methods the total
number of degrees of freedom is jointly defined by the number of elements and
by the polynomial order of approximation. Consequently, for example,
considering a $3^{\mathrm{rd}}$ polynomial order approximation on a $16\times
16$ grid, the nomenclature will read: $16\times 16$p$3$ (_i.e._ $16$ elements
times $16$ elements for a polynomial approximation of degree $3$). The total
number of degrees of freedom is defined as $(16\times 4)^{2}=64^{2}$. In
comparisons of simulations involving different orders of approximation, a
common choice is to match the total number of degrees of freedom (for example,
$16\times 16$p$3$ and $32\times 32$p$1$ simulations).

Time integration was performed using a classical $4^{\mathrm{th}}$ order
Runge-Kutta scheme. Also, in order to test the versatility of the proposed
approach, both quadrilateral and triangular meshes were used. Examples of the
meshes used in this work are shown in figure 2.

(a) $16\times 16$ quadrilateral mesh (mesh $a$).

(b) $16\times 16$ Triangular mesh (mesh $b$).

Figure 2: Examples of computational meshes used in the present work.

In figure 3, the interfacial profiles at half period and after a full period
are shown for different levels of mesh refinement using a $3^{\mathrm{rd}}$
order polynomial approximation based on square elements.

(a) Interface location at $t=T/2$.

(b) Interface location at $t=T$.

Figure 3: Iso-contour $\phi=0.5$ after half period (left) and after one full
period (right) for the $16\times 16$p$3$, $32\times 32$p$3$ and $64\times
64$p$3$ simulations. The solid black line represents the initial location of
the interface.

Two close-up views of figures 3(a) and 3(b) are shown in figure 4 to better
appreciate the convergence of the proposed approach.

In figure 4(a), we can observe that the resolution of the thin filament
improves under mesh refinement. In figure 4(b), instead, the recovery of the
initial condition can be better appreciated.

(a) Interface location at $t=T/2$.

(b) Interface location at $t=T$.

Figure 4: Close-up view of the interface iso-contour at half period near
breakup (left) and after one full period (right) for the $16\times 16$p$3$,
$32\times 32$p$3$ and $64\times 64$p$3$ simulations.

In order to further investigate the proposed numerical strategy, lower order
simulations using the same approach have been performed, keeping fixed the
total number of degrees of freedom. In figure 5 the final location of the
interface is shown for two different polynomial orders ($1$ and $3$) for two
different resolutions ($128^{2}$ and $256^{2}$ DoF).

It can be noticed that the higher order approximation provides more accurate
results. For both resolutions, in fact, the p$3$ simulation outperforms the
p$1$ computation in matching the exact initial solution. In particular, the
upper part of the droplet is better approximated. In fact, this region is
particularly stretched in the first half-period leading to significant under-
resolution.

(a) $128^{2}$ DoF

(b) $256^{2}$ DoF

Figure 5: Iso-contour $\phi=0.5$ after one full period for $128^{2}$ (left)
and $256^{2}$ DoF (right) with $1^{\mathrm{st}}$ and $3^{\mathrm{rd}}$
polynomial order approximations. The solid black line represents the initial
location of the interface.

In order to further examine how under-resolution behaves, the location of the
interface at the half-period is shown in figure 6 for the two different
polynomial orders. Different zones of the interface are highlighted to better
appreciate the differences between the two simulations. It can be seen that in
many locations the higher order approximation is generally smoother and it
better matches the reference solution. In this case, the reference solution is
simply a more refined simulation.

These observations further highlight the significant improvement in
considering high-order discretisations for the proposed phase field equation.
In particular, it is the believe of the authors that low-order approximations
not only introduce excessive amounts of numerical dissipation but also provide
inaccurate evaluation of the normal vectors, therefore deteriorating the
interface capturing technique. These assumptions will be further confirmed
through additional numerical experiments for different orders of
approximation.

Figure 6: Iso-contour $\phi=0.5$ after half period for $128^{2}$ DoF with
$1^{\mathrm{st}}$ ( ) and $3^{\mathrm{rd}}$ ( ) polynomial order
approximations on rectangular elements. Some specific regions of the the
interface are enlarge to highlight the difference between the two simulations.
The simulation for $256^{2}$ DoF is used as reference.

Secondly, in order to test the robustness and generality of the proposed
approach, different types of finite elements were considered. In particular,
in the following part of this section we compare quadrilateral and triangular
elements (see figure 2).

In figure 7, the location of the interface using quadrilateral and triangular
elements is shown for $128^{2}$ DoF ($3^{\mathrm{rd}}$ polynomial order). In
particular, the solution is shown at half-period on the left and after a full
revolution on the right. In figure 7(a) it is interesting to note that the two
type of meshes (rectangular and triangular) provide almost exactly the same
results in well-resolved regions of the domain, whereas noticeable differences
can be observed only in proximity of the tail of the droplet. After one full
revolution, the upper part of the droplet is slightly better approximated by
triangular element as it can be observed in figure 7(b). It is important to
highlight that such results further confirm the robustness and generality of
the present approach. Using different types of elements did not require any
ad-hoc modification of the underlying technique, which emphasises the
simplicity and flexibility of the proposed implementation. The capability of
handling triangular elements without deep modifications in the interface
capturing techniques is a very desirable feature for the simulation of more
complex configurations of engineering interest where unstructured meshes are
often required.

(a) Interface location at $t=T/2$.

(b) Interface location at $t=T$.

Figure 7: Iso-contour $\phi=0.5$ after half period (left) and one full period
(right) for the the $3^{\mathrm{rd}}$ polynomial order simulation with
$128^{2}$ DoF comparing quadrilateral and triangular elements. The solid black
line represents the reference solution (highly resolved simulation for $t=T/2$
and initial solution for $t=T$).

In a more quantitative way, the $L^{1}$-errors were computed and are listed in
table 1 where they are compared with simulations from the literature involving
the same number of degrees of freedom and similar interface capturing
techniques (_i.e._ based on the phase field approach). It can be observed that
the present results are in good agreement with previously published
simulations of the same test case. In the same table, it can also be noted
that the differences between $1^{\mathrm{st}}$ and $3^{\mathrm{rd}}$
polynomial order approximations observed in previous plots, can now be
quantified more clearly: the $L^{1}$-errors in the p$3$ simulations are always
smaller than the p$1$ counterparts (sometimes even twice/three times smaller).
Also, the gap between the two approximation orders tends to grow under mesh
refinement.

Table 1: $L^{1}$-errors for the Rider-Kothe test case compared with respect to simulations from the literature involving the same number of degrees of freedom. Method | $64^{2}$ | $128^{2}$ | $256^{2}$
---|---|---|---
Mirjalili et al. [42], mesh (a) | $1.96\times 10^{-2}$ | $5.43\times 10^{-3}$ | $1.25\times 10^{-3}$
Al-Salami et al. [36], mesh (a), p$3$ | $-$ | $5.50\times 10^{-3}$ | $1.05\times 10^{-3}$
Present mesh (a), p$3$ | $1.75\times 10^{-2}$ | $4.83\times 10^{-3}$ | $1.23\times 10^{-3}$
Present mesh (a), p$1$ | $2.22\times 10^{-2}$ | $8.71\times 10^{-3}$ | $3.31\times 10^{-3}$
Al-Salami et al. [36], mesh (b), p$3$ | $-$ | $5.51\times 10^{-3}$ | $1.25\times 10^{-3}$
Present mesh (b), p$3$ | $1.58\times 10^{-2}$ | $4.03\times 10^{-3}$ | $1.18\times 10^{-3}$
Present mesh (b), p$1$ | $1.99\times 10^{-2}$ | $8.23\times 10^{-3}$ | $2.65\times 10^{-3}$

Finally, another feature of great interest for interface capturing techniques
is mass conservation. Consequently, the relative mass conservation error was
evaluated over time as

$E_{m}=\frac{\int(\phi(\textbf{x},t)-\phi(\textbf{x},0))d\Omega}{\int\phi(\textbf{x},0)d\Omega}.$
(44)

The time evolution of the mass conservation error throughout the whole
simulation is shown in figure 8 for different mesh resolutions and different
orders of approximation.

Since the transport term in the phase field equation is discretised in a non-
conservative form, mass conservation errors are not identically zero.
Considering the p$1$ simulation, the mass conservation errors are relatively
small, ranging between $10^{-8}$ and $10^{-11}$ for different mesh
resolutions, as the grid size decreases. Even if these errors are already
considerably small for a first order polynomial approximation, it can be
observed that the mass conservation errors are even smaller for the p$3$
simulation, as they range between $10^{-10}$ for the coarsest mesh to almost
machine precision for the most refined simulation.

In the work by Al-Salami et al. [36], in the same case for $256^{2}$ DoF, mass
conservation errors were larger by $11$ orders of magnitude.

Figure 8: Mass conservation errors for different resolutions on mesh (a) using
a p$1$ (left) and p$3$ approximation (right) with $64^{2}$, $128^{2}$ and
$256^{2}$ DoF.

###  4.2 Two-phase flows

Following the analysis of the interface capturing technique, we will
considered a series of validation tests involving the simulation of two-phase
flows using the five equation model presented in the previous section.

####  4.2.1 Advection of a water droplet in air

A popular test case to validate the present implementation consists in the
advection of a water droplet in air. The one-dimensional form of the present
test case has been extensively used as a benchmark validation problem in the
multiphase flows community [43, 35, 44, 45].

The goal of this test is to assess the ability of the numerical scheme to
preserve the initial circular shape of the droplet for long time integration
and to quantify the robustness of the solver in dealing with high-density-
ratio interfaces.

Another goal of the water bubble advection case consists in verifying the
capability of the numerical discretisation to exactly preserve isolated
material interfaces. A well-balanced discretisation of the five equation model
should, in fact, maintain pressure and velocities exactly constant during time
integration. For this specific case viscous, gravitational and surface tension
effects are neglected.

The bubble of radius $R=25/89$ is located at the center of a $[0,1]^{2}$
periodic square. The material properties for the air medium for this test case
are $\gamma_{1}=1.4$, $\rho_{1}=1\times 10^{-3}$ and $P_{1}^{\infty}=0$,
whereas for the water medium the fluid properties are $\gamma_{2}=4.4$,
$\rho_{2}=1$ and $P^{\infty}_{2}=6\times 10^{3}$.

The initial conditions read:

$\textbf{u}=(5,5),\quad
P=1,\quad\phi_{1}=\frac{1}{2}\bigg{(}1+\tanh\bigg{(}\frac{r-R}{2\epsilon}\bigg{)}\bigg{)}\quad\mathrm{and}\quad\rho=\rho_{2}+(\rho_{1}-\rho_{2})\phi_{1}.$
(45)

The simulation is carried on different mesh resolutions for a total of five
advection periods. $3^{\mathrm{rd}}$ and $1^{\mathrm{st}}$ polynomial order DG
discretisations were considered and time integration was performed using a
$4^{\mathrm{th}}$ order Runge-Kutta scheme.

In figure 9, the phase field at the end of the simulations is shown for the
most refined simulations (equivalent to $96^{2}$ DoF). It can be noticed that
in the p$3$ simulation the droplet is qualitatively identical with respect to
the prescribed initial condition, indicating that the proposed methodology is
successfully able to preserve the initial shape of the water bubble even after
long time integration. It is worthwhile mentioning that a spurious
deformations of the interface and the grid directions is often encountered for
low-order approximations of the interface normals [46]. It is known that the
use of high-order methods in the advection of the phase field can mitigated
such numerical artefacts [46, 47]. The present simulation confirms this
behaviour, where the p$1$ discretisation shows a much more evident deformation
of the droplet after five advection periods.

These results further confirm the overall benefit of using high-order
discretisations already observed in the previous Rider-Kothe test case.

Figure 9: phase field isocontours after five advection periods for the most
refined simulation. Left, $24^{2}$p$3$; right, $48^{2}$p$1$. The interface is
represented by the solid black line. The initial condition is shown as a
dashed black line.

In order to have a better understanding of the overall behaviour of all the
relevant quantities, in figure 10 the phase field, density, pressure and
$x$-velocity are respectively plotted along the line $y=0$ for the p$3$ most
refined simulation ($24^{2}$p$3$). It can be observed that both the phase
field and the total density vary smoothly across the interface, preserving the
prescribed hyperbolic tangent profile without any spurious oscillation in
proximity of the interface.

Finally, we remark that pressure and $x$-velocity remain uniform during the
simulation as shown in figure 10, indicating that the proposed discretisation
is able to fulfil the interface equilibrium condition.

Figure 10: Top row: phase field (left) and density (right) along the line
$y=0$ after five advection periods. Bottom row: slice of pressure (left) and
$x$-velocity (right) along the line $y=0$ after five advection periods. Red,
numerical simulation; black, exact solution. Notice that the $y$-axis is
centered around the expected exact value and scaled by $10^{-6}$.

####  4.2.2 Rayleigh-Taylor instability

The Rayleigh–Taylor instability occurs when an interface between two fluids
with different densities experiences a pressure gradient opposing the density
gradient. The domain consists in a $[d\times 4d]$ with $d=1$. The interface is
initially defined as the curve $y(x)=2d+0.1d\cos{(2\pi x/d)}$. The Rayleigh-
Taylor instability is characterised by the Reynolds number
$\mathrm{Re}=(\rho_{1}d^{3/2}||\textbf{g}||^{1/2})/\mu$ and the Atwood number
$\mathrm{At}=(\rho_{1}-\rho_{2})/(\rho_{1}+\rho_{2})$ which are respectively
set to $3000$ and $0.5$. The top boundary is treated as a Riemann-invariant
boundary condition with zero velocity and constant pressure, the bottom
boundary is a no-slip wall whereas slip wall boundary conditions are
prescribed on the other lateral sides of the domain. Two different
quadrilateral grids were considered for this study involving $5000$ and
$20000$ degrees of freedom ($10\times 20$p$4$ and $20\times 80$p$4$,
respectively). Finally, a $4^{\mathrm{th}}$ order Runge-Kutta scheme was used
for time integration.

The present test is meaningful in considering a more complex physical set-up
where viscous and gravitational force drive the dynamics of the system.

A series of snapshots of the phase field are shown in figure 11 for the two
grids.

Figure 11: Snapshots of the phase field isocontours for different meshes and
times. Top row: isocontour of the phase field at $t^{*}=0.5$ (left) and
$t^{*}=1.5$ (right). Bottom row: isocontour of the phase field at $t^{*}=2.5$
(left) and $t^{*}=3.0$ (right). In each subfigure: left, $10\times 20$p$4$;
right $20\times 80$p$4$. The interface is highlighted in black.

Clearly, the flow field arising from the Rayleigh-Taylor instability is much
richer than the previous test case: viscous and gravitational forces lead to
complex vortical structures causing non-trivial dynamics of the interface,
including primary breakup. These structures are increasingly better resolved
under mesh refinement as it can be observed from figure 11.

From a more quantitive point of view, the predicted top and bottom locations
of the interface versus the non-dimensional time
($t^{*}=t/\sqrt{d/(||\textbf{g}||\mathrm{At})}$) are shown in figure 12.
Excellent agreement with previous studies can be observed.

Figure 12: The evolution of the top and bottom of the interface for the
Rayleigh-Taylor instability versus non-dimensional time for $10\times 40$p$4$
and $20\times 80$p$4$. Symbols indicated the reference solution by Chiu & Lin
[33].

Similarly to the Rider-Kothe test case, the influence of the order of
approximation was investigated. In particular, a $1^{\mathrm{st}}$ polynomial
order simulation has been performed and compared with the $4^{\mathrm{th}}$
polynomial order computation. A series of snapshots of the phase field are
shown in figure 13 for the two different orders of approximation.

It can be noticed that the p$1$ simulation provides less accurate results,
leading to an asymmetric solution.

Figure 13: Snapshots of the phase field isocontours for different orders and
times. Top row: isocontour of the phase field at $t^{*}=0.5$ (left) and
$t^{*}=1.5$ (right). Bottom row: isocontour of the phase field at $t^{*}=2.5$
(left) and $t^{*}=3.0$ (right). In each subfigure: left, $25\times 100$p$1$;
right, $10\times 40$p$4$. The interface is highlighted in black.

From a more quantitative point of view, the location of the upper and lower
plumes were compared for the two different orders of approximations and are
shown in figure 14. It can be noticed that the p$4$ approximation is always
closer to the reference solution throughout the whole simulation, in
particular for the prediction of the top location of the interface.

Figure 14: The evolution of the top and bottom of the interface for the
Rayleigh-Taylor instability versus non-dimensional time using a p$1$
($50\times 200$p$1$) and p$4$ ($20\times 80$p$4$) approximation. Symbols
indicate the reference solution by Chiu & Lin [33].

To better appreciate the differences between high and low orders of
approximation, the locations of the interface at $t^{*}=1.5$ and $t^{*}=3.0$
are shown in figures 15 and 16 at different resolutions and different
polynomial orders. It can be noticed that the p$1$ approximation is
characterised by smoother, over-dissipated profiles along the tangential
direction of the interface. Remember that the smoothness along the normal
direction is governed by the interface capturing technique and it should not
be particularly influenced by the spatial order of approximation, but the
accuracy along the tangential direction of the interface, instead, should be
more sensitive to the approximation order. At later times, when primary beak-
up occurs, more complex small scale structures are captured by the high-order
simulation, in particular for the simulations on the finest grid. It can also
be noticed that the p$4$ simulation is characterised by a more stretched
interface, (_i.e._ bottom and top locations of the interface are respectively
lower and higher with respect to the corresponding p$1$ simulations). This
behaviour is in agreement with what is observed in figure 14.

Figure 15: Iso-contour $\phi=0.5$ with $1^{\mathrm{st}}$ and $4^{\mathrm{th}}$
polynomial order approximations at $t^{*}=1.5$. Left, $50\times 200$ DoF;
right $100\times 400$ DoF. Figure 16: Iso-contour $\phi=0.5$ with
$1^{\mathrm{st}}$ and $4^{\mathrm{th}}$ polynomial order approximations at
$t^{*}=3.0$. Left, $50\times 200$ DoF; right $100\times 400$ DoF.

####  4.2.3 Rising bubble

In this test case the rise of a $2$D bubble of light fluid within a heavier
fluid due to buoyancy is simulated. Initially, a circular bubble of radius
$R=0.25$ is placed at $(0.0,0.5)$ in a $[-0.5,0.5]\times[0,2]$ domain. The
density and viscosity of the fluids are chosen such that
$\rho_{1}/\rho_{2}=\mu_{1}/\mu_{2}=10$. The Reynolds number is set to
$\mathrm{Re}=(\rho_{1}d^{3/2}||\textbf{g}||^{1/2})/\mu_{1}=35$ whereas the
Eötvös number is $E_{0}=4\rho_{1}||\textbf{g}||r^{2}/\sigma=10$. Similarly to
the previous test case, at the upper boundary a Riemann-invariant boundary
condition with zero velocity and constant pressure is imposed and the bottom
boundary is a no-slip wall boundary. Left and right boundaries are treated as
slip walls.

Notice that the reference solution coincides with the results presented in
[48].

Figure 17 shows the evolution of the bubble in time as it rises due to
buoyancy.

Figure 17: Location of the interface at different times for the most refined
simulation ($20\times 40$p$4$). In order, $t=0.0,1.0,2.0,3.0$. At final time,
symbols indicate the reference solution by Hysing et al. [48].

Finally, the location of the center of gravity of the bubble and its vertical
velocity are shown in figure 18. In the work by Manzanero et al. [49], a
similar numerical framework based on the DG method was presented. In their
case, a weakly compressible approach was employed and coupled with the Cahn-
Hilliard equations. It can be noticed that the pressure reflections coming
from the wall boundaries cause visible oscillations in the bubble’s rising
velocity. Our simulation produces a similar dynamics of the bubble without
spurious oscillations. In our case, a Riemann invariant boundary condition was
prescribed on the upper boundary, avoiding undesirable pressure reflections.

Figure 18: Location of the bubble’s center of gravity (left) and bubble’s
vertical velocity (right) over time for the present simulation $10\times
20$p$4$, $20\times 40$p$4$ and results by Manzanero et al. [49].

A more challenging test case with larger density ratio and smaller surface
tension coefficient was also considered. In particular, the density and
viscosity of the two phases were chosen such that $\rho_{1}/\rho_{2}=1000$ and
$\mu_{1}/\mu_{2}=100$. Furthermore, the Eötvös number was increased to
$E_{0}=125$ whereas the Reynolds number remained the same. Under these
conditions, the surface tension effects are not strong enough to prevent
strong deformations of the bubble.

In figure 19, the evolution of the bubble is shown. It can be seen that at
late times, the edges of the bubble tend to elongate in considerably thin
ligaments. This behaviour is caused by the insufficiently small surface
tension forces with respect to gravitational forces.

Figure 19: Location of the interface at different times for the most refined
simulation ($20\times 40$p$4$). In order, $t=0.0,1.0,2.0,3.0$. At final time,
symbols indicate the reference solution from [49].

Similarly with respect to the previous case, the final location of the
interface was compared with a reference solution. It is well known that for
this specific test case, different codes produce slightly different results in
terms of the elongated ligaments. In figure 19, at the final time, the present
simulation is consequently compared with the results by Manzanero et al. [49],
who proposed a similar diffused-interface approach based on the Cahn-Hilliard
equations. From this comparison, it can be seen that the predicted location of
the interface agrees well with the reference solution.

Finally, the location of the bubble’s center of gravity and the mean rising
velocity of the bubble are evaluated over time and are shown in figure 20. The
results are compared with classical reference solution from sharp interface
solvers [48] and the diffused-interface method from Manzanero et al. [49]. It
can be seen that a better agreement is obtained with this latter reference,
due to the similarities between the two approaches.

Figure 20: Location of the bubble’s center of gravity (left) and bubble’s
vertical velocity (right) over time for the present simulation $10\times
20$p$4$, $20\times 40$p$4$ and results by Manzanero et al. [49].

##  5 Conclusions

A high-order numerical approach based on the Discontinuous Galerkin method was
proposed for the simulation of two-phase flows. The interface between the two
phases was modelled using the conservative Allen-Cahn equation, which was
subsequently implemented within a five equation model to describe the motion
of two immiscible compressible fluids.

A series of benchmark cases were considered, including both kinematic tests
involving the numerical resolution of the Allen-Cahn equation only and of
fully coupled five equation model for the simulation of two-phase flows.

In kinematic tests, studying the same approach for different orders of
approximation, it was found that higher polynomial orders provided significant
advantages in terms of overall accuracy, including smaller $L_{1}$ errors and
mass conservation errors. The same simulations were also performed with both
quadrilateral and triangular elements without any ad-hoc modifications of the
model’s parameters, further underlining the robustness and flexibility of the
proposed approach.

Similar observations were made in the simulation of two-phase flows using the
five equation model. Lower orders of approximation were characterised by more
pronounced spurious features such as deformations of the interface and
artificial break-up. Finally, more complex problems involving viscous effects,
gravitational forces, and surface tension were considered. The proposed
methodology successfully recovered accurate solutions in agreement with
existing simulations in the literature.

Overall, the proposed methodology showed the capability of preserving
fundamental features of two-phase flows such as the boundedness of the phase
field, accurate computation of the interface normal vectors, small mass
conservation errors and exact resolution of contact discontinuities, without,
at the same time, giving up the significant benefits coming from high-order
spatial discretisations. All of this is achieved with only minor modifications
of the central core of the DG method, by choosing appropriate numerical fluxes
and parameters of the model. At the same time, because of this, the
methodology herein presented did not show the need of ad-hoc modifications for
different type of elements or orders of approximation further emphasising its
robustness, flexibility and generality.

## Acknowledgments

Financial support from NSF (Award 1909379) and Daikin is greatly appreciated.

## References

  * [1] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, D. Mavriplis, CFD vision 2030 study: A path to revolutionary computational aerosciences, Tech. rep., NASA Langley Research Center, Hampton, Virginia 23681-2199, NASA Contractor Report 218178 (2014).
  * [2] J. S. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer, 2008.
  * [3] B. Cockburn, C. Shu, The local Discontinuous Galerkin finite element method for convection-diffusion systems, SIAM Journal of Numerical Analysis 35 (1998) 2440–2463.
  * [4] B. Cockburn, C. Shu, The Runge-Kutta Discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems, Journal of Computational Physics 141 (1998) 199–224.
  * [5] C. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. D. Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, B. Jordi, H. Xu, Y. Mohamied, C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R. Kirby, S. Sherwin, Nektar++: An open-source spectral/hp element framework, Computer Physics Communications 192 (2015) 205 – 219.
  * [6] H. T. Huynh, A flux reconstruction approach to high-order schemes including Discontinuous Galerkin methods, in: AIAA Paper, 2007-4079, 2007.
  * [7] P. E. Vincent, P. Castonguay, A. Jameson, A new class of high-order energy stable flux reconstruction schemes, Journal of Scientific Computing 47 (2011) 50–72.
  * [8] F. D. Witherden, A. M. Farrington, P. E. Vincent, PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach, Computer Physics Communications 185 (2014) 3028–3040.
  * [9] D. A. Kopriva, J. H. Kolias, A conservative staggered-grid Chebyshev multidomain method for compressible flows, Journal of Computational Physics 125 (1996) 244–261.
  * [10] G. Karniadakis, S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford University Press, 2013.
  * [11] I. Bosnyakov, S. V. Lyapunov, A. Troshin, V. V. Vlasenko, A. V. Wolkov, C. Hirsch, A high-order Discontinuous Galerkin method for external aerodynamics, in: AIAA Paper, 2014.
  * [12] Y. Lv, M. Ihme, Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion, Journal of Computational Physics 270 (2014) 105–137.
  * [13] G. Mengaldo, D. Moxey, M. Turner, R. C. Moura, A. Jassim, M. Taylor, J. Peiró, S. Sherwin, Industry-relevant implicit large-eddy simulation of a high-performance road car via spectral/hp element methods, SIAM Review 63 (2021) 723–755.
  * [14] N. Tonicello, G. Lodato, L. Vervisch, Analysis of high-order explicit LES dynamic modeling applied to airfoil flows, Flow, Turbulence and Combustion 108 (2022) 77–104.
  * [15] N. Tonicello, G. Lodato, L. Vervisch, Turbulence kinetic energy transfers in direct numerical simulation of shock-wave–turbulence interaction in a compression/expansion ramp, Journal of Fluid Mechanics 935 (2022) A31.
  * [16] E. J. Ching, Y. Lv, P. Gnoffo, M. Barnhardt, M. Ihme, Shock capturing for Discontinuous Galerkin methods with application to predicting heat transfer in hypersonic flows, Journal of Computational Physics 376 (2019) 54–75.
  * [17] D. Moxey, C. D. Cantwell, Y. Bao, A. Cassinelli, G. Castiglioni, S. Chun, E. Juda, E. Kazemi, K. Lackhove, J. Marcon, et al., Nektar++: Enhancing the capability and application of high-fidelity spectral/hp element methods, Computer Physics Communications 249 (2020) 107110.
  * [18] Y. Lv, X. I. Yang, G. I. Park, M. Ihme, A Discontinuous Galerkin method for wall-modeled large-eddy simulations, Computers & Fluids 222 (2021) 104933.
  * [19] T. Dzanic, F. D. Witherden, Positivity-preserving entropy-based adaptive filtering for Discontinuous spectral element methods, Journal of Computational Physics 468 (2022) 111501.
  * [20] E. Ferrer, G. Rubio, G. Ntoukas, W. Laskowski, O. Mariño, S. Colombo, A. Mateo-Gabín, H. Marbona, F. M. de Lara, D. Huergo, et al., HORSE3D: A high-order Discontinuous Galerkin solver for flow simulations and multi-physics applications, Computer Physics Communications 287 (2023) 108700.
  * [21] S. K. Lele, Compact finite difference schemes with spectral-like resolution, Journal of Computational Physics 103 (1992) 16–42.
  * [22] C. Bogey, C. Bailly, A family of low dispersive and low dissipative explicit schemes for flow and noise computations, Journal of Computational physics 194 (2004) 194–214.
  * [23] K. Van den Abeele, C. Lacor, Z. J. Wang, On the stability and accuracy of the spectral difference method, Journal of Scientific Computing 37 (2008) 162–188.
  * [24] P. E. Vincent, P. Castonguay, A. Jameson, Insights from von Neumann analysis of high-order flux reconstruction schemes, Journal of Computational Physics 230 (2011) 8134–8154.
  * [25] R. C. Moura, S. J. Sherwin, J. Peiró, Linear dispersion–diffusion analysis and its application to under-resolved turbulence simulations using Discontinuous Galerkin spectral/hp methods, Journal of Computational Physics 298 (2015) 695–710.
  * [26] J. Vanharen, G. Puigt, X. Vasseur, J.-F. Boussuge, P. Sagaut, Revisiting the spectral analysis for high-order spectral Discontinuous methods, Journal of Computational Physics 337 (2017) 379–402.
  * [27] G. Mengaldo, R. Moura, B. Giralda, J. Peiró, S. Sherwin, Spatial eigensolution analysis of Discontinuous Galerkin schemes with practical insights for under-resolved computations and implicit LES, Computers & Fluids 169 (2018) 349–364.
  * [28] N. Tonicello, G. Lodato, L. Vervisch, A comparative study from spectral analyses of high-order methods with non-constant advection velocities, Journal of Scientific Computing 87 (2021) 1–38.
  * [29] R. DeBar, Fundamentals of the KRAKEN code, Technical Report UCIR-760.
  * [30] B. Nichols, C. Hirt, Methods for calculating multidimensional, transient free surface flows past bodies, in: Proceedings of the First International Conference on Numerical Ship Hydrodynamics, Vol. 20, Naval Ship Research and Development Center Bethesda, MD, USA, 1975.
  * [31] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, Journal of Computational Physics 210 (2005) 225–246.
  * [32] Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, Journal of Computational Physics 220 (2007) 626–653.
  * [33] P.-H. Chiu, Y.-T. Lin, A conservative phase field method for solving incompressible two-phase flows, Journal of Computational Physics 230 (2011) 185–204.
  * [34] J. U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational Physics 100 (1992) 335–354.
  * [35] G. Allaire, S. Clerc, S. Kokh, A five-equation model for the simulation of interfaces between compressible fluids, Journal of Computational Physics 181 (2002) 577–616.
  * [36] J. Al-Salami, M. M. Kamra, C. Hu, A high order flux reconstruction interface capturing method with a phase field preconditioning procedure, Journal of Computational Physics 438 (2021) 110376.
  * [37] F. H. Harlow, A. A. Amsden, Fluid dynamics. A LASL monograph, Tech. rep., Los Alamos National Lab (LANL), Los Alamos, NM (United States) (1971).
  * [38] J. Cheng, F. Zhang, T. Liu, A quasi-conservative Discontinuous Galerkin method for solving five equation model of compressible two-medium flows, Journal of Scientific Computing 85 (2020) 1–35.
  * [39] E. J. Ching, B. Bornhoft, A. Lasemi, M. Ihme, Quail: A lightweight open-source Discontinuous Galerkin code in python for teaching and prototyping, SoftwareX 17 (2022) 100982.
  * [40] R. Hartmann, P. Houston, An optimal order interior penalty Discontinuous Galerkin discretization of the compressible Navier–Stokes equations, Journal of Computational Physics 227 (2008) 9670–9685.
  * [41] W. J. Rider, D. B. Kothe, Reconstructing volume tracking, Journal of Computational Physics 141 (1998) 112–152.
  * [42] S. Mirjalili, C. B. Ivey, A. Mani, A conservative diffuse interface method for two-phase flows with provable boundedness properties, Journal of Computational Physics 401 (2020) 109006.
  * [43] R. Saurel, R. Abgrall, A simple method for compressible multifluid flows, SIAM Journal on Scientific Computing 21 (1999) 1115–1145.
  * [44] E. Johnsen, T. Colonius, Implementation of WENO schemes in compressible multicomponent flow problems, Journal of Computational Physics 219 (2006) 715–732.
  * [45] R. Saurel, F. Petitpas, R. A. Berry, Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures, Journal of Computational Physics 228 (2009) 1678–1712.
  * [46] A. Tiwari, J. B. Freund, C. Pantano, A diffuse interface model with immiscibility preservation, Journal of Computational Physics 252 (2013) 290–309.
  * [47] M. Aslani, J. D. Regele, A localized artificial diffusivity method to simulate compressible multiphase flows using the stiffened gas equation of state, International Journal for Numerical Methods in Fluids 88 (2018) 413–433.
  * [48] S.-R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska, Quantitative benchmark computations of two-dimensional bubble dynamics, International Journal for Numerical Methods in Fluids 60 (2009) 1259–1288.
  * [49] J. Manzanero, G. Rubio, D. A. Kopriva, E. Ferrer, E. Valero, Entropy–stable Discontinuous Galerkin approximation with summation–by–parts property for the incompressible Navier–Stokes/Cahn–Hilliard system, Journal of Computational Physics 408 (2020) 109363.