An acoustic-transport splitting method for the barotropic Baer-Nunziato two-phase flow model

. This work focuses on the numerical approximation of the barotropic Baer-Nunziato two-phase ﬂow model. We propose a numerical scheme that relies on an operator splitting method corresponding to a separate treatment of the acoustic and the material transport phenomena. In the subsonic case, this also corresponds to a separate treatment of the fast and the slow propagation phenomena. This approach follows the lines of the implicit-explicit schemes developed in [8]. The operator splitting enable the use of time steps that are no longer constrained by the sound velocity thanks to an implicit treatment of the acoustic waves, while maintaining accuracy in the subsonic regime thanks to an explicit treatment of the material waves. In the present setting, a particular attention will be also given to the discretization of the non-conservative terms that ﬁgure in the two-phase model. We prove that the proposed numerical strategy is positivity preserving for the volume fractions and the partial masses. The scheme is tested against several one-dimensional test cases including ﬂows featuring vanishing phases.


Introduction
We are interested in the computation of compressible two-phase flows with a two-velocity two-pressure model derived after the Baer-Nunziato (BN) model [3].The (BN) model was initially dedicated to modeling a gas-solid flow for the simulation of detonation phenomena [3].Its range of applications has since been extended to a wider range of two-material flows like liquid-gas flows (see for example [11,14,18,20,22,23,25]).In the present work, we adopt a simplified version of the (BN) model that involves only mass and momentum equations equipped with a barotropic pressure law (as in [15,16]): the pressure p k of each component k = 1, 2 is assumed to only depend on the phasic density ρ k .A specific feature of this model that is inherited from the full (BN) model lies in the presence of non-conservative terms that drive the coupling between both fluids.These terms involve socalled interfacial velocity and interfacial pressure.The definition of these parameters require additional closure relations.This choice is not straightforward and has important consequences on the mathematical properties of the model, like the ability to define weak-solutions or to equip the model with an entropy evolution equation (see for example [19]).
From a numerical point of view, the barotropic (BN) model raises several issues.A first difficulty comes from the presence of non-conservative terms and how they may interact with conservative fluxes.A second difficulty occurs when the volume fraction of one of the component tends to zero in some regions of the computational domain.This situation is referred to as a vanishing phase and may cause severe problem as numerical methods often fail to compute the phasic quantities associated with the vanishing phase [15].
One may find in the literature several papers devoted to the numerical resolution of two-fluid two-pressure models and the question of how to discretize the non conservative terms.Most of them deal with the non barotropic case.The reader is referred to [1,2,26,29] for the numerical methods relying on time-explicit, exact or approximate Riemann solver and the references therein.We also mention some other finite volume techniques that have been used.In [19], the authors extend Rusanov's scheme and the VFRoe method to the context of non conservative systems.Other schemes rely on relaxation techniques (see for instance [1,15]).
The discretization of the full (BN) model has been investigated by many authors in the literature (see for example [1,6,[14][15][16]19,25,26,29]) but in most contributions, for stability reasons, the time steps ∆t is subject to a CFL condition that depends on the material velocity and the sound velocity of each material.For configurations where the sound velocities are much larger than the material velocities, this can lead to very small time steps although the acoustic waves are not driving phenomena in the flow.A CFL condition based on the most influent waves, the contact waves associated to the material velocities would be more adapted.The idea is then to propose a time-implicit treatment of the acoustic waves, in order to get rid of a too restrictive CFL condition, together with an explicit treatment of the contact waves in order to preserve accuracy.This strategy has already been exploited earlier within the framework of Euler equations and Shallow Water equations (see for instance [8,10,17]) using a Lagrange-Projection approach, but also for two-phase flows models (see [6,9,16,24,28,30]).
The paper is organized as follows.In Section 2, we present the set of partial differential equations (PDEs) of barotropic (BN) model, and we recall its main mathematical properties.In Section 3, we propose an operator splitting method for approximating the solutions of the model.In Section 4, we describe the numerical treatment of each step.In Section 5, we build explicit and semi-implicit numerical schemes based on this operator splitting method for the solution of the barotropic (BN) model.Finally, Section 6 is devoted to the numerical experiments, where one-dimensional test cases are presented.

The barotropic Baer-Nunziato model
Let us consider two compressible materials k = 1, 2, equipped with a barotropic Equation Of State (EOS) of the form ρ k → e k (ρ k ), where e k and ρ k respectively denote the specific barotropic potential energy and density of the fluid k.The pressure p k is defined by p k (ρ k ) = ρ 2 k de k /dρ k .We assume the sound velocity c k of the fluid k to be real-valued and defined by c We denote α k the volume fraction associated with the component k and impose that α 1 + α 2 = 1.If u k is the velocity of the fluid k = 1, 2, for one-dimensional problems the barotropic BN model reads: where u I and p I are the interfacial velocity and pressure for which one must provide closure laws.This system involves two evolution equations for both the partial mass and momentum of the fluid k = 1, 2 supplemented by an equation for the volume fraction The eigenstructure of (1) consists of five real eigenvalues The barotropic BN model is hyperbolic under the following conditions: When ( 1) is hyperbolic, one can easily check that similarly to the classical gas dynamics equations, the characteristic fields associated with the eigenvalues u k ± c k are genuinely nonlinear.
The choice of the closure law for u I has an important impact on the evolution of α k as it may not be possible to obtain a maximum principle for α k if the field associated with u I is not linearly degenerate.Moreover, the choice of both u I and p I affects the availability of weak solutions and a companion entropy evolution equation for the system.For these reasons, as in [15,16] we restrict our study to the following closure relations that issued from [12,19].This choice ensures that the field associated with u I is linearly degenerate, it also provides a complete set of jump relations that enables a full definition of weak solutions and an entropy evolution equation.As a consequence, smooth solutions of (1) with (2) also verify the following conservation equation Let us underline that for the numerical results presented in the last section, we will restrict our choice to the case χ = 0, µ = 1 so that (u I , p I ) = (u 2 , p 1 ).

Operator splitting acoustic-transport
Before going further on, we recall that this idea has already been used earlier within the framework of Euler equations (see for instance [17]) using a Lagrange-Projection approach, but also for the Baer-Nunziato model (see [6,16]).In [6], the authors propose an operator splitting with three steps for the Baer Nunziato model with energy equations.Both first steps correspond to the acoustic and transport parts of the Euler equations for each phase.The third step gathers the non conservative terms that couple both phases.In [16], the authors propose an operator splitting in two steps where the coupling terms are distributed in the two subsytems.In this case, the transport part arises a weakly hyperbolic system and a careful treatment has to be made for the numerical discretization.In this section, we introduce an operator splitting for the barotropic BN model in two steps where the two subsystems account for the coupling terms like [16] while preserving the hyperbolicity of each subsystem as in [6].
The first step corresponds to the propagation of acoustic waves: The second step considers the propagation of both material waves and the evolution of α k :

Properties of the acoustic sub-system
If we note τ k = ρ −1 k the specific volume, the acoustic system (8) takes the form Granted that c 1 = c 2 , the acoustic system (10) is strictly hyperbolic and the eigenstructure of the system is composed of five fields associated with the eigenvalues {±c 1 , 0, ±c 2 }.The waves associated with ±c k are genuinely nonlinear.The wave associated with 0 is linearly degenerate.
We propose further approximations for the acoustic system (8) by freezing the time dependence of ρ k in front of ρ k ∂ t (•) in (10) and by using a Suliciu-type relaxation approximation of (10) (see [7,13,27]).We proceed by introducing new independent variables π k and r k that respectively act as surrogate pressure and densities for k = 1, 2. If one considers smooth solutions of (8), the pressure p k verifies The variables π k are evolved according to their own partial differential equations, whose purpose is to implement a linearized version of (11).The new variables r k are associated with a stationary wave.Therefore, within the time interval t ∈ [t n , t n + ∆t[, we propose to consider the following relaxation system where ν is a relaxation parameter and a k is a constant approximation of the acoustic impedance of the fluid k.
In the limit regime ν → ∞, we formally retrieve the solution of (8).In order to provide a stable approximation of ( 8), the constants a k must be chosen in agreement with the Whitham subcharacteristic condition where the max is taken over all the density values ρ k in the solution of (12).We adopt the classic method that allows to reach the regime ν → ∞: at each time step, we enforce the equilibrium relations and solve (12) with ν = 0.For ν = 0, the relaxation system can take the compact form: with the unknown . Let us discuss a few properties of (14).Straightforward computations (see Appendix A) provide the following property on the characteristic fields of the relaxation system.Proposition 3.1.For any state vector (α 1 , W 1 , W 2 ), system (14) has the following characteristic speeds: where 0 is of multiplicity five.Moreover, all the characteristic fields are linearly degenerate and system ( 14) is hyperbolic.

Properties of the transport sub-system
We now consider the time evolution corresponding to the transport step.Starting from the output of the first step U , we want to compute the updated data U n+1 at time t n+1 by approximating the solution of The eigenstructure of ( 9) is derived thanks to straightforward computations.Proposition 3.2.For any state vector (α 1 , α 1 ρ 1 , α 1 ρ 1 u 1 , α 2 ρ 2 , α 2 ρ 2 u 2 ) such that ρ 1 > 0 and ρ 2 > 0, system (9) has the following characteristic speeds: Assuming that the interfacial velocity u I is defined by the closure law (2) then the characteristic field associated to the (simple) eigenvalue u I is linearly degenerate.The characteristic fields associated to the (double) eigenvalues u k are genuinely non linear and system ( 9) is weakly hyperbolic.Remark 3.3.Let us emphasize that the nature of the characteristic field associated with the eigenvalue u I in (9) depends on the chosen closure law for u I in the complete model (1).For example, assuming a different closure law as in [21] by setting: would lead to a genuinely non linear characteristic field associated to the eigenvalue u I in both ( 1) and ( 9).
Remark 3.4.The transport system (9) shares a similar structure with the system of pressureless gases (see for example [4,5]).However, in our case, its sole purpose is to yield a discretization means for the transport terms of the barotropic BN model (1).Let us underline that in the discretization strategy presented in Section 4.2, the propagation velocities (9) will be frozen at values u I and u k computed at the acoustic step.This will enable a discretization of ( 9) that verifies a local maximum principle for α k , α k ρ k and α k ρ k u k .Moreover, the overall discretization is consistent with the barotropic BN model (1).In this sense, possible well-posedness issues involved with the transport system (9) are not a genuine concern here and will not be investigated.

Discretization of the acoustic and transport sub-systems
In this section, we use the operator splitting method in order to derive an implicit-explicit numerical scheme, the aim being to approximate the solutions of the barotropic BN model ( 6).Let ∆t be the time step and ∆x the space step, which we assume here to be constant for simplicity in the notations.The space is partitioned into cells )∆x are the cell interfaces.At the discrete times t n = n∆t, the cell average of the solution of ( 6) is approximated on each cell C j by a constant value denoted by In the following two sections, we describe the discretization strategy associated with the operator splitting method in order to calculate the values U n+1 j j∈Z of the approximate solution at time t n+1 from those at time t n .Section 4.1 displays the numerical treatment of the Lagrangian step (8) while section 4.2 deals with the material transport step (9).

Treatment of the first step
We need to propose a discretization strategy for (14).The solution of a Riemann problem for (14) consists in six constant states separated by five contact discontinuities.We choose to build an approximate Riemann solver for the relaxation acoustic system (14).We will use in the sequel a discretization of the non conservative product that is consistent with the term π I ∂ x α k to derive an approximate Riemann solver for (14).Let ∆x L > 0, ∆x R > 0. We consider a piecewise initial data defined by: where the left and right states are defined by: Note that (π k ) L , (π k ) R , (r k ) L and (r k ) R are at equilibrium.Let us now build an approximate Riemann solver for the relaxed acoustic system (14).We look for a function (W k ) RP composed of six states separated by discontinuities as follows: The function (α 1 ) RP is defined as follows: The intermediate states in (16) are such that the following properties hold: (1) (W k ) RP is consistent in the integral sense with the barotropic Baer Nunziato model.More specifically in our context, if ∆t is such that where {π I ∂ x α k } is consistent with the non conservative term, in the sense: (2) We impose that (W k ) L and (W k ) L (resp.(W k ) R and (W k ) R ) verify the Rankine Hugoniot jump conditions accross − a k (r k ) L -wave (resp. (3) Similarly, across the discontinuity of velocity 0 we impose that: Relations ( 20) and ( 21) do not provide enough information to determine the intermediate states (W k ) L and (W k ) R .Indeed, they provide only seven independent relations while we need eight quantities, namely We choose to add another jump relation accross the stationary discontinuity of (W k ) RP , we impose where M k is a function to be specified.Relations (18), ( 20), ( 21) and ( 22) lead to: We now only need to determine M k such that the conditions 1), 2) and 3) are satisfied.The integral consistency requirement of Condition 1) imposes A simple mean to comply with ( 24) is to choose: where ) is a consistent approximation of π I in the sense that: At last, we choose: This yields that By construction, the approximate Riemann solver defined by ( 23) and ( 26) verifies the three conditions 1), 2) and 3).We end up with the following update at the acoustic step: where

Treatment of the second step
We now consider the numerical treatment of the time evolution corresponding to the second step.Starting from the output of the first step U j , we want to compute the updated data at time t n+1 , U n+1 j .Denoting φ k ∈ {m k , m k u k }, we use a standard time-explicit upwind discretization for the transport step by setting where Let us note that the transport update (28) equivalently reads: Let us note that the interface value of the velocity (u k ) j+1/2 coincides with the one proposed in the first step, which is actually crucial in order for the whole scheme to be conservative.

Two-step numerical method
In this section, we now give the details of the two-step process proposed in Section 3 for solving the barotropic Baer Nunziato model.Let us briefly recall that this two-step process is defined by • Update U n j to U j by approximating the solution of (8).• Update U j to U n+1 j by approximating the solution of (9).
We begin with a fully explicit discretization of the Baer Nunziato model, which means that both steps of the process are solved with a time-explicit procedure, and we will go on with a mixed implicit-explicit strategy for which the solutions of (8) are solved implicitly in time and the solutions of ( 9) are solved explicitly.The latter strategy allows to get rid of the strong CFL restriction coming from the acoustic waves in the subsonic regime and corresponds to the motivation of the present study.

Time-explicit discretization
Let us begin with the time-explicit discretization of the acoustic system (8), or equivalently (10).The acoustic update is achieved thanks to the proposed relaxation approximation and the corresponding approximate Riemann solver detailed in Section 4.1.More precisely, we propose to simply use a Godunov-type method based on this approximate Riemann solver.If we focus on the conservative variable: 27) yields the following formula for the acoustic update: Overall Discretization After injecting (30) in ( 29) one obtains the complete update procedure from t n to t n+1 for the conservative variables: The next statement gathers the main properties satisfied by our explicit in time and two-step algorithm.
(4) Under the Whitham subcharacteristic condition (13) and the CFL conditions the scheme preserves the maximum principle for the phase fractions: 0 < α k < 1 and positive values of the densities ρ k > 0. Proof.

Semi-implicit discretization
Let us now consider the last algorithm of this paper.It consists in considering a time-implicit scheme for the acoustic step and keeping unchanged the transport step.This strategy will allow us to obtain a stable algorithm under a CFL restriction based on the material velocity u k and not on the sound velocity c k .In order to derive a time-implicit scheme for the acoustic step, we follow the following standard approach where the numerical fluxes are now evaluated at time t , that gives here the same update formulas as in the explicit case but where the numerical fluxes now involve quantities at time t apart from the term consistent with {π I ∂ x α k }, which writes: Let us observe that we suggest here to keep on evaluating the interfacial pressure source term at time t n .It is interesting to see that it is equivalent to the following system written in characteristic variables: where the new variables These quantities are the Riemann invariants associated with the characteristic speeds ±a k τ k of the relaxation system (14).We firstly compute Remark 5.2.In the semi-implicit Lagrange projection, we need to solve a linear system to update the solution at time t .
where the block matrices are bidiagonal and the block matrices The implementation of this sparse matrix is made by the use of the Python library "scipy.sparse"and we solve the linear system (34) with a direct method from "scipy.sparse.spsolve".
Proposition 5.3.The implicit-explicit scheme (33)-( 28) satisfies the following: that our Lagrange-Projection schemes correctly capture the intermediate states even for this rather coarse mesh of 1000 cells.It appears that the contact discontinuity is captured more sharply by the explicit Lagrange projection method than by Rusanov scheme for which the numerical diffusion is larger.The two Lagrange projection schemes are comparable to the results given by the relaxation scheme introduced in [15] (see page 34).Nevertheless, we observe that the ImEx method generates an overshoot at the contact discontinuity where the phase fraction α k jumps.This oscillation is reduced on a finer mesh of 10000 cells in Figure 2 and does not generate an instability.To have a better understanding of this behavior at the contact discontinuity, we study a Riemann problem with a stationnary contact discontinuity in the last test case.The appearance of all the waves in this first test case and of their numerical diffusion can make the interpretation of the results difficult.Figure 3 also shows that the approximate solution computed thanks to the Lagrange Projection schemes converges towards the reference solution.The errors converge towards zero with a rate between ∆x 1/2 and ∆x.
The expected order of ∆x 1/2 could be otained by implementing the calculation on much more refined meshes in order to recover this asymptotic order of convergence.Table 1 gives for each test case the number of iterations needed to perform the computations.As expected, the gain is important when using the proposed implicit-explicit algorithm and the corresponding CFL restriction based on the material waves (instead of the acoustic waves as for the explicit scheme).We now consider a Riemann problem in which one of the two phases vanishes in one of the initial states.That means that the corresponding phase fraction α 1 or α 2 is equal to zero.This configuration poses a difficulty in the two-fluid model owing to its independent velocities.The singularity arises when one computes the absent phase velocity using the conservative variables u k = α k ρ k u k α k ρ k .For this kind of Riemann problem, the u I -contact separates a mixture region where the two phases coexist from a single phase region with the remaining phase.Assuming for instance that α L 1 = 1 and 0 < α R 1 < 1, the right state is a mixture of both phases while the left initial state is composed solely of phase 1.We consider the following initial data [16], V L = (1, 1.8, 0.747051068928543, 3.979765198025580, 0.6) if x < 0.5, V R = (0.4,2.081142099494683, 0.267119045902047, 5.173694757433254, 1.069067604724276) if x > 0.5.
The exact solution is composed of a (u 1 − c 1 )-shock wave in the left-hand side region where only phase 1 is present.This region is separated by a u I -contact discontinuity from the right-hand side region where the two phases are mixed.In this RHS region, the solution is composed of a (u 2 + c 2 )-rarefaction wave followed by a (u 1 + c 1 )-rarefaction wave.In practice, the numerical method requires values of α L 1 and α R 1 that lie strictly in the interval (0, 1).Therefore, in the numerical implementation, we take α L 1 = 1 − 10 −9 .The aim here is to give a qualitative comparison between the numerical approximation and the reference solution.Moreover, there is theoretically no need to specify left initial values for the phase 2 quantities since this phase is not present in the LHS region.For the sake of the numerical simulations however, one must provide such values.We choose to set ρ L 2 and u L 2 to the values on the right of the u I -contact discontinuity.As for the first test case, we can see in Figures 4,5 that for the same level of refinement, the explicit Lagrange projection scheme is more accurate than Rusanov scheme.That can be seen especially for phase 1.As regards the region where phase 2 does    not exist, we can see that the three numerical schemes develop some oscillations when it comes to divisions by small values of α 2 .This behavior is also observed with the relaxation scheme in [15] (see page 37).The ImEx method also generates an oscillation at the location where the phase fraction α 2 jumps and is more diffusive for the capture of the (u 1 + c 1 )-rarefaction wave (see Phase 1 density in Figure 4).The oscillation and numerical diffusion are lower on a finer mesh of 10000 cells in Figure 2 and ImEx method converges towards the reference solution.
Figure 6 also shows that the approximate solution computed thanks to the Lagrange Projection schemes converge towards the reference solution.The errors converge towards zero with a rate between ∆x 1/2 and ∆x.The expected order of ∆x 1/2 could be otained by implementing the calculation on much more refined meshes in order to recover this asymptotic order of convergence.

Test case 3: a pure contact discontinuity
In the last test case, we seek to analyze the behavior of the explicit and ImEx Lagrange Projection (LP) schemes and the Rusanov scheme for the capture of a pure stationnary contact discontinuity.This Riemann problem is built following the procedure explained in Appendix B. Here, in the exact solution, all the physical quantities are transported with the constant velocity u I = 0.The initial data is defined as V L = (0.8, 7.0710678118654755, 0.4448746176198241, 3.7907146169832258, 0) if x < 0.5, V R = (0.2, 3.1622776601683795, 3.979079545848609, 3.6840314986403864, 0) if x > 0.5. (38) Figures 7 and 8 gather the curves of (α 1 , ρ 1 , u 1 , ρ 2 , u 2 ) as functions of x and Figure 9 gathers the L 1 relative error compared to the initial state for the variables (α 1 , ρ 1 , u 1 , ρ 2 ), and the L 1 absolute error for u 2 (as |u 2 | = 0 everywhere).The approximation of a pure contact discontinuity represents a challenging test case.Indeed, for all schemes considered in the present article, including the Rusanov scheme several pathological phenomena can be observed.First, none of the numerical schemes succeeded in perfectly preserving the zero-velocity u I = u 2 profile.Such an error on u 2 also induces an error on ρ 2 as it its profile cannot be kept stationary.Nevertheless, u 2 seems to converge to zero (in L 1 ) with both Acoustic-Transport schemes, but not with the Rusanov scheme.
For the ρ 1 profile, both Acoustic-Transport schemes produce an overshoot at the location of the discontinuity.This results from the computation ρ 1 = (α 1 ρ 1 )/α 1 .This overshoot reduces when using a finer mesh.It is not observed with the Rusanov scheme that is more diffusive.Moreover, on the right-side of the contact discontinuity, ρ 1 seems to evolve as a regular solution over an important part of the computational domain and then finally reaches a plateau value on the right side of the domain.Even for finer meshes, the solution does not seem to converge towards the initial state when ∆x → 0. The same problem appears on the u 1 profile.Moreover, this very pathological behavior was also obtained with the relaxation scheme [15] (the simulations with the relaxation scheme [15] were performed by K. Saleh).Let us mention that the problem of capturing correct profiles for the coupling wave in BN models has been successfully investigated by [1] by the mean of an approximate Riemann solver.

Conclusions
We proposed in this work an extension of the Lagrange-Projection decomposition for the barotropic BN model by means of an operator splitting strategy.We studied both a full explicit time discretization and an implicit-explicit method that enables the use of time steps constrained by CFL conditions that do not involve the sound velocity of the fluids.Both the full-explicit and implicit-explicit methods are shown to ensure maximum principle for the volume fractions α k and positive values of the densities ρ k .Future work include entropy analysis through the Lagrange projection splitting, extension of the strategy to multi-dimensional problems and to the full Baer-Nunziato model with energy equations.The authors would like to thank K. Saleh for his advices and his help in performing the numerical simulation of the pure contact discontinuity test case with the relaxation scheme.

Figure 1 .
Figure 1.Test case 1: Space variations of the physical variables at the final time T = 0.14.Mesh size: 1000 cells.

( a )Figure 2 .
Figure 2. Test case 1: Space variations of the physical variables at the final time T = 0.14.Mesh size: 10000 cells.
(a) Explicit Lagrange Projection scheme (b) IMEX Lagrange Projection scheme

Figure 4 .
Figure 4. Test case 2: Space variations of the physical variables at the final time T = 0.1.Mesh size: 1000 cells.

( a )Figure 5 .
Figure 5. Test case 2: Space variations of the physical variables at the final time T = 0.1.Mesh size: 10000 cells.

Figure 7 .
Figure 7. Test case 3: Structure of the solution and space variations of the physical variables at the final time T = 0.05.Mesh size: 1000 cells.

( a )Figure 8 .
Figure 8. Test case 3: Structure of the solution and space variations of the physical variables at the final time T = 0.05.Mesh size: 10000 cells.

Table 1 .
Number of time-iterations for each test case with 1000 cells 6.2.Test case 2: vanishing phase