A pressure-based method for weakly compressible two-phase flows under a Baer-Nunziato type model with generic equations of state and pressure and velocity disequilibrium

Within the framework of diffuse interface methods, we derive a pressure-based Baer-Nunziato type model well-suited to weakly compressible multiphase flows. The model can easily deal with different equation of states and it includes relaxation terms characterized by user-defined finite parameters, which drive the pressure and velocity of each phase toward the equilibrium. There is no clear notion of speed of sound, and thus, most of the classical low Mach approximation cannot easily be cast in this context. The proposed solution strategy consists of two operators: a semi-implicit finite-volume solver for the hyperbolic part and an ODE integrator for the relaxation processes. Being the acoustic terms in the hyperbolic part integrated implicitly, the stability condition on the time step is lessened. The discretization of non-conservative terms involving the gradient of the volume fraction fulfills by construction the non-disturbance condition on pressure and velocity to avoid oscillations across the multimaterial interfaces. The developed simulation tool is validated through one-dimensional simulations of shock-tube and Riemann-problems, involving water-aluminum and water-air mixtures, vapor-liquid mixture of water and of carbon dioxide, and almost pure flows. The numerical results match analytical and reference ones, except some expected discrepancies across shocks, which however remain acceptable (errors within some percentage points). All tests were performed with acoustic CFL numbers greater than one, and no stability issues arose, even for CFL greater than 10. The effects of different values of relaxation parameters and of different amount equations of state -- stiffened gas and Peng-Robinson -- were investigated.

• A general thermodynamic description is assumed during the derivation of the governing equations and the numerical method, which are thus valid for different EOSs, such as stiffened gas, Peng-Robinson, and more complex, multi-parameter EOSs.
• Two different pressure variables are defined according to the scaling proposed by Bijl and Wesseling 38 , so that the acoustics is filtered out from the model.
• Staggered grids are used to prevent stability issues related to the checker-board problem.
• A robust, but scheme-dependent, discretization of the non-conservative terms has been derived following the pressure non-disturbance condition 69 , to avoid spurious oscillations across multi-material interfaces.
These choices are explained and justified better in the next sections.

Paper structure
The next section concisely reviews some key features of pressure-based formulations proposed for single-phase flows, which are important to support the modeling choices made in this work. Section 3 begins with the presentation of the underlying BNtype model, continues with the derivation of the dimensionless pressure-based model, and ends with a short digression about the thermodynamic models. Section 4 presents the numerical method developed to solve the resulting BN-type model and it is split in four subsections: Sec. 4.1 introduces the semi-implicit time integration of the hyperbolic operator, Sec. 4.2 defines the organization of the variables over the grids, Sec. 4.3 details the spatial discretization of each equation of the hyperbolic part of the model, and Sec. 4.4 explains how the relaxation terms are treated. To have an organized framework, the results are presented in three different sections. First, Sec. 5 concerns the verification of the proposed numerical method for single-phase flows, and compares some possible alternatives in the solution strategy. Section 6 moves to two-phase simulations but still on a verification level, as it presents some water-air problems without relaxation, to validate the behavior of the hyperbolic operator. Finally, Sec. 7 presents the results of the complete model. Sections 7.1 and 7.2 give also an illustration of the effects of different values of finite relaxation parameters, Sec. 7.4 considers almost pure fluids, and Sec. 7.5 compares the results obtained with two different thermodynamic models. Lastly, in Sec. 8 we draw the conclusions of our works and we discuss future development and potential exploitation. The manuscript includes also A and B, which contain some passages omitted in the derivation of the model and the numerical discretization for reason of space.
B. Re and R. Abgrall

2 SOME KEY FEATURES OF STANDARD PRESSURE-BASED APPROACHES FOR SINGLE PHASE FLOWS
Many researchers have proposed diverse pressure-based formulations for the Euler equations, able to address the challenges of the low-Mach limits. These studies serve as a precious basis for our work, in which we attempt to blend together some key features of these previous studies and to extend them to multiphase flows. For this reason, before explaining our work, we briefly review here, without claiming to be exhaustive, some fundamental concepts widely used in numerical methods for weakly compressible single-phase flows.
One of the most challenging aspects of low Mach flows is that the governing equations change their character: the system of equations of the compressible gas-dynamics is purely hyperbolic, while its incompressible counterpart has a mixed hyperbolicelliptic character with infinite propagation speed. As a consequence, pressure and density are weakly coupled, so the problem of retrieving the pressure from the density becomes ill-conditioned. This explains, at least partially, the misbehavior of standard density-based methods at low Mach 33 and has motivated the widespread of pressure-based formulations. These latter strategies reflect, in general, the weak pressure-density coupling by solving the governing equations in a segregate approach: the velocity is first predicted using a pressure approximation in the momentum equation, then a correction step is carried out to update the pressure and, finally, the velocity is corrected for the new pressure 38,70,71 . Segregate solution strategies prompt the use of staggered spatial discretization where thermodynamic variables are stored at cell centers, while velocity variables are stored at cell faces 44,72,73 . Contrary to co-located formulations, staggered formulations filter spurious pressure modes providing an improved stability, at a similar level of efficiency and conservation properties 74 .
Asymptotic analyses have suggested the use of multiple pressure variables, which account for the different physical roles played by the different orders of the pressure in the low Mach limit 46,75 . Performing a single time scale/multiple space scale asymptotic analysis of the Euler equations, in which the pressure for small Mach numbers ( ) is expressed as Klein 75 showed that a scheme for low Mach flows should take into consideration at least two pressure variables: the leading order (0) which plays the role of the thermodynamic variable, and the second-order term (2) which is the "standard pressure" that accounts for local force balancing and, for → 0, satisfies the Poisson equation. Instead, the first-order term (1) is associated with long wave acoustics and it should be taken into account when pressure waves of order () are important. The pressure decomposition (1) allows the compressible equations to converge toward the correct limit, namely to the solution of the incompressible ones for → 0 46,71 . Standard numerical methods for compressible flows use on a non-dimensionalization based on a single reference velocity (e.g. computed from a set of reference pressure, density, and length), which introduces, in the low Mach limit, a singularity in the momentum equation, due to the presence of the term 1 2 in front of the pressure gradient 71 . The pressure decomposition (1) cures this singularity. As a final remark, this asymptotic analysis highlights that the divergence condition for incompressible flows, that is ⋅ = 0 (where is the flow velocity), results from the energy equation, not from the continuity equation, which, in the zero Mach limit, simply describes the advection of density fluctuations.
The idea of multiple pressure variables can be implemented in several ways 38,46,76 . Here, we follow the strategy proposed by Bijl and Wesseling 38,71 , who defined the pressure scaling as where is the density, is a scalar velocity, the tilde indicates dimensional variables, and the subscript r denotes reference quantities. The scaling (2) is characterized by the parameter r (reference Mach number), defined as which expresses the overall compressibility of the flow field. This strategy does not take into account the first order pressure (1) . In the low-Mach limit, the system of differential equations describing the evolution of the flow becomes stiff. Hence, the explicit time stepping schemes, routinely used for highly non-linear compressible flows, become inefficient 46,54 , because the CFL (Courant-Friedrichs-Lewy) condition imposes a severe limitation of the maximum admissible time step. To circumvents the most stringent time step limitation, the acoustic terms should be integrated implicitly, while the convective and diffusive terms can be treated explicitly, since they impose only a mild stability limitation of the time step based on the flow velocity. This strategy is called semi-implicit time integration 75,77,78 , and it is a common feature of pressure-based schemes, shared especially by the schemes that are able to represent all-Mach numbers-i.e., from very small to Mach numbers of order one. The concept can been easily implemented in a semi-discrete, fractional step projection method, in which the equations to be solved sequentially contain both implicit and explicit terms 52,55,57,70,76,79 . This is the strategy adopted in the present work, but, alternatively, a similar idea can be enforced by splitting the fluxes in two parts, advective and non-advective 47,54,80 . Both strategies can be used to derive asymptotic preserving schemes 49,55,81 .

A BAER AND NUNZIATO TYPE MODEL FOR NON-EQUILIBRIUM MULTIPHASE FLOWS AT LOW MACH NUMBER
In this section, we derive the set of equations that is the basis of the proposed numerical method, explaining the starting point and the modeling choices. We start presenting the BN-type model and the relevant notation, then we derive the pressure formulation and we apply the scaling (2).

The Baer and Nunziato type model
The non-equilibrium multiphase model derived by Saurel and Abgrall in 11 assumes that each phase is compressible and evolves with its own pressure, temperature, and velocity. The model does not consider heat or mass transfer and it tends to the Euler equations far from the interfaces. The system of 7 governing equations can be written in the following compact non-conservative form: where is the vector of evolution variables, is the flux function (the pure conservative part), contains the non-conservative part, and are vectors of source terms modeling, respectively, pressure and velocity relaxation. According to the standard notation, the variables are: the volume fraction (note that 1 + 2 = 1), the density, the velocity, = the momentum, the pressure, and the total energy, defined as = + 1 2 2 , where is the internal energy. The numerical subscript of each variable denotes the phase to which it refers 1 . The pressure I and the velocity I model the average interface values over the two-phase control volume, and they are estimated as while and are relaxation parameters that express, respectively, how fast the pressure and velocity equilibrium is reached 11,14 . These variables depends on the nature of each fluid as well as on the topology of the multiphase flow. For this reason, in this work, they are user-definite finite and positive parameters. All variables in Eqs. (4) and (5) are dimensional. Even if presented for two phases, this model can be adapted to three or more phases, provided a definition of interface and relaxation terms is given. Moreover, considering 1 = 1 and 2 = 0, the model simplifies to the classical Euler equations for single-phase flows. Furthermore, if we sum the equations per variables, we have the Euler equations for the mixture conservative variables, that is for the mixture density, i.e.,̄ = 1 + 2 , the mixture momentum, and the mixture total energy. Thermodynamic models are required to close the model. Each component obeys its own EOS as a pure material, so for each fluid or phase, we consider a generic EOS in the shape = ( , ), where is the internal energy per unit of volume, namely = with the specific internal energy. We remind here some thermodynamic definitions and relations which are of interest 7 in the next sections. First, we introduce the following thermodynamic derivatives Accordingly, the definition of the speed of sound reads for each fluid Definitions (6) and (7) are valid for each fluid separately, although we have omitted the subscript denoting the phase to simplify the notation. For later convenience, we define also, for each phase, an interface speed of sound as

Pressure-based formulation
To formulate a pressure-based BN-type model, we need to derive an equation for the pressure evolution from the conservative form (4). Here, we describe only the main steps and the results, while the step-by-step derivation is given in A. The first step consists in expressing the total energy in terms of pressure, density, momentum, and energy. Given an EOS in the form = ( , ), we can express the partial derivatives of with respect to = { , } as Thus, we insert this into the energy equation for phase = {1, 2} in (4), which can be re-written as 2 where we have introduced the operator Δ which takes the difference between the phase and the opposite one, i.e., Δ 1 = 1 − 2 and Δ 2 = 2 − 1 .
As detailed in A, we replace the temporal derivatives of , , and according to the respective equations in (4). Then, re-arranging the terms and recalling the definitions of the speed of sound (7) and interface speed of sound (8), we write the pressure formulation of the BN-type model (4) as

Dimensionless pressure-based BN-type model
In this section, we proceed to make the system of governing equations dimensionless, according to the pressure scaling (2) proposed by Bijl and Wesseling 38 . For the sake of clarity, we re-write here the volume fraction, density and momentum equation of system (4), along with the pressure equation (10), for one phase only, highlighting the dimensional variables with a tilde. We   (11)- (14). In particular, we show the combination of reference quantities (̃ r ,̃ r ,̃ r , and̃ r ) required to express each dimensional variable in terms of its dimensional counterpart. The first two columns report the variables that are not affected by the pressure, whose scaling is standard. The last column refers to the thermodynamic variables that require particular care because of the presence of the pressure in their definition.̃ remind that the volume fraction is, by definition, a dimensionless variable.
The scaling procedure requires the definition of the set of (dimensional) reference variables. The first entries in this set are: a densitỹ r , a length̃ r , and a velocitỹ r . Conventionally, we define dimensionless density, length, and velocity as Combinations of these three reference variables are sufficient to make dimensionless all the variables in Eqs. (11)- (14), as shown in Tab. 1. However, as anticipated in Sec. 2, we adopt a special scaling for the pressure to filter out the long-wave acoustics and to cure the singularity in the momentum equations in the zero Mach limit. Indeed, we introduce also a pressure reference variablẽ r , and we define the dimensionless pressure as =̃ Consequently, the pressure derivatives defined in (6) are scaled as More care is required for the speed of sound, which depends explicitly on the pressure. We want to preserve the definitions (7) and (8) where we have introduced the reference Mach number defined in (3). A similar expression is found for the interface speed of sound, which is reported in Tab. 1. As we show in the next paragraph, the additional term 1 2 r plays a fundamental role in the scaling of the pressure equation (10).
Following the definitions given above and in Tab. 1, we express all variables in Eqs. (11)-(14) in terms of their dimensionless counterpart. By using a verbose notation to show all substitutions, we obtaiñ Clearly, the previous equations can be simplified. Noting that in Eq. (20) the two terms involving̃ r cancel out, we can immediately simplify Eqs. (18)-(21) by deleting the factors comprising̃ r ,̃ r , and̃ r . Then, we multiply Eq. (21) by 2 r and we re-arrange the terms, gathering those with this factor together. In summary, the final system of equations expressing the pressure-based formulation of the BN-type model defined in Eq. (4) reads We highlight that the adopted pressure scaling, by means of the additional term proportional to in the definition of 2 and 2 I, , is directly responsible for the peculiar expression of Eq. (25), in which terms proportional to 0 r and 2 r coexist. The fundamental benefit of this choice is expressed by the following Remark.

Thermodynamic models used in this work
The system of governing equations presented above and the numerical method described in the next section are derived without any specific assumption on the thermodynamic models, as long as the EOS of fluid can be expressed as = ( , ). Since this requirement is pretty easy to meet, most of the EOSs used for academic and industrial purposes can be adopted to model the behavior of each component within the proposed pressure-based BN-type model. To introduce the nomenclature used in the following sections, we briefly introduce here the models used in for the simulations presented in result sections 5, 6 and 7, namely the stiffened gas model 82,83 and the polytropic Peng-Robinson EOS 84 . A complete thermodynamic model of a pure fluid at equilibrium can be obtained from two independent EOSs, the thermal and caloric one. For the stiffened gas, their dimensional expressions read (omitting the tilde to lighten the notation) where is the temperature and = is the specific heat capacity at constant volume, which is constant in the stiffened gas approximation. The parameters (ratio of specific heat capacity), ∞ , and depend on the material and can be determined by fitting experimental data, e.g., the saturation curve 22,85 . The expressions of other thermodynamic variables can be found in 86,87 .
Stiffened gas model can be considered an extension of the polytropic ideal gas (which is recovered when ∞ = = 0) able to take into consideration the repulsive effects present in all states of matter (modeled by the term − −1 ) and the cohesive forces typical of liquid and solid states (thanks to the term ∞ ) 85 . This capability together with its simplicity accounts for its wide use in the research activities focused on the development of models and numerical tools for two-phase flows, as the present one. When the focus is the study of complex two-phase flow behavior, e.g., for the investigation of water cavitation problems or in process simulation of renewable energy technologies, more accurate EOSs may be required. An answer to this demand may come from cubic EOSs, which are widely used also in industrial applications because they combine a decent accuracy with computational efficiency. A popular instance in this class is the Peng-Robinson 84 EOS, which is used in this work to model liquid and vapor CO 2 in Sec. 7.5. The expression of thermal and caloric EOSs for this model can be found in 86,88 .
All thermodynamic variables are made dimensionless following the scaling rules defined in Tab. 1 and a standard scaling for the temperature according to a reference dimensional valuẽ r , that is̃ = ̃ r .

NUMERICAL METHOD
For the aim of this work, a numerical method that is first order accurate, both in time and in space, is considered. The system of governing equations (22)- (25) is solved according to the Strang splitting approach, as in 11,25,26,89 . Hence, given the solution at a initial time , the solution +1 after a time interval Δ is obtained by the sequence of operators where hyp is the operator that solves the hyperbolic part of the system over a time step Δ , while relax is the relaxation operator that solves the system of ordinary differential equations (ODEs) considering only the relaxation terms for the velocity and the pressure. We describe hyp in the sec. 4.1-4.3, while relax in sec. 4.4.

Temporal discretization of the hyperbolic operator
For the numerical discretization of the hyperbolic operator hyp , we start from the time integration, keeping the spatial derivatives continuous. To mitigate the time step restriction imposed by the CFL constraint, we use a semi-implicit temporal discretization where the acoustic effects are treated implicitly. This requires to integrate implicitly the pressure gradient in the momentum equations and the divergence of the velocity in the pressure equations. To easily handle the first task, we adopt a time splitting in which the momenta (and the velocities) are first estimated by treating explicitly the pressure gradient, then they are corrected according to the updated pressure values, as done for instance in 47,54 . Moreover, at the end of the time step, we recompute the density with the current advection velocity, to have a better accuracy when the Mach number is particularly low. The semidiscretization of the governing equations per each phase reads where * , * = ( ) * * and * are the predicted density, momentum and velocity. The superscripts and + 1 indicate, as usual, variables at the previous time step, , and at the end of the hyperbolic operator, +1 . The double star in the momentum correction equation (33) highlight that it is related to the predicted density, i.e. * * = ( ) * +1 . We can interpret the previous set of equations also in the framework of multiple-pressure variables. Indeed, this semi-implicit discretization computes the convective and thermodynamic effects in a predictor step, composed by Eqs. (29)-(31), while the high-order pressure effects, that is the ones due to the term (2) in Eq. (1), are corrected implicitly through the Eqs. (32) and (33) 75 . Moreover, if we sum (31) and (33), we get the equation for the momentum * * with implicit pressure gradient, i.e., * * − The final momentum is computed after solving (34) as On the other hand, we use a different approach for the density equations (30) and (34): the results of the former one, i.e., * , are used only while solving (31)- (33); but, at the end of the time step, while solving (34), the densities +1 are computed starting from , discharging * . This re-computation allows the use of the most updated advection velocity, +1 , which is particularly important in flow problems close to the incompressibility limit, where the density equations simplify to transport equations. We compare the results obtained with and without density re-computation in the first numerical test, in Sec. 5.1.
A final remark concerns the divergence of the velocity in (32), which needs to be treated implicitly to overcome acoustic CFL limitations 54 . However, since Eqs. (29)-(34) are solved in a segregate approach in the order they appear, the value of the velocity +1 is not known while solving the pressure equation. Inspired by the use of the momentum equation to derive an implicit pressure equation 54,55 , we use Eq. (33) to approximate the value of +1 . Indeed, to a first approximation, the difference in the convective terms can be neglected 38 , so the final velocity can be approximated as As explain better in the following, this choice, and in particular the term +1 I = ∑ +1 , couples the pressure equations for all phases together. On the contrary, Eqs. (30), (31), (33), and (34), are solved per each phase independently.
Remark 4 (Alternative formulation). Considering the definition of * * and that ( ) * is already known, instead of the momentum correction (33), we could also correct directly the velocity by solving Spatial discretization and variables positioning. At the bottom, the computational domain Ω = [ 0 , ] is drawn, along with the position of the grid nodes and the boundary domains sketched with grey dashed lines. The upper part of the picture illustrates the primary and staggered grids. They are drawn with a vertical development and separately from the computational domain only for greater clarity, but all grids and cells here defined should be considered as one-dimensional. In the upper part of the picture, blue color refers to the primary cells  and to the quantities , , and , which are stored at their centers, represented by blue square marks. Green color refers to the staggered cells, at the center of which the kinetic variables and are stored; these cells are sketched by dashed lines and their centers are represented by green circle marks. In red, the boundary values for the primary grid complement the picture. To lighten the notation of variables, phase indication is omitted and the numerical subscripts refer simply to the spatial cells.

Variables positioning: primary and staggered grids
For the spatial discretization of the pressure-based BN-type model, we consider two finite-volume schemes based on staggered grids: one for the thermodynamic variables, and one for the kinematic variables. FIGURE 1 shows how the staggered grids are defined. We split the computational domain Ω = [ 0 , ] in intervals, defined by the equidistant grid nodes , with = 0, … , . Then, the grid for the finite-volume discretization of the thermodynamic quantities (hereafter called primary grid) is built by defining each cell  corresponding to the grid element [ −1 , ]. Conversely, the grid for the finite-volume discretization of the kinematic variables (hereafter called staggered grid) is built by centering each cell around the grid nodes , between the centroids of the adjacent grid elements (or the boundary). In summary, the cells on primary and staggered grids are defined as primary: As it appears from the given definitions, starting from grid nodes equally spaced by a distance Δ , all primary cells have the same size | |  | | = Δ , while the staggered cells have the same size | | | | = Δ only far from boundary. Indeed, the first and last staggered cells are half the size, i.e., | | 0 | | = | | | | = Δ ∕2. A cell-centered finite-volume discretization over the primary grid is used to solve the volume fraction, density, and pressure equations, that is Eqs. (29), (30), (34), and (32). So, the thermodynamic variables (sometimes called "scalar" in contrast to the kinematic, vectorial variables) are approximated over the cell  as On the other hand, the momentum and momentum update equations, (31) and (33), are discretized over the staggered grid. Using a finite-volume scheme, we define the cell value of the momentum as These are the finite-volume cell values, illustrated also in Fig. 1. However, it is often required to map variables from "their" grid to the other. In this case, we perform a weighted average, which, being the grid nodes equidistant, simply results in an arithmetic mean.
Remark 5 (Notation and mapping). To have a clear notation, we use the subscript for quantities over the primary grid, and for quantities over the staggered one. Accordingly, a thermodynamic variable with a subscript refers to its mapped value over the staggered cell , and vice versa. To clarify this point, consider the following example. The notation ( ) indicates the mapped density over the cell , computed as ( ) = [( ) + ( ) +1 ]∕2 for = . We can then use this mapped density to estimate the velocity in the cell as ( ) = ( ) ∕( ) .

Spatial discretization of the hyperbolic operator
Each hyperbolic differential equation in the model is integrated in time for the interval Δ = +1 − and in space over all cells ( or ). The spatial derivative of the convective fluxes is approximated through numerical evaluations of the fluxes at the cell interfaces. In particular, we use a first-order approximation based on the Rusanov flux, as in 11,19 . This choice is motivated by simplicity, as it avoids the complexities related to the solution of local Riemann problems with several waves [90][91][92] . The use of staggered grids makes the discretization of some specific terms easy and natural.
• The convective velocity to be used in the flux computation on the primary grid is directly the velocity defined over the staggered grid. For instance, the Rusanov flux for at the interface between  and  +1 is computed as • The pressure gradient in the momentum equation is readily approximated by a centered difference scheme, since the values of the pressure at the faces of the staggered cells are available: • The divergence of the velocity in the pressure equation is easily discretized through a centered difference scheme: A major complexity of the spatial discretization of Eqs. (29)-(34) concerns the presence of non-conservative terms involving the gradient of the volume fraction. This is a challenge common to all BN-type models, which include the term ( ) 1 that models the momentum and energy transfer among phases but prevents to write Eq. (4) in divergence form. This means that it is not possible to define weak solutions in the standard sense of distribution and to determine unique wave speeds. From a numerical point of view, these non-conservative products have to be integrated as source terms, rather than as fluxes. Since a naive discretization may introduce spurious oscillations across material interfaces between phases with different specific heat ratios, we seek a robust discretization of non-conservative terms involving the volume fraction gradient by explicitly enforcing that uniform velocity and pressure profiles are maintained 69 . Honestly, different strategies can be followed to integrate the nonconservative terms associated to the linearly degenerate fields, as, in particular, path conservative schemes 93 . However, this approach does not guarantee to always converge to the correct weak solution of non-conservative hyperbolic problems 94 . In addition, a primitive formulation of the governing equations, such as the pressure one here considered, facilitates preserving pressure equilibrium near material interfaces 36 . All in all, for weak discontinuities, as the ones considered in the framework of weakly compressible flows, any consistent and accurate enough method would be adequate to achieve a satisfactory solution 95 .

Volume fraction and density equations
Without any non-conservative term, the density equations (30) and (34) are easily discretized in space as where the expression for the Rusanov fluxes is given in Eq. (38), and the superscript ⋄ corresponds to * and to + 1 in the spatial discretization of (30) and (34), respectively. The discrete volume fraction equation is where ( +1 , I ) ≈ ∫  I +1 d is a suitable approximation of the non-conservative term. To define this operator, we follow the idea that starting from a uniform pressure and velocity, no variations in these variables should be generated 11,19,69 , see also Remark 3.
If we assume a uniform velocity field, e.g. ( ) = ( ) +1 = ( I ) = , the discrete mass equation reads where we have dropped the superscripts in the left hand side to lighten the notation. Let us consider now the special case when also the density field is uniform 19 . If ( ) = ( ) −1 = ( ) +1 , the mass equation reads If velocity and density are uniform, the density should remain constant, i.e., ( ) +1 = ( ) . So, in order to make Eq. (42) compatible with Eq. (44) in this specific case, we need that From this, we define the following non-conservative operator : wherê Rus We use the notation̂ Rus to highlight that that the resulting discretization of depends on the discretization for the convective flux in the mass equation, but, at the same time, the⋅ indicates that it is not a proper flux, as ( I ) is the mapping of the interface velocity over the primary cell  , not an interface velocity. This choice guarantees that = 0, if the volume fraction is uniform, as expected by the integration of I .

Momentum equations
The spatial discretization of the momentum equations (31) and (33) requires the integration of three terms: the convective flux, for which we adopt a Rusanov flux; the pressure gradient, discretized by the central finite difference defined in (39); and the non-conservative term, for which we define the operator ( +1 , I ) ≈ ∫ I +1 d exploiting the non-disturbance pressure and velocity condition, as explained in the following. Accordingly, the discrete equations of the predicted momentum and of the corrected momentum read where we have used the operator to identify the jump in the kinetic and pressure variables between the prediction and correction step. More precisely, The Rusanov fluxes are defined, as usual, as where . The same expression but with ( ) * * instead of ( ) * is used for B. Re and R. Abgrall

15
The discretization of the non-conservative term ( , I ) is derived, similarly to ( , I ) , by imposing the nondisturbance pressure and velocity constraint. The whole process is detailed in B. For conciseness, we report here only the final definition of the non-conservative operator : where ( I ) = 1 2 ( I ) + ( I ) +1 is the interface pressure mapped at the staggered cell . Finally, we highlight that Eqs. (46) and (47)

Velocity correction equation
If we consider the velocity correction equation (37), its discretization is straightforwardly derived from the spatially discrete equation of momentum correction (47) and it reads where the Rusanov fluxes are defined, similarly to (48), as

Pressure equation
We develop now the discrete version of the non-conservative pressure equation (32). First, we observe that, in the considered finite volume context, the thermodynamic variables and the volume fraction are constant within the primary cell, as in 38 . So, integrating Eq. (32) over a cell  , we can write where, to have a more compact expression, we have introduced the two coefficients (K ) = 2 r ( ) ( 2 ) + ( ) , (K ) = 2 r ( ) ( 2 I, ) + ( ) , which are known, because the variables ( 2 ) , ( 2 I, ) , and ( ) are computed using the thermodynamic state at cell  and at time . For instance, ( ) = ( ) , ( ) , according to definition (15).
The second step concerns the discretization of the first integral term, which, thanks to the product rule, is re-written as To approximate the first term in the previous expression, we define the following flux (similar to (38)) , while for the second one, we rely on the central approximation scheme for the divergence of the velocity given in (40). We obtain:

B. Re and R. Abgrall
A third aspect to be considered is the approximation of the non-conservative term involving the gradient of the volume fraction. Given the similarities with the non-conservative term in the volume fraction equation, we adopt the same operator defined in (45), but for the velocity jump. Thus, The remaining integral term in Eq. (52) is easily approximated by a central difference scheme, but it requires an expression for the velocities at the time step +1 . This latter is derived from the discretization of the velocity update, Eq. (51), discharging the differences in the convective terms. It reads In conclusion, the discrete version of the pressure equation is (cfr. Eqs. (36) and (53)). Recalling the definition of ( I ) and the mapping from the primary to the staggered, we have from which it appears evident the involvement of the pressure of both phases in the definition of the velocity ( ) +1 . Consequently, we need to solve the pressure equations (54) for both phase together, i.e., in a coupled way.

Boundary conditions
To impose boundary conditions, we distinguish between primary and staggered grid. For the primary grid, we use a standard method based on two ghost states defined outside the computational domain. With reference to Fig. 1, these states are denoted by subscripts 0 and + 1, on the left and right boundary, respectively, and are defined as According to the physical boundary condition we need to model, the value of the variables in +1 B mirrors the state of the adjacent internal cell ( 1 or  ), or it is directly imposed as boundary value (for the details about this selection process, see for instance 96 ). The boundary state +1 B is then used in the discrete equations (41), (42), (46), (54), and (47) to evaluate the fluxes, the non-conservative terms, and the central difference schemes at the boundary interfaces.
For the staggered grid, we use a different strategy, because the first and the last staggered cells ( 0 and ) are boundary cells. In addition, the velocities ( ) 0 and ( ) are already stored at the boundary interfaces (see Fig. 1). Thus, the momentum and velocity in these two cells are not computed by solving Eqs. (46) and (47), but they are computing according to the physical boundary condition. In particular, we distinguish two cases: if the boundary velocity ( ) B is known, its value is imposed; otherwise the velocity is extrapolated from the two closest internal cells. For instance, considering the left boundary: where ( ) 0 , ( ) 1 are the density values on the primary cells. This definition applies also to the implicit velocity in the pressure equations, i.e., while solving the equation (54) for ( ) +1 1 and ( ) ) the expressions for the velocity ( ) +1 0 and ( ) +1 are the ones given above, instead of (53).

Solution of implicit system
The implicit treatment of some terms in the discretization of the hyperbolic operator makes the equations coupled between adjacent cells. Indeed, the structure of mass, volume fraction, and momentum equations, e.g. (41), (42), (46), and (47), can be approximately represented as where is the unknown of a specific phase in the cell (over the primary or staggered grid), Δ is the cell volume, L and R are fluxes across the left and right cell face, R and L refer to the cell centered discretization and the term represents the discretization of the non-conservative terms (involving different values of +1 according to the equation we are considering). The superscripts and + 1 indicates, generally, known and unknown values, respectively. Obviously, not all right hand side terms are present in every equation, but there is always at least one term that generates the cross-coupling.
The previous expression can be further simplified as where Φ includes the fluxes or the non-conservative terms that are function of the unknowns themselves and R includes all the known terms. We use a first order Taylor expansion to approximate Φ as where +1 = +1 − . Since we use Rusanov fluxes, the derivatives of fluxes and the non-conservative term (required only in the volume fraction equation) can be easily computed analytically. Hence, for Eqs. (41), (42), (46), and (47), for each phase separately, we need to solve a set of equations in the form which is comprised of or + 1 equations, depending on whether we are considering the primary or the staggered grid. The resulting systems are linear and they can be written, in a compact form, as [ ] = , where [ ] is a tridiagonal matrix including the derivatives of Φ, is the vector of unknown, and is the known term. These systems are solved through the Generalized Minimal Residual (GMRES) algorithm provided by the PETSc library 97 .
Similar observations can be drawn also for the pressure equation (54), which however is solved for both phases together. In this case, the generalized term Φ includes also the unknown terms deriving from (53), that is Consequently, the first-order Taylor expansion involves six different unknowns, which can be organized in a vector in this order: … ,

Relaxation operator
According to the Strang splitting introduced in (28), the solution of the hyperbolic operator described in Secs. 4.1-4.3 provides a known set of variables that are used as initial data to solve the system of ODEs associated with the relaxation terms. For this reason, in this section, we re-define the notation to distinguish the intermediate solutions after the hyperbolic operator hyp , and the relaxation operator relax as follows In practice, in this subsection the superscript • denotes what in subsections 4.1 and 4.3 was denoted by + 1, and the superscript • refers to the variables computed during the relaxation processes. The relaxation operator relax plays a fundamental role in driving phasic velocities and pressures toward the equilibrium, close to interfaces. The characteristic time of these processes depend on many factors, as the fluids features and the multiphase flow topology. For instance, the parameter , which expresses the velocity of the pressure relaxation, may depend on the compressibility of the fluid and the parameter , which governs the rate of the velocity homogenization, may depend on fluid viscosity 11 .
In general, pressure and velocity relaxation are much faster than the dynamics associated to the wave propagation, to the point that they are something modeled as instantaneous phenomena, by assuming infinite and 11,25,27 . However, in this work, we use finite relaxation parameters, as in 2,98 , to allow wider modeling possibilities. Indeed, we could define the relaxation parameters in terms of the average interfacial area of bubbles 99 , or, if we had experimental data about different multiphase flow topologies, we could tune the relaxation parameters in our model to match the data.
Assuming a characteristic time much shorter than the one characterizing the hyperbolic operator, the ODE system associated to the relaxation operator is derive from the continuous governing equations (22)-(25) neglecting convective and transport terms. It reads This system is characterized by a high degree of stiffness, so we use the implicit Backward Euler scheme for the time integration.
Equations (60) give immediately +1 = • . If we use this result in (61) and we integrate in time, we have where the only unknowns are the velocities. The solution of this system, expressed in term of Δ , is which, as expected, gives Δ → 0 when → ∞, so that In the opposite case, for = 0, we have The remaining part in the ODE system comprises the volume fraction equation (59) and the two pressure equations (62). After the discretization of the time derivatives, these equations can be re-written as where K = 2 r 2 I, + , as in (52), and K = 2 r ( I − ). Reminding that the velocities • are given by (63), the last terms in (65) and (66) are known. For the discretization of the coefficient K , we approximate the thermodynamic variables and the interface pressure by using the values at the end of the hyperbolic operator. This choice is a simplifying assumption, which slightly mitigates the non-linearity of the system (64)- (66), and it is motivated by the absence of differences noted in 24 while solving the pressure relaxation system approximating the integral value of the interface pressure by • or • .
From a numerical point of view, the non-linear system (64)-(66) presents some unfavorable features, such as the simultaneous presence of very small and very large terms, which could cause a loss of accuracy, the stiffness and the non-linearity. To tackle these aspects, we rely also for the relaxation operator on the PETSc non-linear solver 97 and, in particular, on the trust-region Newton-based solver.

VERIFICATION FOR SINGLE-PHASE FLOWS
Since both the model and the numerical method proposed in this work are new, before focusing on two-phase simulations, we present in this section some single-phase tests, to verify the numerical method and, in particular, the low-Mach treatment and the core part of the hyperbolic operator. The governing equations (29)-(34) (in their fully discrete versions given in Sec. 4.3) are here solved only for one phase, but the numerical solution algorithm is kept unaltered. And by that, we mean that the volume

Low Mach Riemann problem for a perfect gas
We start with a Riemann problem test at particularly low Mach number, presented in 100 . The pipe is filled with a perfect gas, i.e. air with the parameters given in Tab. 2, at very low pressures and the left and right chambers features weak pressure and velocity jumps, according to the data given in the row lmAir in Tab. 3. The solution is represented by two rarefaction waves, plus a central contact discontinuity which moves at = 4.7 10 −3 m∕s. The Mach number is lower than 0.012 all over the domain. Figure 2 displays the results at the final time obtained with the standard formulation described by (29)-(34) considering only one phase. Six simulations are run imposing six different time steps Δ , defined as Δ = ∕ where is the final time and is the number of integration steps used to reach the final time . The simulations are labeled in the picture according to the number , which goes from 500 to 15 (from the smallest to the largest time step). As reported in the caption, some of them lead to an acoustic CFL greater than one, in short, CFL(| | + ) > 1. The contact discontinuity, which has moved by only one cell, is sharply represented, while the rarefaction waves are smeared because of the first-order accuracy of the Rusanov flux.
We use this test also to show the role played by the density correction and the velocity formulation. Figure 2 compares the results obtained with and without density re-computation. In particular, we compare three formulations: : solves only the mass equation (30)  From the density profile, we can notice that solving the mass equation only at the beginning of the time step, so that using the convective velocity , leads to some oscillations across the contact discontinuity, which are amplified if the CFL number increases. The old value of the velocity does not account for the pressure correction, which is responsible for enforcing the incompressibility condition (see Remark 1). On the other hand, from the velocity profile, we can notice that the computation of the density only at the end of the time produces slightly worse results than the standard formulation with density re-computation. Since the density equation (41) does not present numerical difficulties, e.g. it does not include non-conservative terms, and its computational effort is almost negligible with respect to the solution of the other equations, we adopt the re-computation as the standard formulation. A further open question in the development of the numerical method here proposed concerns the momentum or velocity correction, that is whether (33) can be substituted by (37). For this reason, we have re-run the simulations presented in Figures 2  and 3 with the velocity correction, where (33) is solved instead of (37). No notable differences are detected, and, in particular, the same conclusions about the density re-computation are drawn, as shown by Fig. 4.

Low Mach Riemann problem for a stiffened gas
In this section, we address the simulation of a water pipe flow under the stiffened gas model, proposed in 100 . The thermodynamic parameters are given in Tab. 2 and the initial data are reported in the row lmWater in Tab. 3. The Riemann problem is characterized by a weak pressure ratio and has the structure of the test presented in Sec. 5.1, with the contact discontinuity that moves at = 8.04 m∕s. Although we observe a higher speed in this test with respect to the test considering air, the Mach number is  even lower, below 0.01, due to the high speed of sound. Figure 5 displays the results obtained by solving (29)-(34) for one single phase. We have considered different time steps, which correspond also to acoustic CFL numbers greater than one (up to 12), as indicated in the caption by CFL(| | + ). The smaller is the time step, the better is the agreement with the analytical solution. Additionally, we test the capability to capture travelling material waves over a long simulation, by repeating the same Riemann problem test but over a longer time, that is = 0.095 s, as proposed in 100 . All test information are reported in the row lmWaterLong in Tab. 3. Figure 6 displays the detail of the solution field close to the contact discontinuity, which at the end of this test has reached = 0.7638 m, so it has crossed 7 grid cells. The quality of our results compares well with the ones reported by Abbate et al. 100 for their implicit scheme, so we can state that our scheme is able to correctly compute the position and the velocity of the moving material wave. refers to the standard formulation with density correction. In all cases, the momentum correction equation is used. refers to the standard formulation with density correction. Differently from Fig. 3, here the velocity correction equation (37) is used. Since the results do not differ notably from the ones obtained with the momentum correction, we report here only the zoomed regions, highlighted by rectangles in the first row of Fig. 3.

Lax problem
Finally, we end this single-phase section presenting the results for the Lax shock-tube test 101 , routinely used to validate standard compressible schemes, to investigate the behavior of the proposed scheme at Mach numbers between 0.3 and 1, so not so low. This shock test is characterized by an initial discontinuity also in the velocity, which is not null in the left chamber. Initial conditions and test data are given in Table 3, in the row Lax. When the diaphragm bursts, the initial discontinuity evolves in a leftward moving rarefaction waves and a rightward moving shock wave, with a contact discontinuity in between. The results at = 0.12 s are shown in Fig. 7 for different time steps Δ . In the right part of the domain, where the Mach number is higher, the numerical solution does not agree well with the analytical one. However, this discrepancy is an expected manifestation of the non-conservation of the total energy, but beyond that, the numerical results of the proposed scheme show an acceptable agreement with the analytical ones although we are not operating within the target regime of weakly compressible flows.

NUMERICAL RESULTS FOR THE HYPERBOLIC OPERATOR FOR TWO-PHASE FLOWS
In this section, we present two-phase flow results computed by using only the hyperbolic operator, without any relaxation process. Results of the complete numerical methods are shown in next section. Taking into consideration the conclusions of the previous section, we use the standard formulation with density re-computation, i.e. (29)- (34). We organize our analysis in subsequent steps, starting from the numerical validation of two fundamental properties: the behavior of the hyperbolic operator without mixing in Sec. 6.1, and the fulfillment of the pressure non-disturbance condition in Sec. 6.2. Then, we present the results of the proposed method on some reference Riemann problems available in the literature about BN-type models, in Sec. 6.4, and, finally, on a water-air mixture problem in Sec. 6.5.

No mixing water-air test
The first two-phase test involves liquid water and air governed by the stiffened gas model with the parameters listed in Tab. 2. As initial condition, the phases are uniformly dispersed with equal volume fraction 1 = 2 = 0.5, in a shock-tube where a mild pressure jump is imposed between the two chambers: = 100 bar at the left, and = 50 bar at the right. A null velocity and the temperature = 270 K are applied uniformly in the domain. The initial position of the discontinuity is at d = 0 and the grid spacing is Δ = 0.001 m.
Being the volume fraction uniform and given the absence of relaxation terms, in this test, each phase evolves independently from the other one. Thus, the exact solution can be computed by solving the Riemann problem for the Euler equations. Figure 8 shows the results at the final time of 0.16 ms computed with two different time steps: the smallest one corresponds to an acoustic CFL number slightly above 1 only for the liquid, while the biggest time step results in acoustic CFL number greater than 2 for 24 B. Re and R. Abgrall

Pure advection water-air problem
Here, we investigate a pure advection problem: a column of water-air mixture with a liquid volume fraction 1,c = 0.9 is transported at a velocity of 100 m∕s in a uniform pressure field at = 1 bar, involving a mixture with 1,L = 1,R = 0.1. The initial temperature is 270 K for both phases. The parameters of the stiffened gas model for the fluids are the same as in the previous test, and are listed in Tab. 2. Initially, the column is located at 0. This test is performed considering different discretizations, all imposing the convective CFL = 0.5. The results at time = 3 ms over three grids (with 400, 800, and 1600 cells) are shown and compared to the exact solution in Fig. 9. From the second row of the picture, we can appreciate that the no pressure or velocity oscillations arise and the initially uniform fields are correctly preserved during the time evolution. This achievement is crucial for a correct discretization of the non-conservative terms 69 . Beyond pressure and velocity, a good agreement between the numerical and the exact solution is observed also for the volume fraction and mixture density variables, for which the smearing of the contact discontinuity decreases with the grid refinement. To confirm this behavior, we have performed also a grid convergence study, presented in Fig. 10, computing the discrete 1 error between the numerical and the exact mixture density at the final time, normalized by the 1 norm of the initial mixture density, as The numerical error converges with the order of Δ 1∕2 , as expected when using a first-order scheme for BN-type models 89 .
B. Re and R. Abgrall

FIGURE 9
Pure advetion test of a column of water-air mixture at uniform velocity, results at = 3 ms, with CFL(| |) = 0.5. The initial position of the column is shown as a dashed line in the top row, which displays also the water volume fraction 1 on the left, and the mixture densitȳ = 1 + 2 on the right; the results at obtained over three different grids (with the number of primary cells) are compared with the exact solution shown by a solid black line. The second row presents the pressure and velocity fields at for both phases; only the results computed over the coarsest grid are shown for brevity, as no differences are observed using finer grids.

FIGURE 10
Pure advetion test: grid convergence study. The considered grid spacings correspond to 100⋅2 [1∶7] points and the 1 norm of the error on the mixture density is computed according to (67). The dashed line displays the convergence rate of Δ 1∕2 .

Verification using a manufactured solution
In this subsection, we perform a grid convergence study for a test with non-uniform pressure and velocity. To do that, we switch off all relaxation terms and add to the right-hand side of the equations a source ( , ) to be determined by a manufactured solution. We follow the strategy proposed by Hennessey et al. 102 , and we express the exact solution in terms of the primitive variables = 1 , 1 , 1 , 1 , 2 , 2 , 2 T as ex ( , ) = + 1 + 1, + 2, where , , 1 We run this test over different grids representing the domain 0 ≤ ≤ 1, and we compute the final solution at time = 0.5. The exact solution is given by evaluating Eq. (68) at the final time. The -infinity norm of the error on the mixture density, is shown in Fig. 11. The order of convergence is close to 1, as expected while using the Rusanov scheme, as here.

Reference Riemann problems with perfect gases
The goal of this section is to validate the proposed approach through some tests commonly used in the research community devoted to the development of one-dimensional numerical schemes for the BN-type models. Neglecting tests involving strong shock waves or vanishing phases, we have selected from the literature three Riemann problems for which the analytical solution is given: the first two, i.e., the sonic point and the 123-problem, are taken from 103,104 (named there Test 3 and Test 4, respectively); the third one reproduces the Test-case 1 in 105 , and it is called solid contact in the following.
Before describing each one, let us remark that these tests are not properly representative of low-Mach problems, but they provide anyhow an important contribution for the verification of the hyperbolic operator. Moreover, in these three tests, both fluids follow the perfect gas model, with the air parameters given in Tab. 2, so we prefer to use the notation phase 1 and phase 2, rather than solid and gas. Finally, for the sake of completeness, we report here the estimate for the shock speed 103 we have used to draw the analytical solutions: where the subscript pre and post refer to the pre-and post-shock states and the plus and minus sign is used for a right and a left traveling shock, respectively.

Sonic point
This test was presented in 103 to assess the correct resolution of a sonic rarefaction. The two phases have initially the same pressure, density and velocity, but the mixture composition differs between the left and right states:  103 . The convective CFL with respect to the shock speed, which is ≈ 2 ∕ , is CFL( ) = 0.24, which is pretty similar to the acoustic one since CFL(| | + ) = 0.32 for both phases. some discrepancies are expected as we are using a pressure-based, that is non-conservative, solver. Moreover, the asymmetry between phases in the density profiles is simply due to the smearing of the volume fraction discontinuity. More important, in this test, is the correct resolution of the sonic rarefaction, without any non-physical entropy glitch at the sonic point.

Two-phase 123-problem
This test involves a region close to vacuum, so it is useful to assess the pressure positivity. Initially, the fluids are at uniform pressure = 0.4 Pa and density = 1.0 kg∕m 3 . A discontinuity is imposed in the middle of the domain, = 0: on the left, the volume fraction is 1 = 0.8 and the velocity is −2.0 m∕s, on the right, the volume fraction is 1 = 0.5 and the velocity is +2.0 m∕s. The solution consists in two symmetric rarefactions and a stationary contact discontinuity in between, where, at the final time = 0.15 s, the pressure and density are extremely small: ( , ) = 0.0019 Pa and ( , ) = 0.0219 kg∕m 3 103 . The numerical results are displayed in Fig. 13. The pressure and the density are computed accurately, preserving the positivity. However, the discontinuity in the volume fraction appears to be very diffused, but a similar behavior for the Rusanov's scheme is reported also by Coquel et al. 105 .

Solid contact
The last Riemann problem we present was proposed by Coquel et. al. 105 and, differently from the previous ones, it does not involve an initial symmetry between the two phases. The initial field is described in Tab. 4, and its evolution encompasses seven different types of waves: for phase 1, a left-traveling shock, a material contact discontinuity moving at velocity 1 , a phase fraction discontinuity moving with velocity 2 and a right-traveling rarefaction wave; for phase 2, a left-traveling rarefaction fan, the phase fraction discontinuity, and a right-traveling shock. To be able to compare our results with the analytical and numerical solution in 105

Water-air mixture test
In this section, we reproduce the test proposed in 19 under the name Smooth shock tube test case. The fluids are water and air but, differently from Sec. 6.1, an initial discontinuity is imposed also in the volume fraction. The water is modeled under the stiffened gas model, using ∞ = 6.0 ⋅ 10 8 Pa as in 19 , while the remaining EOS parameters are the same given in Tab. 2. Initially, the fluids are at rest and the densities, w = 1050 kg∕m 3 for the water and a = 1.2 kg∕m 3 for the air, are uniform along the tube.  Fig. 15. In the first case, the maximum acoustic CFL is smaller than one (max ( + | |) ≈ 0.5 for both phases), whereas in the second case, the acoustic CFL is greater than 2 for both phases, with convective CFL reaching max a CFL(| |) = 0.75 for air and max w CFL(| |) = 0.0005 for water. An estimation based on the values of the solution variables across the volume fraction discontinuity leads to a value for the interface velocity ≈ 0.6 m∕s; since at the final time step it has moved only about Δ ∕10, its displacement cannot be distinguished in Fig. 15. As a reference, Fig. 15 displays also the results for the air reported in 19 . The match is reasonably good in proximity to the rarefaction wave and the contact discontinuity. On the contrary, the shock position is not captured correctly by the present method. This is an inherent limitation of the adopted pressure-based formulation and we are aware that the introduced error increases with the shock strength, however this test is on the boundary of the target application area, as the maximum Mach number of the air is well above one. Concerning the water results, we cannot compare them with the results in 19 , as their model is based on different assumptions which make the dispersed phase-water in this test-invariant across the shock. Nevertheless, as pointed out also in 19 , such a large pressure disequilibrium between water and air in this test is not much physically reasonable. Indeed, this test serves mainly as a further validation for the hyperbolic operator, especially when dealing with different scale velocities between the phases and acoustic CFLs greater than one: even with the largest time step, the scheme is stable and no spurious oscillations appear in the solution.

NUMERICAL RESULTS FOR TWO-PHASE FLOWS WITH RELAXATION
In this section, we finally present the results obtained with the full numerical scheme, that is including velocity and pressure relaxation. In Sec. 7.1, the results of the BN-type model with pressure and velocity relaxation are compared to the analytical results of Kapila's model for two-phase flows in mechanical equilibrium, and the role of the finite relaxation parameters is investigated. Then, the water-air mixture test of Sec. 6.5 is re-run in Sec. 7.2 with pressure relaxation to compare the results of the proposed model to the ones achieved through a different numerical method for a BN-type model. The last three subsections

FIGURE 15
Water-air mixture test: results at the time = 350 μs, over a uniform grid with 650 cells, in absence of relaxation terms. Two different time steps are considered: the smallest one (resulting from = 500 steps, and displayed in blue for water, and yellow for air) corresponds to max CFL(| | + ) ≈ 0.5 for both phases; the largest one (resulting from = 125 steps, and displayed in orange for water, and violet for air) corresponds to max CFL(| | + ) ≈ 2.2 for both phases. The initial condition exhibits a pressure and volume fraction discontinuity at = 0, but its displacement (approximately 0.21 mm) is too small to be observed in these graphics. The left column illustrates the water volume fraction and the pressures; the right column displays the partial densities , the velocity and the Mach number of water (scales on left axis) and air (scales on the right axis). The pressure and the velocity of the air, as well as the liquid volume fraction, given in Saurel et al. 19 are displayed as reference.

TABLE 5
Stiffened gas parameters for the pure fluids used in the two-phase numerical tests with relaxation involving a wateraluminum mixture (Sec. 7.1), a water-air mixture (Sec. 7.2), and almost-pure water and air flows (Sec. 7.4).

Water-aluminum mixture test
The test presented in this section was proposed by Furfaro et al. 106 and involves a mixture of two condensed phases, water and aluminum. The parameters of the stiffened gas models used in this test are summarized in Tab. 5. Initially, the phases are uniformly dispersed with equal volume fraction 1 = 2 = 0.5 in a shock-tube with two chambers: the left one at high pressure , and the right one at low pressure ( = 10 5 Pa). The densities are uniform: w = 1000 kg∕m 3 for water and al = 2700 kg∕m 3 for aluminum. The initial velocity is zero everywhere. Due to relaxation processes, as the time evolves, the volume fraction changes across the expansion and compression waves, besides the contact discontinuity originating from the initial pressure jump. To build a reference solution, we consider the mechanical equilibrium model of Kapila 12 , which is the limit model of BN model considering instantaneous pressure and mechanical relaxation, and the exact Riemann solver proposed in 107 . The computed reference solution is given in Tab. 6. The maximum Mach number is 0.06 (achieved by water) and this motivates the choice of this test to illustrate the capabilities of the proposed pressure-based method.
We start the illustration of the results considering a set of relaxation parameters sufficiently large to drive water and aluminum toward mechanical equilibrium, namely = 10 9 kg∕(m 3 s) and = 10 5 m s∕kg. The results at the time = 111 μs obtained with the proposed method correlate well with the reference solution, as Fig. 16 shows. A good match is reached also considering acoustic CFLs higher than one, indeed in the simulation carried out using = 200 time steps we have max CFL(| |+ ) w = 1.5 and max CFL(| | + ) al = 3.0. Considering these parameters, we conducted also a grid convergence study, illustrated in Fig. 17. As expected, a coarse grid leads to smoother results, but with grid refinement, the solution converges toward the reference one.
After the preliminary verification of the results, we use this test to illustrate now the effects of finite relaxation parameters. In the previous tests, we have seen how finite, large values can successfully replicate the mechanical equilibrium. Intuitively, smaller relaxation parameters leave a higher level of dis-equilibrium between phases. This behavior is confirmed by Fig. 18, which displays the results obtained with different sets of relaxation parameters: • the strong set with = 10 9 kg∕(m 3 s) and = 10 5 m s∕kg used also in the previous tests, • an intermediate set with = 10 8 kg∕(m 3 s) and = 10 4 m s∕kg, • a mild set with = 10 7 kg∕(m 3 s) and = 10 3 m s∕kg.
These tests are performed over a grid with Δ = 0.001 m at acoustic CFLs higher than one. From Fig 18, we note a substantial dis-equilibrium in the phasic velocity, in particular while using the mild set, but it can be observed also for the intermediate set.
To have a quantitative idea, the maximum velocity difference in these sets is Conversely, no disequilibrium in the pressure can be distinguished, but the effect of decreasing relaxation parameters is a smoothing of the expansion and compression waves, which causes the disappearance of the intermediate uniform regions between them (that is the regions called left* and right* in the mechanical equilibrium solution in Tab. 6). This behavior is reflected also by the volume fraction profile, where this smoothing diminishes also the degree of the mixing, as with the mild set the water volume fraction w reaches only a maximum of 0.519 after the rarefaction and a minimum of 0.475 after the compression, instead of max( w ) = 0.522 and min( w ) = 0.470 as in other cases. Finally, we better investigate the effects of pressure relaxation only, without relaxing the velocities (i.e., using = 0). Figure 19 compares the results obtained without pressure relaxation (so full disequilibrium, considering the hyperbolic operator only), with a weak relaxation parameter = 10 −1 m s∕kg and with a strong one = 10 5 m s∕kg. To allow for the faster waves in the first  The top-left panel shows how the initially uniform volume fractions change due to the relaxation processes. The top-right panel displays the phasic densities: please, note the different scales for water (left axis) and aluminum (right axis). In the bottom line, we have the pressure and the velocity of each phase: for each simulation, it is impossible to distinguish between water and aluminum as the relaxation processes drive them toward the equilibrium. In these plots, the reference pressure and velocity are only one, displayed as solid line. test, we consider a longer domain, Ω = [−0.8, 0.8] m, but the same grid spacing and time step as before. The first observation we can make from Fig. 19 is that the weak parameter is sufficient to drive the phasic pressures to the equilibrium, as the lines for water and aluminum appear to be overlapped also in the zoomed view. These lines overlap also to the ones corresponding to the strong relaxation and the Kapila's model, confirming that the pressure equilibrium is achieved in both tests. Similar observations between the two different relaxation parameters can be drawn also for the volume fraction and the velocity fields, which, conversely, differ significant from the Kapila's model because of the velocity dis-equilibrium. Furthermore, comparing the results with and without pressure relaxation, we can observe that wave speeds in the relaxed condition are similar to the wave speeds in non-equilibrium water, that is much slower than the aluminum ones. The equilibrium pressure is also pretty similar to the non-equilibrium water, whereas the equilibrium velocity is an intermediate value between the ones of water and aluminum.

Water-air mixture test with pressure relaxation
In this section, we reconsider the Smooth shock tube test case proposed in 19 , but differently from Sec. 6.5, we use now pressure relaxation, and we compare our results again with the ones shown in 19 under stiff pressure relaxation. The initial conditions and the stiffened gas parameters of air and water are those given in Sec. 6.5 and in Tab Figure 20 shows the results computed at = 350 μs, with two different time steps-the larger one resulting in an acoustic CFL of about 1.85 for both phases-and a relaxation parameter = 10 5 m s∕kg. The agreement with the reference results is not excellent, probably because the maximum Mach number for air is higher than 1.2, so our non-conservative scheme is not able to correctly capture the shock speed and this inaccuracy affects also the pressure profiles, where we have a difference of about 0.3 bar (i.e., 4%) after the rarefaction and of 0.2 bar (i.e., 3%) after the shock. We also remind that the model used in 19 is not symmetric and it is based on different modeling hypotheses for the water and we are not able to estimate whether and how this variation in the models may justify the discrepancy in the results. Nevertheless, the general behavior of the pressure and velocity of both phases is consistent with the reference results, and this serves as a confirmation of the correctness of the pressure relaxation scheme.
To better investigate the effect of the finite pressure relaxation parameter, we have re-run this test considering the largest time step and different values of . Figure 21 compares the pressure profiles and the differences in the phasic pressure obtained with two values of . The value = 10 5 m s∕kg is large enough to drive the phasic pressures toward the equilibrium, and the phase difference is null everywhere. On the other hand, the small value = 10 −1 m s∕kg allows a certain degree of phase The solution of Tab. 6 is also displayed for comparison with previous figures, but here it is labeled as "Kapila's model" as it is not a reference solution for this non-equilibrium test. Markers distinguish lines which, especially for pressure, overlap between water and aluminum (filled vs. empty symbols) and between the two tests with pressure relaxation (differentiated by yellow squares vs blue downward-pointing triangles). A zoomed view of the pressure field close to the rarefaction end is shown in the top-left corner. disequilibrium, especially across the shock, where it reaches the maximum (206.7 Pa), and across the material discontinuity, where there is also a sign change. Indeed, due to the smearing of material discontinuity, the water on its left undergoes an expansion (the increase of the volume fraction has an effect similar to an expanding nozzle) and the water on its right is slightly compressed, and the weak pressure relaxation is not sufficient to overcome this phenomenon.

Two-phase water expansion tube
The benchmark presented in this section was originally proposed in 22 and it has been investigated also in 25,87 , with and without phase transition. Here, we consider it without phase transition. This test involves a tube filled with liquid water at pressure = 1 bar and density liq = 1150 kg∕m 3 , to which a weak volume fraction of vapor vap = 0.01 is added. As in 87 , we compute the vapor initial conditions from the pressure and temperature of the liquid, considering the stiffened gas parameters given in Tab. 7. An initial velocity discontinuity is located at the center of the tube ( = 0): the left velocity is −2 m∕s, the right one is 2 m∕s.
The solution is computed at time = 3.2 ms and it consists of two symmetric rarefaction waves moving outwards. At the center of the domain, the volume fraction increases because of the vapor mechanical expansion and a gas pocket is dynamically generated 22 . As for the water-aluminum test in Sec. 7.1, a reference solution is computed in the limit of stiff mechanic relaxation according to the Riemann solver proposed by Petitpas et al 107 . The numerical results computed using a grid spacing Δ = 10 −4 m are shown in Fig. 22. First, we compute the solution using = 160000 time steps, corresponding to a maximal acoustic CFL of 0.3 (for the liquid). With this choice, we achieve a good match with the reference solution, especially in terms of pressure, capturing correctly the constant region between the rarefaction waves, while a little overshoot affects both the density and the volume fraction, but this behavior has been exhibited also in the previous works. Considering the small velocities and the severe rarefactions involved in this test, the capability of the proposed low-Mach scheme to correctly compute the solution at mild CFL numbers is a notable result. Remarkably, we perform a second test enforcing a 100 times bigger time step, leading to maximal acoustic CFL number of 28.6 for the liquid and 15.8 for the vapor. As we can see from In the left column, we see the water volume fraction and the pressures, which overlap between phases because of relaxation; in the right column, we have the partial densities and the velocity of each phase (note the different scales: for water on left axis, for air on the right axis). Thin, black lines display the pressure, the velocity and the volume fraction in the limit of stiff pressure relaxation given by Saurel et al. 19 .

FIGURE 21
Water-air mixture test with different pressure relaxation parameters: results at the time = 350 μs, over a uniform grid with 600 cells, in absence of velocity relaxation ( = 0), considering = 150 time steps. The compared values are 10 5 (blue) and 10 −1 m s∕kg (orange). On the left, the water and air pressure profiles (in bar) are displayed, but it is impossible to distinguish any difference. On the right, the pressure disequilibrium between water and air (in Pa) is plotted over the domain.
behavior of the fluids is well captured. By comparison, Zein et al. 25 used a CFL equal to 0.03 to obtain a stable solution, while Han et al. 87 used a maximum CFL of 0.9.

Almost pure fluids test
Here, we present a test involving almost pure fluids, water and air, with material parameters given in Tab. 5. The configuration involves a shock-tube with a left chamber ( < 0) filled with air at = 100 bar and the right chamber ( > 0) filled with water at = 50 bar. Since the proposed method is not able to deal with null volume fractions, we define the air volume fraction in the left chamber as a = 1 − and the water one in the right chamber as w = 1 − , where = 10 −4 . The air and water densities

FIGURE 22
Water expansion tube test: results at the time = 3.2 ms, over a uniform grid with 20000 cells, with pressure and velocity relaxation ( = 10 7 kg∕(m 3 s) and = 10 5 m s∕kg). The results obtained with = 160000 and with = 1600 time steps are displayed, along with the reference solution, computed assuming mechanical equilibrium between phases according to 107 . The displayed variables are: on the top, vapor volume fraction (initially, vap = 0.01 everywhere) and mixture densitȳ = liq + vap (it follows the color of Vapor in the legend); on the bottom, pressure and velocity of each phases (Water in the legend refers to liquid water). are a = 100 kg∕m 3 and w = 1000 kg∕m 3 , respectively. Starting from a state of rest, the almost pure liquid on the right is set into motion by the almost pure gas at higher pressure on the left. Considering pure fluids separated by a material discontinuity, it is possible to compute the analytical solution according to the Euler equations. It comprises a left-traveling rarefaction for air and a right-traveling shock wave for water, while the intermediate state is characterized by a pressure ⋆ = 98.887 bar and a velocity ⋆ = 2.989 m∕s. This velocity confirms the low-Mach regime, as the maximum Mach number is 0.008, and the shock speed is 1636 m∕s (1.0025 times the pre-shock speed of sound).
The numerical results at = 0.8 ms, obtained with pressure and velocity relaxation, are displayed in Fig. 23 and they reach an excellent agreement with the analytical, pure-fluid solution, except the shock wave which is considerably smeared. Figure 23 shows also how the solution varies for higher values of , that is for higher levels of mixing. The figure contains some zoomed frames to highlight the details close to the most important flow structures. The level of mixture given by = 10 −2 , that is the same as the one considered for the liquid-water expansion test in sec. 7.3, is already enough to significantly depart from the pure-fluid behavior. Indeed, the intermediate velocity is higher, resulting in a faster material discontinuity and a slower shock. On the other hand, = 10 −3 is already sufficient to capture qualitatively the behavior of the pure fluids, but it leads to some discrepancies in the values of the intermediate state. These are overcome by using the value = 10 −4 .
To investigate better the capability to correctly capture the interface between almost pure fluids, we repeat this simulation, for a long time. We consider a longer domain, Ω = [−20, 60] m, with the same grid spacing (Δ = 10 −3 m) and we compute the solution at = 0.03 s, using the time steps Δ = 0.2 and Δ = 1.0 , corresponding to a maximal acoustic CFL for water of 0.33 and 1.64, respectively. The results are displayed in Fig. 24. The final position of the material interface is ⋆ = 0.0897 m = 1200 time steps, and imposing pressure and velocity relaxation. On the top-left, the water volume fraction is displayed, along with a detailed view close to the material discontinuity, where the dashed vertical line illustrates the initial position of the material discontinuity. On the top-right, the density, or more precisely the mixture densitȳ , is displayed, with a zoom of the region close to the shock; the disagreement among the initial densities (clearly visible for the blue line at = 1.4) is due to the different values of , which lead to different weights in the average of the densities. On the bottom line, we have the pressure profiles, with a zoomed view of the rarefaction, and the velocities. In all frames, the solid black line displays the pure-fluid solution computed analytically. and it is computed correctly, as it can be observed in the zoomed view in the top-right corner of Fig. 24. Similarly, also the intermediate state and the rarefaction fan are in excellent agreement with the analytical solution. On the contrary, the position of the shock is not computed accurately, because the non-conservativeness of the scheme introduces an error in the shock velocity, which, after a long time, results visible in the position. However, the relative error in the shock position is about 1.0%, so it is acceptable.

Two-phase carbon dioxide test in saturation conditions, with accurate EOS
In this final section, we compare the results computed under the stiffened gas approximation with the ones computed using a more accurate thermodynamic model for CO 2 , based on the Peng-Robinson (PG) EOS 84 . In particular, we take advantage of the implementation of this EOS tailor-made for CO 2 flows, provided by an in-house thermodynamic library developed at SINTEF Energy 108 , that exploits the concept of the corresponding states to enhance the accuracy of specific properties, such as the density and the speed of sound, for the liquid phase.
The set-up in which this comparison is carried out consists in a pipe, 600 m long, filled with two-phase CO 2 mixture at saturation conditions, that is the initial conditions of the fluids, represented in the thermodynamic plane, lie along the saturation, or vapor-liquid equilibrium, curve. 3 In the left part of the domain, the CO 2 mixture is composed by 75% of liquid ( liq = 0.75) and 25% of vapor ( vap = 0.25), at the saturation temperature = 260 K. In the right part, the composition is inverted and the saturation temperature is 20 K higher, that is liq = 0.25, vap = 0.75 and = 280 K. We remind that, at saturation conditions, defining the temperature unequivocally defines also the pressure and the density of each phase. Hence, in the test with the Peng-Robinson EOS, no other inputs are defined. Conversely, in the test with the stiffened gas model, as we do not have a saturation We compute the flow field for = 1 s from the moment the diaphragm separating the two mixtures, initially at rest, is removed. We consider = 2000 time steps and a grid spacing Δ = 0.1 m, which result in max (| | + ) liq = 2.3 and max (| | + ) vap = 1.2. Higher values of CFL could be enforced without preventing stability, but we found that these values leads to a good compromise between accuracy and efficiency. Figure 25 illustrates the profiles of the liquid volume fraction, the pressure and the velocity, which are driven toward the equilibrium by relaxation processes with = 10 7 kg∕(m 3 s) and = 10 3 m s∕kg. The different thermodynamic models lead to a difference in the minimum velocity reached in the intermediate region between the shock and the rarefaction wave. This is reflected on the position of the material discontinuity: as highlighted by the zoomed view on the top-right box in Fig. 25, the stiffened gas predicts a larger displacement (in the negative -direction) of the contact discontinuity due to faster fluids velocities. Quantitatively, the position of the contact discontinuity is −11 m using the Peng-Robinson EOS, while it is −11.3 m using the stiffened gas.
The differences introduced by the thermodynamic models are clearer in the density profiles, displayed in Fig. 26. Given as input the same temperatures and pressures, the initial conditions entail already discrepancies in the density: ,vap < Stif f ,vap . Moreover, Peng-Robinson EOS predicts stronger initial density jumps for both phases. This leads to a notable disagreement in the behavior of the CO 2 vapor after the diaphragm rupture, when the two thermodynamic models predict opposite (in sign) density jumps across the material discontinuity.
Finally, we plot the results obtained with the Peng-Robinson EOS in the pressure-volume thermodynamic plane, in Fig. 27. This illustration gives a clear idea of how the proposed full non-equilibrium BN-type model allows each phase to evolve independently according to its own thermodynamic model, although the phasic pressures (and velocity) are immediately driven toward

FIGURE 25
Two-phase CO 2 test: results at the time = 1.0 s, over a uniform grid with 6000 cells (Δ = 0.1 m), considering = 2000 time steps. Initially, the left ( < 0) and right ( > 0) parts of the tube contain a saturated mixture with different composition and at different temperatures. The results obtained with the Peng-Robinson EOS (blue and triangular markers, PG EOS in the legend) are compared to the ones obtained with the stiffened gas model (orange and square markers, Stiff Gas in the legend). On the top, the liquid volume fraction is displayed, with a detailed view close to the material discontinuity in the box on the right. On the bottom, the pressure and the velocities are displayed: the liquid (full markers) and vapor (empty markers) CO 2 are driven to the equilibrium by means of relaxation. the equilibrium. The resulting flow states lie in close proximity to one side (liquid on the left, vapor on the right) of the dome, and only the evolution of the mixture density develops in the fully two-phase region. However, no special treatment, such as the definition of a speed of sound for the mixture, is required for this situation.

SUMMARY AND CONCLUSIONS
Starting from the symmetric variant proposed by Saurel and Abgrall 11 , we derived a pressure-based BN-type model for weakly compressible two-phase flows, which, as illustrated by Remark 1, encompasses the divergence-free condition for multiphase flows in the zero-Mach limit, recovering the correct scaling of the pressure. Although the single components of the finite-volume solver we built to solve the 1D pressure-based BN-type model are not novel by themselves, the resulting numerical method is unique and includes notable features. For instance, it preserves by construction the non-disturbance condition on pressure and velocity, it overcomes the stringent limitation on the time step imposed by the acoustics, it allows for different thermodynamic models, and it provides a wide variety of modeling choices thanks to the relaxation terms with finite parameters, which can be used to control whether and how the phasic pressure and/or velocity are driven toward the equilibrium.
The main drawback of the proposed method relates to the non-conservative character of the pressure-based formulation, which prevents to resolve accurately shock waves. However, this behavior was expected and numerical tests showed that, for moderate intensities, the error in the position of the shock is between 1 and 2%, consistently with the predictions made by Karni 36 . Notwithstanding the acceptable impact on weakly compressible flows, a remedy for this limitation is imperative to FIGURE 26 Two-phase CO 2 test: comparison of the density ( liq and vap ) at the initial time (black dashed and dotted lines) and at = 1.0 s (colored lines with markers) computed with Peng-Robinson EOS (black dashed; blue and yellows lines with triangular markers) and with stiffened gas model (black dotted; orange and violet lines with square markers). The upper part of the plot refers to the liquid CO 2 , while the bottom part refers to the vapor CO 2 . Note that the -axis is discontinuous, but the tick spacing is preserved.

FIGURE 27
Two-phase CO 2 test: visualization of the flow field computed using the Peng-Robinson EOS in the thermodynamic plane pressure-specific volume. The saturation curve is plotted on the background with a grey line. The initial states for the liquid (blue) and vapor (orange) phase, as well as the CO 2 mixture (yellow), are displayed by full triangles: left-pointing triangles refer to the left state (saturated temperature = 260 K) and right-pointing ones refer to the right state (saturated temperature = 280 K). The results at = 1.0 s are displayed by lines. The evolution of the mixture (dashed line) is plotted as 1∕̄ versus the mixture pressurē = liq + vap , where however the phasic pressures liq and vap are equal thanks to relaxation. move toward an all-speed scheme, so this matter will be given a high priority in the future development of the method. To mitigate the effects of non-conservativeness, we will consider strategies already proposed for single-phase flows-for instance, an additional scheme-dependent viscous term that provides a correction at the leading order of the discrete approximation 111 , or the use of the pressure equation only as a prediction for variables in the total energy equation 54 -and for reduced two-phase flows, e.g., considering a supplementary, explicit correction step 60 or moving to residual distribution schemes 61 .
As told in the introduction, this work describes the first stage of a longer project, which aims to develop a reliable and robust tool for multi-phase simulations of weakly-compressible flows. To reach this aim, the numerical scheme underlying the solution strategy described in this work should be enhanced: higher-order discretization techniques will be considered, in addition to the previously mentioned correction for non-conservativeness. Afterwards, we will extend the proposed numerical tool to multidimensional grids. In particular, to account for complex geometries, we will look at unstructured grids, for which a face-based staggered discretization has been recently proposed by Bermúdez et al. 52 in the framework of weakly compressible single-phase flows.
In conclusion, we propose here the first numerical tool based on the BN-type model that includes simultaneously all the features outlined in the first paragraph of this conclusive section. As identified above, the proposed method is not free from defects, but the main limitations could be overcome by adapting the solution strategy to include some of the techniques already proposed in the literature for other models. Given these points, we think this work could pave the way for a wider application of BN-type models, especially to investigate two-phase flows involving liquid components, which generally exhibit low Mach number and require specific, accurate thermodynamic models. The ingredients of the solution strategy, taken individually, are not complex, so the main ideas presented in this work could be applicable also in software used in applied research and industry, where simplicity and efficiency are of utmost importance. Similarly, the modular nature of the presented method could facilitate its integration in other BN solvers, mainly used in academic research, so far. Furthermore, the use of finite relaxation parameters may enhance the modeling capabilities of BN-type models. Indeed, if experimental data were available, parametric analyses could be carried out to define the best values (or range of values) for each multiphase regimes, contributing to increase the fidelity of numerical results.

APPENDIX A DERIVATION OF THE PRESSURE FORMULATION FOR THE BN-TYPE MODEL
In this appendix, we show step by step how we have achieved the pressure formulation of Eq. (10) Finally, recalling the definitions of the speed of sound (7) and (8)