Hybrid Explicit Residual Distribution Scheme for Compressible Multiphase Flows

The aim of this work is the development of a fully explicit scheme in the framework of time dependent hyperbolic problems with strong interacting discontinuities to retain high order accuracy in the context of compressible multiphase flows. A new methodology is presented to compute compressible two-fluid problems applied to the five equation reduced model given in Kapila et al. (Physics of Fluids 2001). With respect to other contributions in that area, we investigate a method that provides mesh convergence to the exact solutions, where the studied non-conservative system is associated to consistent jump relations. The adopted scheme consists of a coupled predictor-corrector scheme, which follows the concept of residual distributions in Ricchiuto and Abgrall (J. Comp. Physics 2010), with a classical Glimm’s scheme (J. Sci. Stat. Comp. 1982) applied to the area where a shock is occurring. This numerical methodology can be easily extended to unstructured meshes. Test cases on a perfect gas for a two phase compressible flow on a Riemann problem have verified that the approximation converges to its exact solution. The results have been compared with the pure Glimm’s scheme and the expected exact solution, finding a good overlap.


Introduction
Due to the high interest in phenomena of shock propagation concerning areas such as engines, high explosives and energy, the domain of multiphase flows counts continuously increasing contributions. The main focus is on the ability of models to approximate the interfaces of mixed compressible flows. In particular, as outlined in [1], two main streams of approximations have been developed in the past years: one investigates flows using single pressure and velocity, while the other focuses on flows with multiple pressures and velocities. One drawback of treating multiple pressures and velocities is the fact that to each phase a velocity and pressure is associated, meaning that each phase has its own momentum and energy equations. This results in very large systems. Moreover, these systems are in the overall non-conservative and therefore it is not easy to determine explicit jump relations. One example of such non-equilibrium models is given by the seven equation system of Baer and Nunziato [2]. To simplify the model, it is possible to assume the shock propagation with a large interfacial area, meaning that single pressure and velocity can be assumed. This stiff mechanical relaxation has been widely investigated in [3,4,1] and many more. The following paper treats a five equation reduced model, as also many contributions as [1,4,5,6,7,8,9,10] do. This system originates from Kapila et al. [11], where they take the formal limit [12] of the Baer and Nunziato seven equation model when the relaxation parameters tend simultaneously to infinity, and thus permitting to assume a single pressure and velocity. The non-conservative character of the scheme due to the equation for the volume fraction keeps the question open how to provide compatible jump relations [13] to be able to obtain a sufficient accurate and robust approximation that converges to the exact solution. The possibility to simply consider conservative equations has been shown to induce large errors in the pressure and velocity computation as outlined with emphasis in [14,15,16]. Many methods have been studied in order to be able to treat the non-conservative system's jump relations (see e.g. [4]). Colella [17] presented the idea to apply Glimm's scheme [18] for continuous parts of the flow, and switch to Godunov scheme when large pressure jumps occur. Similarly, even though in contrast to the approach of Colella, Abgrall et al. [19] observed the ability of Glimm's scheme to capture a shock and the weakness of that scheme to propagate some numerical error in smooth areas, and coupled Glimm's scheme to the Roe scheme. In the following paper, we have coupled Glimm's scheme with an explicit residual distribution (ERD) scheme, and the switch between one scheme and the other is performed by means of a shock detector, capturing large momentum jumps. We examine the set of equations of Kapila et al. [11] in one dimension. The ERD scheme is a predictor-corrector method, following the concept of residual distributions in Ricchiuto and Abgrall [20]. The advantage of this formalism is that the time step is dictated by a CFL-like condition, contrarily to the original Baer-and Nunziato-like system. In general, a further advantage of the residual distributions is given by the absence of a Riemann solver, which results particularly practical for conservative set of equations as there is no need in the exact understanding of the Riemann solution. Moreover, the residual distribution formulation allows to see the component of the discretization as an error between two different approximations. In particular, the residual distribution is computed considering a Lax Friedrich's scheme [20]. For higher order accuracy in space, a limiter, designed similarly to a positive stream-wise invariant [21] with a filtering term taken from Harten and Yee [22], is applied. This paper is divided into four sections. In section 2 we give some definitions and assumptions for the problem set and recall the Kapila's system of equations. In section 3 we describe the explicit residual distribution scheme and in the next section how it is coupled with Glimm's scheme. Finally, we compare our approximation to the original Glimm's method with respect to an exact solution and present our conclusive remarks.

Mathematical model
The considered modeling equations are designed for hyperbolic compressible two-phase flows in one dimension. Among the multiphase models, Bear and Nunziato [2] gave a generic formulation of a seven equation model with heat and mass transfer, which consists of three conservation equations for each phase k that read and a transport equation of the volume fraction ∂α 1 ∂t + u I The terms α k , ρ k , u k , P k , E k represent respectively the volume fraction, the density, the velocity, the pressure and the energy of each phase. Quantities with the subscript I correspond to interface terms, while those denoted by tot are given by mixture terms. In particular, the following mixture relations hold for the volume fraction k α k = 1 and the total density k α k ρ k = ρ tot . Moreover, in (1)Ẏ represents the mass transfer term, while Q is associated to the heat transfer term. The coefficients µ and λ represent the dynamic compaction viscosity and the relaxation velocity parameter respectively [9].

The five equation reduced model
In this study we consider calorically perfect flows, s.t. terms associated to the heat transfer Q are neglected. Due to its high complexity we make the assumption to consider a very large interface and let µ and λ tend simultaneously to infinity. We retrieve a stiff mechanical relaxation of (1) as performed in [11] and [12]. The mass transferẎ is in the following neglected. The system shown in (1) is reduced to the well known five-equation model with the total energy given by the internal and kinetic energy E tot = e + 1 2 ρ tot u 2 , where the total internal energy is defined as e = k α k e k . In particular, the internal energy of phase k is defined as e k = ρ k k , and k corresponds to the specific internal energy of phase k. Further, we introduce a parameter K s.t. K = , where c 2 k = ∂P k ∂ρ k Entropy in k is the squared speed of sound associated to each phase. The total speed of sound is given by the Wood velocity, which is a weighted harmonic mean of the speed of sound of each phase Note that due to the relaxation of the parameters µ and λ, we retrieve u 1 = u 2 = u as well as P 1 = P 2 = P .
In order to close the system (2), we consider in the following two equations of state for a stiffened gas. The pressure is thus defined as P (ρ k , e k ) = P 1 (ρ 1 , e 1 ) = P 2 (ρ 2 , e 2 ) = (γ k − 1)e k − γ k P ∞,k for each phase k, with P ∞,k the reference pressure (see [9] for further details).

Jump relations for shock waves
The relation between the initial and final state across a shock (see also [4]) for (2) are given by the shock jump conditions of Rankine Hugoniot. In particular, given a left state V L and a right state V R we define the generic ∆V = V L − V R and have the relations We have = e ρtot , P = 1 2 (P L + P R ), τ = k Y k τ k , with the mass fraction of each phase defined as Y k = α k ρ k ρtot and the specific volume as τ k = 1 ρ k . The relative velocity is given by m = ρ tot (u − σ), where σ is the shock speed. These relations guarantee the conservation of the mixture.
3. An Explicit Residual Distribution (ERD) scheme for second order accuracy The adopted scheme is a predictor-corrector method, with the same spirit as the one of Ricchiuto et al. [20]. In the following we recall briefly the scheme.
Given an initial solution V 0 h at t 0 = 0, we consider the treated system to be in the form The first row of A reads A 1,s = [−u, − K u ρtot , − K u ρtot , K ρtot , 0] T and A r,s = 0, where s = 1, .., 5 and r = 2, .., 5.
We consider the time interval [0, T ] divided into N intervals s.t 0 = t 0 < t 1 <...< t n <...< t N = T and denote the numerical approximation as V n h = V h (t n ). Further, we divide the interval [t n , t n+1 ] in P sub-intervals, s.t. t n = t n,0 < t n,1 <...< t n,m <...< t n,P = t n+1 . Let ∆t n,m = t n,m+1 − t n,m and let V n,m h denote the approximation to V h (t n,m ). We discretise the spatial domain Ω with a grid denoted by Ω h , and denote by Q a generic element of the mesh and h the generic characteristic mesh size. Following the standard derivation of the Galerkin finite element method, we consider a spatial domain Ω divided in cells, with ∆x the spatial step between two neighboring nodes, and on this grid we consider the space of all the piecewise linear polynomials spanned by the basis functions ϕ i corresponding to node i ∈ Ω h . In particular these basis functions have to verify the condition of ϕ i (x j , y j ) = δ i,j , ∀i, j ∈ Ω h , and j∈Q ϕ j (x, y) = 1, ∀Q ∈ Ω h . The approximation of V will be denoted by V h and given by We define the residual to be r h = In general, β Q j determines the distribution strategy and, therefore, the applied residual distribution scheme. The reader may refer to [23] for further details on the generic residual distribution scheme.
We choose now as specific case the sub-intervals to be P = 2, s.t. the scheme is designed as a predictor-corrector method, and update the time step via an initialization at t 0 = 0 of ∆t = CFL · ∆x which is then updated by computing ∆t = CFL · min i ∆x i |u i |+c i . For the predictor we have As for the corrector, it contains a fluctuation computed on the V n+ 1 2 i and reads where |C i | is the area of the median dual cell C i obtained by joining the gravity centers of the cells with the midpoints of the edges meeting in i (see [20]).

The nodal fluctuation φ Q i is computed by a Lax-Friedrichs scheme
with with j = i, i + 1 and θ(V ) = 2 |u| + c tot where V assumes accordingly V n in case of the predictor, while V n+ 1 2 for the corrector. The quantities in (8) are defined as arithmetic averages with V = 1 2 j 1 2 V n j + V j and V j = 1 2 V n j + V j for j = i, i + 1. Vector ψ is the difference of fluxes occurring between nodes on the same cell given by with the flux F = [α 1 ρ 1 u, α 2 ρ 2 u, ρ tot u 2 + P, (E tot + P )u] T . Note that quantities such as u, K and c tot are defined the same way as V .
In general, (8) can also be expressed as where l s (φ j ) are the eigenforms corresponding to φ j for each eigenvector R s . Since by definition in (5), we have to distribute the residual to the different nodes which are part of a cell, we rewrite φ j as φ j = β j φ.
Recalling (10), we introduce a sort of positive stream-wise invariant (PSI scheme) with some stabilization, which has to enforce invariance along the characteristics as explained in [20,21].
For each eigenvector R s , with s = 1, .., 5 and reads for j = i, i + 1 the eigenform of (10) is , while η is a parameter that switches on and off the stabilization and accepts values between 0.05 ≤ η ≤ 0.125 according to [22] and, further, λ s the eigenvalues associated to R s .
The function ϕ (λ s ) provides an entropy fix and is designed according to the least dissipative function among the three point entropy satisfying TVD schemes of Harten-Yee [22], given by In particular we refer the reader to [20] for the proof of second order accuracy of the scheme applied to a different problem.

A Hybrid Glimm-Explicit Residual Distribution (GERD) scheme 4.1. Glimm's scheme
Glimm's scheme [18,24] is implemented according to the description in Colella [17] where the wave propagation is modeled by means of solving a Riemann problem, with consistent jump relations as shown in Eq. (3) (see also [4]). The approximation at each step is obtained by sampling the solution with the van der Corput random generator.

Coupled scheme
The ERD scheme, as also shown in Figure 4, approximates accurately the solution everywhere, besides the shock front, where the numerical dissipation affects the solution. On the other side, as shown in the limit test case in Figure 1-4 and as also outlined by many authors in this field (see [17], [19]), Glimm's scheme works perfectly for shock fronts, whereas gets some numerical oscillations in smooth regions. Therefore, we couple both schemes, so that the approximation is computed with the ERD scheme everywhere, with exception to the shock front. The criterion with which the switch is chosen is inspired partly by [17], partly by [19], but is intrinsically different and reads The idea is to check when the momentum becomes strictly greater than zero, since in that case a shock front is occurring. The solution n is first approximated by the ERD scheme to get , and then (13) is evaluated. In case θ n i > 0 Glimm's scheme is applied.

Numerical results
To validate the presented scheme we consider a typical very severe test case of a shock tube problem with two different mixtures and very high pressure jumps across the interface. This benchmark has been chosen due to its popularity among many authors as in [1,7,4,3], etc.
We are dealing with a Riemann problem, where a tube of 1 meter length is divided into two chambers by an interface at x = 0.6m and initial conditions as outlined in Table 1.
For the solution V n+1 i , Glimm's scheme is applied in case θ n i > 0, else the approximation is given by the ERD scheme. On the left and on the right of the interface a mixture of epoxy and spinel is set up, with on the left and right u = u 1 = u 2 = 0 [ m s ], while on the left a pressure P L = P 1L = P 2L = 2 · 10 11 [P a] and on the right P R = P 1R = P 2R = 1 · 10 5 [P a]. Further, we choose to display the solution at t = 29µs.   We consider the expected solution to be given by Glimm's scheme for a very fine mesh with ∆x = 0.0001m. This expected solution is compared on a grid with ∆x = 0.002m with the pure Glimm's scheme as described in [17] and as briefly recalled in this paper. The pure ERD scheme and the GERD scheme are then compared on the same coarse mesh as Glimm. Finally we display the comparison between the pure Glimm's scheme and the GERD one on a fine mesh with ∆x = 0.0001m. In these simulations the ERD and the GERD scheme have the parameter η = 0.05. We consider the CFL=0.45. In Figure 1 the total density, velocity and pressure as well as the volume fraction of phase 1 with respect to the position in the shock tube for ∆x = 0.002m are shown. The pure Glimm and the GERD scheme show in the overall a good overlap also w.r.t. the expected solution. Clearly, also from the details of Figure 1 displayed in Figure 2 and Figure 3, it is possible to observe that indeed the pure Glimm results to have some oscillations in the smooth area after x = 0.3m, while captures sharply the shock front even with a coarse spatial discretisation. The GERD scheme is, on the contrary, able to approximate without oscillations the expected solution and the error is just due to the coarse mesh. Indeed, following this observation, in Figure 5 it is possible to see for the volume fraction of phase 1 that the GERD method converges to the expected solution.
Finally we see that there are no delays or time displacements and the presented method reveals itself as very robust, being able to capture the correct approximation for extremely large jump conditions. As a last comparison, we show in Figure 4 the pure Glimm and ERD scheme with ∆x = 0.002m and compare them. While the approximation along the shock tube, away from the shock front, results to be well approximated by the ERD, and more subject to oscillations for the pure Glimm, we observe that Glimm is able to approximate very well the shock, while ERD scheme fails due to the fact that arithmetic averages are being considered. For the other quantities, as the pressure, displayed in Figure 4 (d)-(f), we see that there is not a relevant error between the Glimm and the ERD scheme. This confirms our initial statement that the non-conservative equation of the volume fraction represents a difficulty for shock capturing.

Conclusions
We have presented an ERD method for the approximation of a compressible, one dimensional hyperbolic multiphase problem given by the Kapila et al. five equation reduced model.
Since the system is non-conservative, arithmetic averages of the volume fraction on the shock front do not make sense and result in a local numerical dissipation. On the other hand, Glimm's scheme permits to compute the exact solution of Riemann problems, but results to be less accurate then the ERD method where the solution is not displaying a shock. The idea of the paper has been thus to investigate the possible advantages of a coupled Glimm's scheme with the ERD scheme. To our knowledge the first attempt to check a hybrid scheme with Glimm and Roe on the five equation model has been performed only by Abgrall et al. in [19]. To check the robustness of the scheme a very severe test case has been considered, showing good overlaps with the expected solution. With respect to other contributions in that area, we are, moreover, currently extending the methodology to that particular non-conservative system, using a specific Roe-type averaging that respects all the invariants of the original system.