High-order ADER Discontinuous Galerkin schemes for a symmetric hyperbolic model of compressible barotropic two-fluid flows

This paper presents a high-order discontinuous Galerkin finite element method to solve the barotropic version of the conservative symmetric hyperbolic and thermodynamically compatible (SHTC) model of compressible two-phase flow, introduced by Romenski et al., in multiple space dimensions. In the absence of algebraic source terms, the model is endowed with a curl constraint on the relative velocity field. In this paper, the hyperbolicity of the system is studied for the first time in the multidimensional case, showing that the original model is only weakly hyperbolic in multiple space dimensions. To restore strong hyperbolicity, two different methodologies are used: i) the explicit symmetrization of the system, which can be achieved by adding terms that contain linear combinations of the curl involution, similar to the Godunov-Powell terms in the MHD equations; ii) the use of the hyperbolic generalized Lagrangian multiplier (GLM) curl-cleaning approach forwarded. The PDE system is solved using a high-order ADER discontinuous Galerkin method with a posteriori sub-cell finite volume limiter to deal with shock waves and the steep gradients in the volume fraction commonly appearing in the solutions of this type of model. To illustrate the performance of the method, several different test cases and benchmark problems have been run, showing the high-order of the scheme and the good agreement when compared to reference solutions computed with other well-known methods.

1 Introduction original BN model and which describe the so-called lift forces.Recently, an exact solution for the corresponding Riemann problem in the barotropic case was found in [31], and an all-Mach number flow solver was developed in [32].
The model proposed by Romenski et al. is strongly hyperbolic in one space dimension, and its homogeneous part without algebraic source terms is endowed with a curl involution on the relative velocity field.However, as we will show in this paper, the original model is only weakly hyperbolic in multiple space dimensions.In order to restore strong hyperbolicity, two different strategies can be followed: the first one consists in using the hyperbolic generalized Lagrangian multiplier (GLM) curl-cleaning methodology introduced in [33,34,35,36], which is a natural extension of the original ideas presented by Munz et al. in [37,38] on hyperbolic GLM divergence cleaning for the Maxwell and MHD equations, which contain the well-known divergence-free condition of the magnetic field.Curl constraints can also be found in many other first-order hyperbolic models, such as hyperbolic models for surface tension [39,34,40], first-order hyperbolic reformulations of the Navier-Stokes-Korteweg equations based on an augmented Lagrangian approach [36,41], or first-order reductions of the Einstein field equations of general relativity [33].The second methodology to recover the hyperbolicity involves the use of some extra terms in the momentum equation that symmetrize the system and is therefore directly based on the theory of Symmetric Hyperbolic and Thermodynamically Compatible (SHTC) systems, following [42,43,44,45].In both cases, we will show that the system in the multidimensional case becomes again strongly hyperbolic.
In the setting of this model, the solutions are often discontinuous in space, that is, solutions to Riemann problems.Therefore, it is necessary to consider an approach that is robust and accurate even in the presence of shock waves or discontinuities.To address the sharp gradients of the numerical solutions, a high-order ADER Discontinuous Galerkin (DG) finite element framework with a posteriori sub-cell finite volume (FV) limiter is considered, see [46,47,48] for further details.The proposed method is high-order in space and time thanks to the use of the ADER approach of Toro and Titarev [49,50,51].To deal with spurious oscillations that may appear in the presence of discontinuities or shock waves, it makes use of an a posteriori subcell finite volume limiter for high-order fully discrete one-step ADER-DG schemes presented in [47,48], which follows the MOOD approach of Clain and Loubère [52,53,54].
This paper is organized as follows.Section 2 recalls the set of governing partial differential equations and shows the EOS that we will use in this work.In Section 2.2, the hyperbolicity of the system is studied in the multidimensional case, showing that the original model is only weakly hyperbolic, and two different strategies to recover strong hyperbolicity are presented.Section 3 introduces the high-order ADER discontinuous Galerkin scheme used in this paper to solve the model numerically.Section 4 shows the results of several test cases and benchmark problems computed using the proposed method.Finally, Section 5 concludes the paper and summarizes the contributions of this work as well as its possible future extensions.

Governing equations
The governing equations of the barotropic compressible two-velocity, two-pressure twofluid model of Romenski, written in terms of the specific total energy potential E = E(α1, c1, ρ, w k ) are given by where α1 and α2 are the volume fractions of the first and second phases, verifying the saturation constraint α1 + α2 = 1, ρ1 and ρ2 are the mass densities of the first and second phase, ρ = α1ρ1 + α2ρ2 is the mixture mass density, cj = αjρj/ρ with j = 1, 2 are the mass fractions of phase j which satisfy c1 + c2 = 1, and g = (g 1 , g 2 , g 3 ) T is the gravity acceleration.If u1 = (u 1 1 , u 2 1 , u 3 1 ) T and u2 = (u 1 2 , u 2 2 , u 3 2 ) T are the velocity vectors of each phase, then the mixture velocity u = (u 1 , u 2 , u 3 ) T is computed as u k = c1u k 1 +c2u k 2 and the relative phase velocity, which is a primary evolution quantity in this model, is given by w The source terms of equations (1a) and (1e) are proportional to thermodynamic forces, with τ the rate of pressure relaxation and ζ the inter-phase friction coefficient.The PDE system (1) is formed by two conservation laws for the volume and mass fractions, Eq. ( 1a) and (1b), respectively, conservation of total mass, Eq. (1c), conservation of the mixture momentum, Eq. (1d), and a balance law for the relative velocity, Eq. (1e).The algebraic source terms appearing in Eq. (1a) and (1e), describe the interaction between the phases and are pressure relaxation and interfacial friction.The mixture equation of state is given by, see [31], where the specific internal energy of the mixture reads as where ej(ρj) is the specific internal energy of the phase j.Since in this work an isentropic process is considered, the derivatives of e can be computed as where hj(ρj) = ej(ρj) + p j (ρ j ) ρ j , j = 1, 2 is the specific enthalpy of the phase j.Then, the derivatives of the specific total energy E are given by Considering the computations in [31], and taking into account ( 2) and (3), the following identities are obtained The PDE system (1) can now be rewritten in the following form, which is more convenient for numerical discretization: The derivation of the PDE system ( 4) is based on the principles of thermodynamically compatible systems [26], and it consists of nine equations: the balance law for the volume fraction, the conservation laws of the two-phase masses, conservation of mixture momentum and the balance law for the relative velocity field.Note that in the absence of source terms in (4e), i.e., for ζ = 0, the relative velocity is curl-free in the sense

Equation of state (EOS)
In order to close the two-phase model (4a)-(4e), it is necessary to define an equation of state for each phase: throughout this paper, we will make use either of the ideal gas law or of the stiffened gas equation of state, which will be used to model a liquid phase.For ideal gases, the EOS is defined as where γ is the adiabatic index or the ratio of specific heats, c0 is the adiabatic sound speed, s is the specific entropy (which in our case will be constant), cv is the specific heat capacity at constant volume, and the pressure is given by For stiffened gases, the EOS reads as where ρ0 and p0 are the reference density and pressure, respectively, and c0 is a constant reference sound speed.In this case, the pressure is computed as

Hyperbolicity analysis
In this section, we will study the hyperbolicity of the model (4a)-(4e).The onedimensional case was already addressed in [31], and here we briefly recall some results for the one-dimensional case before moving to the more general multidimensional case.In particular, we will prove that the original system (4a)-(4e) is only weakly hyperbolic in the multidimensional case and show how the strong hyperbolicity can be restored considering two different strategies: on the one hand, using an extension of the hyperbolic Generalized Lagrangian Multiplier (GLM) curl-cleaning approach and on the other hand, modifying the original system of governing equations by adding the symmetrizing terms that allow rewriting the model in symmetric hyperbolic form, which is natural within the framework of SHTC equations.Defining the vectors of conserved and primitive variables as Q = (α1, α1ρ1, α2ρ2, ρu T , w T ) T and V = (α1, ρ1, ρ2, u T 1 , u T 2 ) T , respectively, the PDE system (4) can be written as where S(Q) contains the source terms, F(Q) is the nonlinear flux tensor and B(Q)•∇Q contains the non-conservative terms.Then, the quasi-linear form of the PDE in terms of the conserved variables Q is given by where If the vector of primitive variables V is considered, the system can be written as where for the sake of readability and since it does not contribute anything to the study of the hyperbolicity of the system, the source terms are set to 0.

One-dimensional case
The hyperbolicity analysis of the system (4) in 1D has been done in detail in [31], so here, only a summary of the main points will be made that will be useful for a better understanding of the multidimensional case.If u = c1u1 + c2u2 is the mixture velocity, with uj, j = 1, 2 the velocity of the phase j and w = u1 − u2 is the relative velocity, the one-dimensional system results Then, the matrix C in ( 9) is given by , with aj the sound speed of phase j, that is defined as a 2 j = ρj It is easy to show that C(V) admits five eigenvalues, whose expressions are given hereafter In this case, all five eigenvalues are real, and a complete set of five linearly independent eigenvectors exists, which means that the system in one space dimension is strongly hyperbolic, see [31] for further details about the eigenvalues and eigenvectors in the one-dimensional case.

Multidimensional case
Now we are in the position to study the multidimensional case.In the following, we use of the property of rotational invariance of Newtonian mechanics, hence it is enough to consider the matrix C(V) only in the x−direction and not all possible space directions.The system under consideration is (4) and, in this case, the matrix C(V) in the x−direction is given by As in the previous section, we make use of the following auxiliary variables to ease the expressions in the matrix The matrix C admits 9 eigenvalues λ1−9 that are given by where the sound speed of each phase aj again is defined as a 2 j = ρi ∂h j ∂ρ j , j = 1, 2. The eigenvalues are all real.To prove whether the system is weakly or strongly hyperbolic, it is necessary to compute the associated eigenvectors.The right eigenvectors are the columns in the matrix below and are given in the same order as the eigenvalues: , where the auxiliary variables that are used to ease the notation are defined as All eigenvalues are real, but two eigenvectors are missing (in fact, if we are working in dimension d there are d − 1 missing eigenvectors), so the system is only weakly hyperbolic.Hereafter we will show two different methodologies that can be used to restore the strong hyperbolicity of the model.

Generalized Lagrangian multiplier (GLM) curl-cleaning approach
To restore the strong hyperbolicity and following the same strategy that can be found in [33,34,35,36], we make use of the GLM curl-cleaning technique, where an evolution equation for a curl-cleaning field ψ k is added to the system (4).Using the abbreviation δ12 = 1 2 (u l 1 ) 2 − 1 2 (u l 2 ) 2 +h1−h2 and denoting the curl-cleaning speed by a ψ this equation is coupled with (4e) via a Maxwell-type sub-system as follows: hence the augmented system with GLM curl-cleaning reads where ψ = (ψ 1 , ψ 2 , ψ 3 ) is the cleaning field, ε = ε klm is the Levi-Civita tensor and a ψ is the curl-cleaning speed.Once the curl-cleaning field has been added, we can compute the eigenvalues and eigenvectors to check the hyperbolicity of the augmented GLM system.The matrix C in x−direction can be written as , where, as in the previous case, we make use of the auxiliary variables (11) to ease the notation of the matrix.In this case, it is easy to compute the twelve eigenvalues of the matrix C that are given by with aj the sound speed of phase j, that is defined as a 2 j = ρj ∂h j ∂ρ j , j = 1, 2. Since the eigenvalues are real, to check if the system is weakly or strongly hyperbolic, it is necessary to compute the associated eigenvectors.Below, we will write the matrix that contains the right eigenvectors in columns.They are listed in the same order as the eigenvalues.
, where the auxiliary variables that have been used to write the twelve right eigenvectors more compactly are defined as Using the GLM curl-cleaning technique, we obtain twelve real eigenvalues and the corresponding twelve linearly independent right eigenvectors, hence the augmented system with GLM curl-cleaning ( 12) is strongly hyperbolic.

Symmetrizing Godunov-Powell terms
The second strategy used to recover the strong hyperbolicity of system ( 4) is based on the theory of SHTC systems and consists in adding terms that are proportional to the curl involution to the momentum equation (4d) so that the system has real eigenvalues and a complete set of linearly independent eigenvectors, see [34].The modified system reads Using (2), system (13) results in As in the previous case, the eigenvalues and eigenvectors are computed to analyze the hyperbolicity of the system (14).The matrix C in the x−direction reads As previously, we have used the auxiliary variables (11) to lighten the matrix.In this case, the nine eigenvalues of the matrix C are given by with aj the sound speed of phase j, which is defined as a 2 j = ρj ∂h j ∂ρ j , j = 1, 2. Since the eigenvalues are all real, to check if the system is weakly or strongly hyperbolic, it is necessary to compute the associated eigenvectors.Below, we will write the matrix that contains the right eigenvectors in columns.They are listed in the same order as the eigenvalues.
, where the auxiliary variables used to write the nine right eigenvectors more compactly are defined as Since we have obtained nine real eigenvalues with a full set of linearly independent eigenvectors, the system (14) with the additional symmetrizing Godunov-Powell terms is strongly hyperbolic.
3 High-order ADER discontinuous Galerkin finite element scheme with a posteriori subcell finite volume limiter As described in Section 2 and references [28,29,31], the system (4) is a hyperbolic system that can be compactly written as where Q is the vector of conserved variables, are the non-conservative terms and S is the vector of algebraic source terms.To solve the system (15), we will use high-order ADER discontinuous Galerkin schemes with a posteriori sub-cell finite volume limiter on uniform Cartesian meshes, see [46,55,56,57,47,48] for further details.

One-Step ADER-DG Schemes
In the following, a description of the method is given for the two-dimensional case (d = 2).The computational domain Ω = [−Lx/2, Lx/2] × [−Ly/2, Ly/2] is discretized with a Cartesian grid composed of Nx × Ny cells.These cells are given by Ωi , with (xi, yi) the barycenter of the cell Ωi, and ∆x = Lx Nx , ∆y =

Ly
Ny the mesh spacing in the x− and y−directions.Let u h (x, t n ) be the discrete solution of (15) in each spatial control volume Ωi at time t n , written in terms of tensor products of piecewise polynomials of degree N , and let V h be the space of tensor products of piecewise polynomials of the degree up to N .Then, the discrete solution u h can be written in terms of the basis functions, ϕ l (x, y), l ∈ [1, (N + 1) d ], in every cell Ωi as where ϕ l = ϕ l (x) are the basis functions associated with V h .We take an orthogonal nodal basis {ϕ l } l∈{0,...,(N +1) d } , generated by the tensor product {ϕ l 1 ϕ l 2 ϕ l 3 } l 1 ,l 2 ,l 3 ∈{0,...,N } where {ϕ l l } l l ,∈{0,...,N } are the Lagrange interpolation polynomials going through the N + 1 Gauss-Legendre quadrature nodes.Multiplying (15) by a test function ϕ l ∈ V h , and integrating the equation over the space-time control volume Ωi × [t n , t n+1 ], the weak formulation can be written as To achieve high-order in space and time, an ADER approach can be used.This methodology was put forward by Toro et al. for the first time in [58] for linear problems on Cartesian meshes, it can be implemented in both the finite volume and the discontinuous Galerkin finite element framework and is uniformly and arbitrarily high-order accurate in both space and time.In this work, an alternative version of the ADER approach is considered, avoiding the use of the Cauchy-Kovalevskaya procedure by using a local space-time discontinuous Galerkin predictor, which is based on a weak space-time formulation of the governing PDE (17).
Using (16) and integrating the term with the time derivative by parts in time and the divergence term by parts in space, then (17) where n is the outward unit normal vector at the cell boundary ∂Ωi, q h is a local space-time predictor, which will be explained below, q + h and q − h are the boundaryextrapolated values of the space-time predictor from within Ωi and its neighbor Ωj.Usually, q h presents jumps across the cell boundaries which can be resolved by the solution of a generalized Riemann problem (see [59,60] for more details).In (18), G denotes the Riemann solver (numerical flux function), which depends on the left state q − h and the right state q + h .In this case, this integral has been approximated by using the Rusanov flux, see [61] where smax = max λ k (q + h ) , λ k (q − h ) is the maximum wave speed at the interface.To deal with the jump terms in the non-conservative product, a path-conservative method is employed, following [62,63].In this setting, a straight-line segment path is chosen Ψ(s, q + h , q − h ) = q − h + s(q + h − q + h ), s ∈ [0, 1], so that the non-conservative terms reduce to the following expression

Local space-time predictor
In the following, we describe the local space-time predictor used to compute the coefficients ûn l,i in Equation ( 18).We consider space-time basis functions θ l , that are obtained as the tensor product θ l (x, t) = ϕ k 0 (t)ϕ k 1 (x), of the same previously introduced Lagrange interpolation polynomials, just that now the basis functions also depend on time.The predictor q h is written in the form as a weak solution to (15).Then, using (20) in (15), multiplying by a space-time basis function θ l and integrating over Ωi × [t n , t n+1 ] yields Integrating by parts only the first term on the left-hand side and taking into account that at time t n we start from the known state u n h allows us to write Equation ( 21) is a system for the unknowns q of the space-time predictor q h (x, t) and can be computed in terms of the spatial degrees of freedom ûn l .It is solved by a fixed point iteration for which convergence was proven in [11].Once the predictor is known, Equation ( 18) allows to compute the polynomial coefficients ûn+1 in each cell by using Gaussian quadrature for the remaining integrals.

A posteriori subcell finite volume limiter
Although the numerical method presented in the previous section is a high-order method, it is linear in the sense of Godunov, which means that spurious oscillations will appear in the presence of discontinuities or shock waves.To overcome this problem, we use the a posteriori subcell limiter for high-order fully discrete one-step ADER-DG schemes presented in [47,48].This subcell FV limiter is based on the MOOD paradigm introduced in [52,53,54] for finite volume schemes.
The scheme described in the previous section is run over the entire domain at each time step, and a so-called candidate solution u * h (x, t n+1 ) is obtained.Then, one checks whether the candidate solution verifies some numerical and physical detection criteria (positivity of the densities ρ1 and ρ2, α1 with values between 0 and 1) and whether the discrete maximum principle (DMP), [47], is verified.If a cell Ωi violates any of the above criteria, that cell is flagged as a troubled cell and for the application of the subcell finite volume limiter.The limiter is denoted as a posteriori because it is applied after the candidate solution has been computed.
The limiter is applied in the following way: all cells Ωi marked as troubled are subdivided into (2N +1) d subcells, which are denoted by Ωi,j where Ωi = j Ωi,j.The discrete solution at time t n is given by the piecewise constant cell averages, denoted by ūn i,j .They are obtained from the high-order DG polynomials u h (x, t n ) by averaging using the definition of the cell average It is worth noting that subdividing a high-order DG element into 2N + 1 FV subcells per space dimension does not reduce the time step size of the overall scheme since the CFL stability condition of explicit DG schemes scales with 1/(2N + 1) in 1D, while the maximum Courant number of finite volume methods is unity in one space dimension.The cell averages (22) are evolved in time using either a second-order MUSCL-Hancock-type TVD finite-volume scheme with minmod limiter, or by making use of a third-order ADER-WENO finite-volume scheme, see [47], which are also both predictor-corrector methods, and thus look almost identical to the ADER-DG scheme, except for the necessary nonlinear reconstruction step.Moreover, in this case, the test function is unity, which implies that the volume integral over the flux term disappears and the volumes computed over Ωi are replaced by the volumes over the subcells Ωi,j, hence |Ωi,j| ūn+1 i,j − ūn i,j + The limited DG polynomial u h at time t n+1 is then obtained by performing a constrained least squares reconstruction and using the averages of all the subcells of Ωi computed using (23).The reconstruction reads ∀Ωi,j ∈ Ωi, with the linear constraint The constraint (24) means conservation of the solution within the element Ωi.In addition to the expansion coefficients ûn+1 i,l of the limited DG polynomial, in all limited DG elements we also keep in memory the averages of the finite volume subcells ūn+1 i,j , as they serve as initial condition for the finite volume limiter of the subcell in case a cell is problematic also in the next time step, see [47].More details about the a posteriori subcell finite volume limiter can be found in [47,48,64].

Numerical results
This section is devoted to showing some test cases to illustrate the high order of accuracy of the proposed method, especially in the presence of steep gradients in the solution.First, some simulations are performed to show the experimental order of convergence (EOC) of the proposed ADER-DG method.Then, our scheme is used to solve some Riemann problems in 1D and 2D.Finally, a dambreak problem is simulated, and the results are compared with those obtained with a reduced Baer-Nunziato model.All test cases have been performed with the system (4) without using the curl-cleaning technique or the symmetrizing Godunov-Powell terms.Although the original system is only weakly hyperbolic, no stability problems have been found in the numerical simulations, contrary to what was reported in [34,40,41] for other weakly hyperbolic systems with curl involutions.Moreover, in all tests, the algebraic relaxation source terms have been neglected, and the gravity g is set to 0, except in the dambreak test case, where it is necessary to consider the gravity.

Accuracy analysis
This section performs a numerical convergence analysis to show the experimental order of convergence of the proposed ADER-DG method.To construct an exact solution, following [65,66], an analytical, stationary, and rotationally symmetric solution of the system (4) is computed considering cylindrical coordinates (r, θ and z), with (u r , u θ , u z ) the velocity vector and (w r , w θ , w z ) the relative velocity vector.The analytical solution is assumed to approach a constant state as the radial coordinate r tends to infinity to be compatible with periodic boundary conditions.
To obtain a steady analytical solution, we first write an equivalent PDE system in the radial direction.For this purpose, the pressure and the velocity relaxation are neglected, the gravity is set to zero and system (4) is rewritten in cylindrical coordinates and assuming no variations in the z−direction (∂/∂z = 0) and considering rotational symmetry (∂/∂θ = 0).The resulting system in radial direction reads where the constraint ∇ × w = 0 has been used in the last two equations.Since we are looking for a vortex-type solution, the radial velocities vanish, that is, we set u r = u r 1 = u r 2 = w r = 0. Also, we are interested in a stationary solution, hence ∂t = 0.With these assumptions, the system (25) reduces to ∂ ∂r With some simple algebra in (26) we get where k is a constant.Now, we prescribe radial profiles for α1 and pi as and then the densities ρi and the velocities u iθ result as , u iθ = 2 3/14 e 1−r 2 r 2 4 − e 1−r 2 5/7 .
With this, we have computed an exact stationary and rotationally symmetric solution of the PDE system (25) that approaches a constant state as r → ∞ to be compatible with periodic boundary conditions.Then, using the principle of Galilean invariance, we can make the test unsteady if we add a uniform velocity field to this solution.After one advection period through a periodic computational domain, the exact solution will be given by the initial condition and we can perform the convergence test.To analyze the convergence order, we compute the solution with the proposed method, using different orders for the DG scheme, and compare it with the exact solution derived above.The computational domain is Ω = [−10, 10] 2 with final simulation time t = 1 and periodic boundary conditions everywhere.Different polynomial approximation degrees have been considered for the DG scheme.The L 2 errors and the corresponding numerical convergence rates for N = 2, 3, 4, 5, are given in Table 1, showing the expected order of convergence.
For analyzing the convergence order using the unsteady solution, we make use of the principle of Galilean invariance of Newtonian mechanics.We add a constant uniform velocity field ū1 = ū2 = 4 to both phases.The computational domain is the same as before Ω = [−10, 10] × [−10, 10] with final simulation time t = 5 and periodic boundary conditions everywhere.The L 2 errors and the corresponding convergence rates for the different degrees N = 2, 3, 4, 5, are given in Table 2, finding the expected convergence order N + 1 of our high-order ADER-DG schemes.

1D Riemann problems
This section is devoted to studying the behavior of the proposed methodology in the presence of shocks.First, we solve one of the one-dimensional Riemann problems (RP) proposed in [31], where a shock in one phase appears inside a rarefaction of the other phase.This RP presents a discontinuity in x = 0 and has the left and right states shown in Table 3.Since the problem that we want to reproduce is a 1D problem, u 2 = u 3 = 0 and w 2 = w 3 = 0.The computational domain is Ω = [−1, 1] and has been discretized using a fourth-order ADER-DG scheme (N = 3) with a posteriori subcell limiter on a mesh with 1024 cells.The simulation is performed up to t = 0.25, and the CFL number is set to 0.25.Two ideal gases are considered for both phases with EOS (6), setting si = 0, γ1 = 1.4,and γ2 = 2, respectively.In Figure 1, the numerical results are shown together with the reference solution computed with a second-order MUSCL-Hancock scheme based on the Rusanov flux as approximate Riemann solver and using a mesh spacing of ∆x = 0.5 • 10 −4 (see [31] for further details).We observe an excellent agreement between our numerical solution and the reference solution.Looking at the density and velocity of each phase, we see that the rarefaction only affects the related phase, but that the interaction of the two rarefactions is observed in the mixing quantities.In addition, in the density of the first phase, we can see that to the right of the contact a rarefaction begins (which does not affect the second phase) until the shock occurs.Once the shock appears, the density jumps according to the jump conditions.Then on the right a plateau of the right state of the shock is observed, and then the rarefaction continues again, see [31] for a detailed discussion of this phenomenon.Table 1: Numerical convergence results for high-order DG schemes of polynomial approximation N = 2, 3, 4, 5 with a uniform Cartesian mesh of N x ×N y elements.The L 2 error norms and the corresponding orders of convergence of the variables α 1 , ρ 1 , ρ 2 , u 1 , and u 2 are computed at time t = 1. 4.2737 2.0532 • 10 −3

2D explosion problems
In this section, we solve the system for multi-phase flows in 2D in a circular computational domain with radius R = 1.The initial condition is given by where QL and QR are described in Table 4.As a reference solution, we will solve the following equivalent (non-conservative) PDE in radial direction with geometric reaction source terms where the parameter d is the number of spatial dimensions minus one.The 2D computations have been performed using a fourth-order (N = 3) ADER-DG scheme with a posteriori subcell limiter.The computational domain is Ω = [−1, 1] × [−1, 1] and has been discretized using a Cartesian mesh with 256 × 256 elements.Following (28), the left state of the RP has been taken as the inner state and the right state of the same RP as the outer state.The reference solution has been computed by solving (29) with 128000 cells using a second-order TVD finite volume method with the Rusanov flux.The simulation is performed up to t = 0.1 with two ideal gases, so for the two phases, the EOS is given by ( 6), with si = 0, γ1 = 1.4,and γ2 = 2, respectively.
Figures 2 and 3 show the numerical results of two circular explosion problems, with the initial conditions of Table 4 and a final time of t = 0.1 for the first one and t = 0.2 for the second one.The numerical solution obtained with the ADER-DG method is then compared with the radial reference solution, showing excellent agreement.Moreover, Figure 4 shows the limiter map of the second explosion problem.The values highlighted in blue are those DG elements where the limiter is not activated, and the red ones are the troubled zones where the a posteriori subcell FV limiter is activated.
The domain Ω has been discretized with a uniform Cartesian mesh with 256x128 cells, using an ADER-DG scheme with N = 3 and a posteriori subcell FV limiter.The simulation has been performed until a final time of t = 0.4, and a slip wall boundary condition is imposed on all boundaries.Following Section 2.1, an ideal gas is considered in Ω1, i.e., the EOS is given by ( 5), with parameters c01 = 1, γ1 = 1.4,ρ01 = 1, α = ε.The initial pressure profile is assumed hydrostatic, p = ρ01g(y − 2).The EOS for the liquid is a stiffened gas EOS given by (7) where c02 = 20, γ2 = 2, ρ02 = 1000, α = 1 − ε, and again a hydrostatic pressure profile p = ρ02g(y − 1) is imposed initially.The simulation was performed with ε = 0, i.e., initially the phase volume fractions are really set to zero and unity, respectively.To obtain the value of the primitive variable ρ k it is necessary to divide by α k , and in this simulation, there exist areas with α k = 0, and it is necessary to apply a filter that avoids division by zero.In this paper, the density variables are filtered as follows, , see also [67], and the filter parameter is set = 10 −12 .The numerical results have been compared with the solution of the reduced barotropic Baer-Nunziato model given in [68].Figure 5 shows the values obtained at time t = 0.4 calculated with the reduced Baer-Nunziato model in the upper plot, the solution calculated with the method proposed in this paper in the center, and a direct comparison between both models in the bottom plot, showing an excellent agreement between both models.Similar results have also been recently obtained with a novel Arbitrary-Lagrangian-Eulerian hybrid finite volume / finite element method applied to the incompressible Navier-Stokes equations on moving unstructured meshes, see [69].

Conclusion
In this paper, the barotropic version of the conservative SHTC model for compressible two-fluid flows of Romenski et al. has been solved for the first time using high-order ADER discontinuous Galerkin schemes in combination with an a posteriori subcell FV limiter.Since the model is only weakly hyperbolic in the general multidimensional case, two different methodologies have been presented to restore the strong hyperbolicity: i) a generalized Lagrangian multiplier (GLM) curl-cleaning approach and ii) the addition of the Godunov-Powell terms to symmetrize the system.We obtain a full set of linearly independent eigenvectors with both methodologies, proving that strong hyperbolicity can indeed be restored.
A high-order ADER discontinuous Galerkin finite element scheme with a posteriori subcell finite volume limiter has been used to deal with discontinuities and steep gradients in the solutions.To validate the model and the proposed method, a numerical convergence analysis has been carried out, and the high order of the method has been confirmed.For this purpose, we have constructed a new exact analytical and stationary equilibrium solution of the PDE system in cylindrical coordinates.Then, several Riemann problems in one and two dimensions have been simulated to show the behavior of the proposed methodology in the presence of shocks.First, a 1D Riemann problem where a shock in one phase appears inside the rarefaction of the other phase has been simulated.The results have been compared with those presented in [31], showing an excellent agreement.Then two 2D explosion problems were solved.Thanks to the radial symmetry of the problem, the obtained results have been compared with an equivalent 1D reference solution, showing the accuracy of the proposed methodology even in presence of sharp gradients in the solution.Finally, a dambreak test case has been considered, where the initial values of the volume fractions are set to α1 = 0 and α2 = 1.The numerical results are compared with those obtained for a reduced barotropic Baer-Nunziato-type model, showing an excellent agreement between both models.
As future work, we plan to extend our methodology to compressible multi-phase

Figure 2 :
Figure 2: 2D circular explosion problem for initial condition CE1 in 4 solved on a Cartesian mesh at time t = 0.1, in comparison with the radial reference solution.Top row: densities of each phase, ρ 1 and ρ 2 .Second row: mixture density ρ and α.Third row: velocities u 1 and u 2 .Bottom row: mixture velocity u (left) and relative velocity w = u 1 − u 2 (right).

Figure 3 :
Figure 3: 2D circular explosion problem for initial condition CE2 in 4 solved on a Cartesian mesh at time t = 0.2, compared with the radial reference solution.Top row: densities of each phase, ρ 1 and ρ 2 .Second row: mixture density ρ and α.Third row: velocities u 1 and u 2 .Bottom row: mixture velocity u (left) and relative velocity w = u 1 − u 2 (right).

Figure 4 :
Figure 4: Left: Limiter map of the explosion problem in 2D.The values in red mean that the limiter is activated.Right: 3D plot with the variable ρ 2 in the z−axis.

Figure 5 :
Figure 5: Dambreak problem at time t = 0.4.Top: reference solution, computed with a third-order ADER-WENO finite volume scheme on a very fine uniform Cartesian grid, solving the inviscid and barotropic reduced Baer-Nunziato model presented in [68].Center: Numerical solution, computed using an ADER-DG scheme with a posteriori sub-cell limiter, to solve the barotropic SHTC model proposed in this work.Bottom: Comparison of the free-surface profile obtained for both models.

Table 2 :
Numerical convergence rates for DG schemes of polynomial approximation N = 2, 3, 4, 5 with a uniform Cartesian mesh of N x × N y elements in the unsteady case.The L 2 error norms and the convergence orders of the variables α 1 , ρ 1 , ρ 2 , u 1 and u 2 , are computed at time t = 5, with ū1 = ū2 = 4.

Table 3 :
Left and right states of the RP1

Table 4 :
Left and right states of the circular explosion problems