Pressure-based algorithm for compressible interfacial ﬂows with acoustically-conservative interface discretisation

A pressure-based algorithm for the simulation of compressible interfacial ﬂows is pre- sented. The algorithm is based on a fully-coupled ﬁnite-volume framework for unstructured meshes with collocated variable arrangement, in which the governing conservation laws are discretised in conservative form and solved in a single linear system of equations for velocity, pressure and speciﬁc total enthalpy, with the density evaluated by an equation of state. The bulk phases are distinguished using the Volume-of-Fluid (VOF) method and the motion of the ﬂuid interface is captured by a state-of-the-art compressive VOF method. A new interface discretisation method is proposed, derived from an analogy with a contact discontinuity, that performs local changes to the discrete values of density and total enthalpy based on the assumption of thermodynamic equilibrium, and does not require a Riemann solver. This interface discretisation method yields a consistent deﬁnition of the ﬂuid properties in the interface region, including a unique deﬁnition of the speed of sound and the Rankine–Hugoniot relations, and conserves the acoustic features of the ﬂow, i.e. compression and expansion waves. A variety of representative test cases of gas–gas and gas–liquid ﬂows, ranging from acoustic waves and shock tubes to shock-interface interac- tions in one-, two- and three-dimensional domains, is used to demonstrate the capabilities and versatility of the presented algorithm in all Mach number regimes. The propagation, reﬂection and transmission of acoustic waves, shock waves and rarefaction fans in interfacial ﬂows are predicted accurately, even for diﬃcult cases that feature ﬂuids with shock impedance matching, transonic shock tubes or strong shocks in gas–liquid ﬂows, as well as on unstructured meshes. © Author(s). the BY license (http://creativecommons.org/licenses/by/4.0/). knowledge of the wave structure to evaluate ﬂuxes. Results for a variety of representative compressible gas–gas and gas–liquid interfacial ﬂows have been presented, including the propagation of acoustic waves, shock tubes and the interaction of shock waves with one-, two- and three-dimensional interfaces. The transmission and reﬂection of acoustic waves at the interface has been shown to be captured accurately by the presented algorithm, including cases with acoustic impedance matching, a capability not previously reported in the literature. Similarly, the interaction of shock waves with the interface in one-, two- and three-dimensional


Introduction
The numerical modelling of compressible flows is associated with a number of difficulties owing to the formulation of the conservation equations, the coupling of hydrodynamic and thermodynamic variables, and high pressure ratios. In addition to the numerical difficulties encountered when modelling compressible single-phase flows, the numerical solution of compressible interfacial flows, in which two (or more) immiscible fluids interact, is further complicated by different fluid properties and speeds of sound (and, thus, Mach numbers) of the bulk phases, complex acoustic behaviour as well as the numerical treatment of the fluid interface. It has proven particularly difficult to define the discrete interface between two compressible fluids in a consistent manner, which retains the main features of the solution, such as the propagation of acoustic waves and shock waves. According to Coralic and Colonius [1], numerical methods able to accurately predict the interaction of interfaces with shock waves, and in extension compressible interfacial flows in general, should satisfy the following criteria: 1. discrete conservation of mass, momentum and energy, 2. avoid the generation of spurious oscillations at the interface or at shock waves, and 3. provide high-order accuracy in smooth regions.
A variety of contemporary algorithms for compressible interfacial flows, e.g. [1][2][3][4][5][6][7], have been shown to discretely satisfy the global conservation of mass, momentum and energy in the computational domain, which is a prerequisite for the accurate prediction of the speed of shock waves [8][9][10], although the conservation in each individual phase is strongly dependent on the applied method. Flux-limited high-order schemes, most notably total variation diminishing (TVD) schemes [11,12], or high-order schemes combined with artificial viscosity models [13], can provide second-order accuracy in smooth regions and avoid oscillations at shocks and discontinuities. Spurious oscillations at the interface as a result of the discontinuity in fluid properties, however, remain a common issue for interface capturing methods [1], in particular when the bulk phases are assumed to be in mechanical equilibrium. Karni [14] and Abgrall [2] were among the first to depart from a fullyconservative discretisation, i.e. a discretisation that is simultaneously conservative with regards to the entire computational domain and each individual phase [15], in favour of avoiding spurious oscillations. The application of a non-conservative discretisation at the interface has also been suggested in the context of incompressible interfacial flows by Brackbill et al. [16], to avoid spurious oscillations as a result of the changing fluid properties at the interface. A non-conservative discretisation of the governing equations is widely, and largely successfully, applied for incompressible interfacial flows, although a naïve non-conservative discretisation does not allow compression and expansion waves to pass the interface. To this end, exact and approximate Riemann solvers are widely applied to allow the exchange of information between the interacting fluids through the interface. Recent studies [17,18] also suggest that a strictly conservative discretisation is important for the robustness of the solution algorithm and the predictive accuracy for incompressible interfacial flows with large density ratios.
State-of-the-art algorithms for compressible interfacial flows are almost exclusively founded on density-based algorithms, which solve the governing conservation equations (momentum, continuity and energy) for the momentum, density and total energy of the flow. In the two-fluid Baer-Nunziato (BN) model [19], each phase is treated as a separate fluid with its own momentum, continuity and energy equations, with an additional transport equation for the volume fraction field of the interacting fluids (e.g. a colour or level-set function), or a topological equation that represents the fluid interface. BN models, also often referred to as seven-equation models, have well-defined hyperbolic properties [20][21][22] and conserve the mass, momentum and total energy of each phase. However, BN models involve a considerable computational complexity due to the seven, nine or eleven governing equations for one-, two-or three-dimensional flows, respectively, which have to be coupled through relaxation terms for pressure and velocity (usually under the assumption of local thermodynamic equilibrium), including non-conservative products [21][22][23][24][25] and a consistently defined interface velocity [6,21], which is typically based on the solution of an exact or approximate Riemann problem in a Godunov-type, Rusanov-type or Roe-type discretisation.
By assuming equilibrium of pressure and velocity in the limit of infinitely fast relaxation to mechanical equilibrium [6,26,27], BN models can be reduced to five-equation models [1,4,[27][28][29][30][31], with a separate continuity equation for each phase, and shared momentum and energy equations. This model conserves the mass of each phase and globally conserves momentum and energy [28]. However, the total energy of each individual phase is not conserved and the relative volume change of the phases due to compression/expansion has to be incorporated in the interface transport equation [6,27]. Abandoning the conservation of a separate density field for each phase, and considering instead only a single density field, leads to a fourequation model [2][3][4] with a single conservation equation for momentum (one for each space dimension), for mass and for energy, plus an interface transport equation. This approach is also frequently referred to as the one-fluid model, because the entire two-phase flow (both phases and mixtures thereof) is numerically treated as one fluid with locally varying properties. It is the simplest model able to simulate interfacial flows, but does not by design conserve discretely the momentum, mass or energy of each phase. Furthermore, the four-equation model is associated with difficulties when recovering the pressure field in density-based frameworks [4,15,27], for which a specific equation of state (EOS), e.g. a stiffened gas EOS, is required. To this end, the five-equation model of Allaire et al. [4], for instance, reduces to the four-equation model of Shyue [3] when both phases are perfect gases.
In the vast majority of methods, such as in the seven-, five-and four-equation models discussed above, a Riemann solver is applied to evaluate the fluxes at cell-faces of the computational mesh, effectively reducing the evaluation of the fluxes to a one-dimensional problem for which a Riemann problem is solved [32]. Assuming an interface coincides with a cell-face, a Riemann problem with different EOS can be solved at the fluid interface to couple the bulk phases in an effective manner [33]. Riemann solvers have a strong mathematical foundation, but exact Riemann solvers are for many practical cases prohibitively time consuming, which motivates the development of approximate Riemann solvers. The Harten-Lax-van Leer Contact (HLLC) formulation, introduced by Toro et al. [34], has become the most widely used approximate Riemann solver for interfacial flows, as it captures contact discontinuities. However, because approximate Riemann solvers rely on an accurate a priori approximation of the wave speeds, their derivation for complex interfacial flows can be cumbersome and the resulting methods are often limited to a specific application or parameter range. Various HLLC-type solvers have been proposed for interfacial flows [1,[35][36][37], including phase transition [38,39] and surface tension [30,38,39]. The Ghost-Fluid Method (GFM), originally proposed by Fedkiw et al. [40], presents an alternative that does not, in general, require the solution of a Riemann problem. The GFM defines ghost cells on either side of the interface for the discretisation of the governing equations, where the real fluid and the ghost fluid coexist, and ghost values are defined in these ghost cells for the fluid on the opposing side of the interface. Due to its conceptual simplicity, the GFM has been used in a variety of numerical frameworks for compressible interfacial flows, notably [41,42]. The GFM has also been extended to interfacial flows with reactions [43,44] and phase transition [44], and has been used in conjunction with discontinuous Galerkin methods [45]. However, the original GFM was found to lack stability for strong shock-interface interactions and compressible gas-liquid flows, a shortcoming which was addressed by coupling the GFM to a Riemann solver to compute the flow states at the interface [46][47][48][49].
Density-based algorithms are arguably the method of choice for flows with considerable compressibility, but are illsuited for low-Mach number flows [50][51][52], requiring sophisticated preconditioning and solution methods, see e.g. [53][54][55][56][57]. The pressure field has to be reconstructed based on the applied EOS, since pressure is not a solution variable in densitybased algorithms. This is generally less of a problem when GFM [40,42,47] or similar so-called sharp interface methods are used, where either phase is present in a mesh cell but never a mixture of them. However, interfacial mesh cells (where the normalised interface indicator function, e.g. a volume fraction, is 0 < < 1) contain a "mixture" of both phases when interface capturing methods or diffuse interface methods are applied, such as the Volume-of-Fluid method (VOF). This fluid mixture, which Abgrall and Saurel [23] aptly called numerical mixture, is of numerical origin and, from a continuum mechanics standpoint, has no physical basis. Reconstructing the pressure in such interfacial cells is problematic [4,15,27], since hydrodynamically and thermodynamically plausible fluid properties as well as a mixture EOS have to be defined, which may misrepresent the interface region by presuming a physical meaning of the finite interface thickness. As a consequence, the speed of sound in the interface region can, for instance, be lower than the speed of sound of either bulk phase [6,[58][59][60], which is inconsistent with the zero thickness of the interface assumed in continuum mechanics.
Pressure-based algorithms, in which the continuity equation serves as an equation for pressure, while density is evaluated explicitly using a suitable EOS, are preferably applied for incompressible flows and yield significant advantages for low-Mach number flows, since density is not required to be a thermodynamic variable and the acoustic degeneration (i.e. vanishing pressure-density coupling) at low Mach numbers does not pose a problem. The success of pressure-based algorithms is facilitated by the unique role of pressure in all Mach number regimes, with the pressure-velocity coupling dominant at low Mach numbers and the pressure-density coupling dominant at high Mach numbers [61,62], as well as the convenient fact that the fully-conservative formulation of the governing conservation laws can still be satisfied accurately, even if non-conserved quantities, such as pressure, are chosen as primary solution variables [61]. However, ensuring a stable numerical solution in the transonic regime, where pressure is strongly coupled with both velocity and density, and formulating consistent shockcapturing schemes for pressure-based algorithms has proven to be difficult [51,63]. Although, since its original conception by Harlow and Amsden [64,65], many pressure-based algorithms have been proposed for compressible single-phase flows [50,61,[66][67][68][69][70][71][72][73][74][75][76], no pressure-based algorithm for compressible interfacial flows has been published yet.
In this article, a fully-coupled pressure-based algorithm for compressible interfacial flows is proposed, based on the finite-volume framework of Xiao et al. [76] for compressible single-phase flows on unstructured meshes. The presented algorithm uses a compressive VOF method [77] for the advection of the interface, leading to a four-equation model. The discretised governing flow equations (momentum, continuity and energy) are solved in a single system of linearised equations for velocity, pressure and specific total enthalpy, with density being evaluated based on the applied EOS. The present study considers the perfect-gas and stiffened-gas models to describe interfacial flows of gases and liquids. A new discretisation method at the fluid interface is proposed, derived from an analogy with a contact discontinuity, that performs local changes to the discrete values of density and total enthalpy in the interface region, assuming local mechanical and thermal equilibrium, whereby the conservative discretisation of the governing equations remains unchanged. Since the proposed algorithm is pressure-based and applies a single pressure field for the entire computational domain, no partial pressures have to be considered and no mixture pressure has to be defined at the interface. Instead, only mixture rules for the fluid properties have to be defined. As demonstrated by representative test cases of compressible gas-gas and gas-liquid flows, the presented algorithm predicts the propagation, reflection and transmission of acoustic waves, shock waves and rarefaction fans in interfacial flows accurately. In particular, the precise simulation of the reflection and transmission of acoustic waves at fluid interfaces as conducted in this study has not been reported in the literature before.
In Section 2, the governing equations are introduced. The numerical framework and the pressure-based algorithm are presented in Section 3 and the applied compressive VOF method is described in Section 4. The new interface discretisation method is proposed in Section 5, followed by a discussion of the iterative solution procedure in Section 6. The results for a variety of representative compressible interfacial flows are presented and discussed in Section 7. The article is summarised and concluded in Section 8.

Governing equations
The considered compressible interfacial flows of inviscid fluids at all speeds are governed, assuming Cartesian coordinates, by the momentum equations the continuity equation and the energy equation where t is time, ρ is the density, u is the velocity vector, p is the pressure and h = c p T + u 2 /2 is the specific total enthalpy, with c p the specific isobaric heat capacity and T the temperature. The enthalpy formulation is chosen for the energy equation, rather than the more common internal energy formulation, as it leads to a straightforward application in the proposed numerical algorithm, since the transient pressure term on the right-hand side of Eq. (3) does not require linearisation. The governing equations require closure by defining the thermodynamic properties of the fluids through an appropriate EOS. In this study, the stiffened-gas model [78] is applied, which provides a very good description of liquids and solids of practical interest [33,79], and reduces to the perfect-gas model for gases. It is, therefore, widely used for interfacial flow modelling. The thermodynamic properties of the fluid are linked via the stiffened gas EOS where R is the specific gas constant, γ = c p /c v is the heat capacity ratio, c v is the specific isochoric heat capacity and is a material-dependent pressure constant, which is = 0 for an ideal gas. The speed of sound for a stiffened gas is (5) from which the specific total enthalpy follows as with the specific isobaric heat capacity c p = c p,0 p + p + γ (7) and c p,0 = R/(1 − γ −1 ). For = 0, the stiffened-gas model reduces to the perfect-gas model.
In order to distinguish the interacting bulk phases, the VOF method [80] is applied to represent the bulk phases by an indicator function ψ(x), typically called colour function, with where = a ∪ b is the computational domain, with a and b the subdomains occupied by fluid a and b, respectively. Consequently, the interface is located in cells with 0 < ψ < 1. Since the density is discontinuous at the fluid interface and fluid does not flow through the interface (assuming no mass transfer), the interface between two immiscible fluids represents a contact discontinuity [23,78]. Thus, the fluid interface is a material front propagating with the flow and, consequently, the material derivative of ψ is Accounting, in addition, for the different acoustic properties of the bulk phases [6,59], Eq. (9) becomes ∂ψ ∂t is a material-dependent compressibility factor [6,59,60], with ρ given by Eq. (4) and a given by Eq. (5). Assuming the bulk phases are in mechanical equilibrium, and since no mass transfer and surface tension are considered in this study, the interface conditions for velocity and pressure are u a · m = u b · m and p a = p b [81], respectively, where m is the normal vector of the interface.

Numerical framework
The presented algorithm is founded on a fully-coupled pressure-based numerical framework, based on the algorithm for compressible single-phase flows of Xiao et al. [76]. The framework is predicated on the finite-volume method with a collocated variable arrangement, and all discretisation methods presented below are applicable to unstructured meshes. Since the discretisation of the governing equations is largely identical for single-phase and interfacial flows, the discretisation of the governing equations presented in this section focuses on single-phase flows. The modifications to this discretisation for interfacial flows, including the appropriate definition of the fluid properties, is described in Section 5 and the iterative solution procedure applied to solve the discretised, linearised governing equations is described in Section 6.

Finite-volume method
The finite-volume method, which forms the foundation for the numerical framework, is based on the integral formulation of the governing conservation laws. It is worth recalling that only the integral form of the governing conservation laws is valid at shocks and discontinuities [32]. Integrating, for example, the continuity equation (2) over the control volume V , the integral form of Eq. (2) is given as follows from the divergence theorem, with S the outward pointing surface vector on the surface dV of control volume V . Assuming the surface of the control volume is constituted by a finite number of flat faces f and applying the midpoint rule [62,82], the surface integral can be expressed with second-order accuracy [82] as where A f is the area of face f and ϑ f = u f n f is the advecting velocity, with n f the outward pointing unit normal vector of face f . Because a collocated variable arrangement is employed, all primary solution variables (u, p, h) as well as the density ρ are stored at cell centres. Therefore, the data required at face centres is obtained by interpolation from the data at adjacent cell centres, as further explained in Section 3.2, while the face velocity u f requires a special interpolation to prevent pressure-velocity decoupling, discussed in detail in Section 3.4. The discretisation of the temporal term is discussed in Section 3.3. The interested reader is referred to the work of Xiao et al. [76] for further details on the applied finite-volume method.

Spatial discretisation
The central differencing scheme is applied for the interpolation from cell centres to face centres of variables that are not advected, given for a general flow variable φ as The geometry coefficient l f is defined based on an inverse distance weighting, given as l f = |r P f |/ s f , where s f is the distance between cell centres P and Q (the cell centres adjacent to face f ), schematically illustrated in Fig. 1a, and r P f is the vector connecting cell centre P with face centre f . Advected variables are interpolated to face centres using the TVD interpolation method for general unstructured meshes proposed by Denner and van Wachem [83]. Considering the velocity vector u as indicated in Fig. 1b, the face value follows as (where ∼ denotes a flux-limited interpolation) where subscripts U and D denote the upwind and downwind cells, ξ f is the flux limiter and L f = s f /|r U f | is a geometry coefficient. Dependent on the studied flow, the first-order upwind scheme (ξ f = 0), the central differencing scheme (ξ f = 1) or the Minmod scheme [84] are applied for the simulations presented as part of this study, but other TVD schemes could also be applied.

Temporal discretisation
The First-Order Backward Euler scheme or the Second-Order Backward Euler scheme are applied for the discretisation of the transient terms of the governing flow equations. The First-Order Backward Euler scheme is readily given for cell P aŝ with t 1 being the current time-step, superscript (t − t 1 ) denotes values of the previous time-level and V P is the volume of mesh cell P . Assuming a varying time-step is applied to solve the governing equations, the Second-Order Backward Euler scheme is given as [85] In what follows, the discretised governing equations are presented using the First-Order Backward Euler scheme, although the Second-Order Backward Euler scheme is also considered in the results presented in Section 7. Irrespective of the chosen scheme, for consistency all transient terms of the governing equations, ∂ρu j /∂t in Eq. (1), ∂ρ/∂t in Eq. (2), as well as ∂ρh/∂t and ∂ p/∂t in Eq. (3), are discretised with the same scheme.

Advecting velocity
A momentum-weighted interpolation (MWI) method is applied to evaluate the advecting velocity ϑ f = u f n f at cell faces f , with ϑ f taking the role of flux-velocity in the discretised advection terms of the governing equations. MWI emulates a staggered variable arrangement by introducing a cell-to-cell pressure coupling through a low-pass filter acting on the third derivative of pressure [51,82,86], which avoids pressure-velocity decoupling as a result of the collocated variable arrangement and provides a robust pressure-velocity coupling for incompressible and low-Mach number flows [76], while preserving the second-order accuracy of the finite-volume method [51,82].
Following previous studies, the definition of the advecting velocity includes various modifications to the MWI formulation originally introduced by Rhie and Chow [87], to account for non-orthogonal meshes [88][89][90], as well as the large density ratios occurring in interfacial flows and the transient nature of the considered problems [90]. Based on the work of Denner and van Wachem [90], the advecting velocity ϑ f = u f n f at face f , see Fig. 1a, is given as (20) The coefficient d f is defined aŝ where e P and e Q are the sum of the coefficients of the primary variable u arising from the advection terms (and, if viscous flows are considered, the shear stress terms) of the momentum equations (see [90] for a detailed derivation). The face density ρ * f is interpolated by a harmonic average, which is necessary for a consistent definition of the coefficient of the pressure term as well as the efficacy of the density weighting [90]. This density weighting of the pressure gradients, with which robust incompressible interfacial flow simulations with density ratios of up to 10 24 [90,91] were previously presented, stabilises the numerical solution of interfacial flows with large density ratios.

Discretised governing equations
Applying the spatial and temporal discretisation techniques described above to the integral formulation of the continuity equation, Eq. (12), the discretised continuity equation for cell P follows as with the advecting velocity ϑ f given by Eq. (20) and superscript (n + 1) denotes variables that are solved for implicitly.
Although the discretised continuity equation is formulated conservative in mass, it is treated as an equation for pressure.
Hence, based on Eq. (4), the implicit density ρ (n+1) P at cell centre P is given by the pressure-implicit formulation where the fluid properties γ , and R are defined based on the most recent available interface position, as explained in Section 5, and the temperature T is deferred, as further explained in Section 6. The advection term of Eq. (23) is linearised with a Newton linearisation, where superscript (n) denotes values from the previous nonlinear iteration, with the implicit formulation of the advecting velocity ϑ (n+1) f given as Applying the same discretisation principles to the momentum and energy equations, the discretised momentum equations (1) become (27) and the discretised energy equation (3) is given as Both the transient terms and the advection terms of the momentum and energy equations are linearised with a Newton linearisation. With primary solution variable χ , i.e. velocity u j in the momentum equations (27) and specific total enthalpy h in the energy equation (28), the transient terms are linearised as where ρ (n+1) P is given by Eq. (24), and the advection terms are linearised with respect to the density and the advecting velocity, following as where ϑ (n+1) f is given by Eq. (26). Both the spatial pressure derivative on the right-hand side of Eq. (27) and the transient pressure derivative on the right-hand side of Eq. (28) are treated implicitly.
The applied Newton linearisation promotes a smooth transition from elliptic/parabolic to hyperbolic behaviour of the governing equations in different Mach number regimes [50,92,93], which is widely understood to be of particular importance for the continuity equation [50,61,76], where it also provides the required pressure-velocity coupling at low Mach numbers [76,90]. Furthermore, the Newton linearisation of all governing equations facilitates an implicit contribution of active flow-dependent variables, i.e. the density and the advecting velocity, in the context of the presented fully-coupled pressure-based algorithm, which improves the performance and stability of the solution algorithm.

A note on accuracy and conservation
The accuracy of both the underpinning finite-volume formulation (Section 3.1) and the MWI used to evaluate the advecting velocity (Section 3.4) is second order [51,82]. Thus, applying a second-order TVD scheme, such as the Minmod scheme considered in this study, for the discretisation of the advection terms and the Second-Order Backward Euler scheme for the discretisation of the transient terms, the discretised governing equations are overall second-order accurate. Moreover, because the velocity u is treated separately from the density ρ, rather than treating volumetric momentum ρu as a single quantity, a second-order discretisation can be applied to ũ f at contact discontinuities (including fluid interfaces), where the density is discontinuous and ρ f requires a low-order discretisation to avoid spurious oscillations [94]. Although Ham and Iaccarino [89] identified a numerical dissipation of kinetic energy associated with the low-pass pressure filter of the MWI, this numerical dissipation has a marginal magnitude and is proportional to x 2 . As a consequence, the numerical dissipation of momentum and energy is negligible, since viscous stresses and heat conduction are neglected in this study, except at nonlinear waves where numerical dissipation is introduced by the TVD scheme to regulate spurious oscillations.
Although the governing flow equations are solved for p, u and h, the governing equations are discretised in conservative form for mass, Eq. (23), momentum, Eq. (27), and energy, Eq. (28). The density is not a solution variable but is defined by the applied EOS as a function of pressure, temperature and the fluid properties. As demonstrated by Van Doormaal et al. [61] in the context of a pressure-based algorithm, choosing primitive variables (h is not a primitive variable but is treated in the same way) as primary solution variables, instead of the conserved variables, in a fully-conservative formulation of the governing conservation laws does not affect the conservation properties, if a consistent and appropriate discretisation is applied (see also [50,68,74,76] for further examples). To this end, the continuity equation acts as a constraint on the pressure field and results in velocities, through the coupling with the momentum equations, and densities, via the applied EOS, that satisfy the conservation of mass in all Mach number regimes [61]. Furthermore, the same discretisation of ρ f and the same advecting velocity ϑ f are applied in the discretised governing equations, ensuring a consistent formulation of the fluxes.
Consequently, when the system of discretised nonlinear governing equations is converged, the conservative form of the governing conservation laws is satisfied on the discrete level within a predefined solution tolerance. To achieve convergence of the nonlinear system of governing equations, the iterative solution procedure described in Section 6 is applied.

Compressive VOF method
A consistent and precise transport of the colour function ψ , which represents the bulk phases and the fluid interface, is critical for the overall accuracy of the numerical algorithm. In order for the applied fluxes to be consistent, the advection term of the advection equation (10) is reformulated using the chain rule, an approach previously applied in [31,95]. Equation (10) then becomes Now both spatial derivatives of velocity can be consistently discretised with the same advecting velocity ϑ f as the discretised governing equations in Section 3.5.

Discretisation
The compressive VOF method introduced by Denner and van Wachem [77], which does not require an explicit interface reconstruction and has been shown to conserve the volume of incompressible interfacial flows on structured and unstructured meshes accurately, is applied to advect the colour function field with Eq. (32). In the applied finite-volume framework, the colour function in cell P is given as thus representing the discrete volume fraction of the fluid occupying subdomain b , as defined in Eq. (8). Applying the Crank-Nicolson scheme to discretise the transient term [77,96], the discretised advection equation (32) follows as where t ψ is the time-step applied to advect the colour function field, the advecting velocity ϑ f is given by Eq. (20) and the material-dependent compressibility factor K is given by Eq. (11). The face value ψ f is interpolated using the CICSAM scheme [96], including an implicit correction of mesh skewness [77]. The CICSAM scheme is able to advect sharp interfaces with very high accuracy, as for instance demonstrated in [77,96,97], by taking into account the orientation of the interface and, contrary to other algebraic interface advection schemes, the available amount of colour function in the upwind cell that can be advected into the downwind cell over the applied time-step. Accounting for the available flux-volume, however, comes at the cost of a very small time-step t ψ to maintain a sharp interface [97], and a dual time-stepping method is applied to optimise the computational performance [77], as further detailed in Section 6. Although a compressive VOF method is applied in this study, the proposed pressure-based algorithm is not limited to a specific method for the representation and transport of the interface. The applied compressive VOF method is chosen for its versatility, as it is also applicable on unstructured meshes. Any other suitable interface capturing method, such as the class of algebraic THINC schemes for compressible flows [49,98], or interface tracking method, e.g. front-tracking methods [41,48], may equally be used in conjunction with the proposed algorithm.

Validation
The compressive VOF method [77] applied as part of this study has previously been successfully used for a variety of different incompressible, i.e. solenoidal (∇ · u = 0), flows and has thereby been shown to be in excellent agreement with experimental data [77,99] and analytical solutions [91,100]. Contrary to incompressible flows, however, the applicability and accuracy of the modified method given above for compressible flows with ∇ · u = 0 is yet to be demonstrated.
Following the work of Raessi et al. [101], a circular and a spherical interface with initial radius r 0 = 0.25 m are simulated in a square and cubical domain, respectively, with edge length 1 m. In the two-dimensional (2D) case, the centre of the circular interface is situated at x c = y c = 0.5 m and the velocity field is prescribed as  In the three-dimensional (3D) case, the centre of the spherical interface is situated at x c = y c = z c = 0.5 m and the velocity field is prescribed as The velocity field is, hence, not divergence free (∇ · u = 0) and is pointing to the centre of the circular or spherical interface with |u| = 1, causing the interface to shrink and the fluid mass enclosed by the interface to reduce. Note that the flow field is fixed and the result is, therefore, independent of the discretisation of the governing flow equations. The change in area (2D) and volume (3D) is predicted accurately and the computational errors converge with second order in both 2D and 3D, demonstrating that the modified compressive VOF method presented above is suitable for compressible flows.

Acoustically-conservative interface discretisation
Recall the earlier observation, discussed in Section 2, that a fluid interface is a contact discontinuity. At a contact discontinuity, velocity and pressure are continuous, while density (and all variables that depend on density) is discontinuous [32,102]. In a compressible single-phase flow, the contact discontinuity is a characteristic of the solution of the governing conservation laws and represents a weak solution to the discretised governing equations; a prominent example is a contact discontinuity in a shock tube. However, the discretised governing equations do not account for the contact discontinuity associated with an interface separating two immiscible fluids.
This problem can be illustrated by considering the simple example sketched in Fig  momentum equation (1) for this one-dimensional problem at cell P of the equidistant mesh shown in Fig. 3, the advection term follows, using standard finite differences with linear interpolation, as and, with a time-step t corresponding to a Courant number of Co = u t/ x = 0.3, the transient term is Hence, in order to satisfy the discretised momentum equation (27), the jump in density at the fluid interface causes an unphysical pressure gradient and corresponding acceleration of the flow [2,14,16]. Thus, without modifications to the discretisation to account for the change in fluid properties at the fluid interface, the discretised governing equations as presented in Section 3.5 yield discontinuous and typically oscillatory velocity and pressure values. The underpinning principle of the proposed acoustically-conservative interface discretisation (ACID) is based on a conservative discretisation of the governing equations in each finite-volume stencil, with the aim of applying the discretisation derived for a single-phase flow in Section 3 and obtaining consistently defined fluid properties in the interface region where 0 < ψ < 1, while still respecting the contact discontinuity associated with the fluid interface. For the numerical framework applied in this study, the finite-volume stencil of a given cell P , schematically shown in Fig. 4a, contains all face-neighbour cells of cell P , as well as cell P itself. In order to satisfy the assumption of a single-phase flow in the finite-volume stencil of a given cell P for the purpose of discretising the governing equations, all cells in the finite-volume stencil of cell P are assigned the colour function value of cell P . This is schematically illustrated in Fig. 4; the colour function values in cell (i, j) and its neighbourhood, shown in Fig. 4a, are taken to be Fig. 4b.
Thus, the colour function is kept piecewise constant in the entire finite-volume stencil, which enables the application of the fully-conservative discretisation scheme presented in Section 3, identical to the one applied for single-phase flows. The relevant thermodynamic properties that are discontinuous at the interface, i.e. density and enthalpy, are then evaluated and discretised based on this piecewise-constant colour function field, assuming mechanical and thermal equilibrium at the interface. Note that the two-dimensional equidistant Cartesian mesh in Fig. 4 is chosen for illustration purposes; the application of this approach on unstructured and/or three-dimensional meshes is straightforward.

Density treatment
In the proposed pressure-based algorithm, density ρ is not a solution variable. Instead, density is evaluated by a linear interpolation of the partial densities of the bulk phases, given as with ρ a and ρ b defined by the applied EOS, Eq. (4). This linear interpolation of the partial densities with respect to the colour function is necessary for the conservation of mass, momentum and energy, and is equivalent to an isobaric closure assumption [4,35]. For the purpose of treating the density implicitly via Eq. (24), the linear interpolation of partial densities in Eq. (37) is interchangeable with defining the specific gas constant R and the stiffened-gas term γ as respectively.
Assuming the colour function ψ is piecewise-constant throughout the finite-volume stencil of cell P , as described above, the density at face f of cell P is given as with the density at the upwind cell U and downwind cell D evaluated based on the colour function value of cell P as and respectively, where the partial densities are evaluated using Eq. (4).
In order for the discrete density field to be defined consistently, the density values at previous time-levels are evaluated using the same procedure, with and, if required, The density values in the bulk phases are unaffected by this discretisation.

Enthalpy treatment
Similar to the treatment of density, the enthalpy in the interface region is reformulated using ACID. To ensure consistency, the total enthalpy H = ρh at face f is assumed to bẽ because both density ρ and the specific isobaric heat capacity c p are properties of the respective bulk phase. The specific isobaric heat capacity is defined by a density-weighted average as with the partial densities given by Eq. (4), the partial specific isobaric heat capacities given by Eq. (7) and ρ given by Eq. (37). Unlike density, the specific total enthalpy h is a primary solution variable of the applied algorithm and, therefore, cannot simply be replaced by a modified value. Thus, a deferred correction defined by a target face enthalpy ĥ f is applied to the implicitly computed value at every nonlinear iteration. Following the assumption defined in Eq. (45) and the density ρ f given by Eq. (40), the target enthalpy is defined aŝ where ρ U is given by Eq. (41), ρ D is given by Eq. (42), and the specific total enthalpy of the upwind and downwind cells are given as respectively. Based on Eq. (46), the specific isobaric heat capacities c p,U and c p,D are evaluated with colour function ψ P as and Based on the target face enthalpy ĥ f , the deferred enthalpy correction is given as and the advection term of the discretised energy equation, Eq. (28), follows aŝ Away from the interface in the bulk phases, the face enthalpy h f is equal to the target face enthalpy ĥ f , with δh f = 0, and the energy conservation in the bulk phases remains unaffected.
In order for the specific total enthalpy to be defined consistently, the specific isobaric heat capacity of previous timelevels is also defined based on the current colour function ψ P , given as with ρ (t− t 1 ) P given by Eq. (43), and similarly for c . The specific total enthalpy at the previous time-level then follows as and, if required,

Thermodynamics properties
The treatments for density and enthalpy defined above directly influence the thermodynamic properties of the interface region. Of particular interest in this context are the speed of sound and the Rankine-Hugoniot relations, as they represent fundamental solutions to the conservation laws governing compressible flows.
With the density ρ given by Eq. (37) and the specific isobaric heat capacity c p given by Eq. (46), the speed of sound is defined based on Eq. (5) as follows from the isobaric closure assumption applied to the density, see Eq. (38), since R = (γ − 1)c v . Thus, the speed of sound is given by a uniquely defined average of the speeds of sound in the bulk phases, a condition that Allaire et al. [4] associated with a well-possedness of the model governing equations. The accurate prediction of the speed of sound by Eq. (57) in the bulk phases and in the interface region is demonstrated in Section 7.3.1.
Considering a flow with a shock wave, the conditions in the post-shock state (I) can be related to the pre-shock state (II) by the Rankine-Hugoniot relations. Based on the Rankine-Hugoniot relations of the stiffened-gas model [33,59] and assuming the pre-shock region is stationary (u II = 0), the pressure ratio for a given Mach number M s = u s /a II of the shock wave, with u s being the shock speed, is where ˆ is the effective material-dependent pressure constant of the stiffened-gas model, which follows from Eq.
with ρ given by Eq. (37), c p given by Eq. (46), and (γ − 1) and γ given by Eq. (58). For an ideal gas, ˆ = 0. The density ratio across the shock wave is then given as with which the post-shock velocity u I is defined as The post-shock temperature T I readily follows from density ρ I and pressure p I via the applied EOS. Hence, the propagation of shock waves is uniquely defined in each bulk phase as well as in the interface region, and the accuracy of the jump relations given in Eqs. (59)-(62) is demonstrated in Section 7.4.1.

Some observations
ACID is formulated using a general finite-volume discretisation and is independent of the numerical schemes applied to discretise the governing equations or the method used to advect the interface (e.g. VOF, level-set or front-tracking). No Riemann solver is required to compute the fluxes, no wave patterns have to be identified and no a priori assumptions about the wave structure have to be made, as for instance required when applying a classical Godunov scheme [32]. The proposed method assumes that the fluid properties are piecewise constant at the interface, i.e. the fluid properties are formally first-order accurate, and the corresponding error vanishes as x → 0. This is consistent with the piecewise-constant definition of the fluid properties typically applied in VOF, level-set and front-tracking methods [103] and does not affect the local first-order accuracy of the discontinuity in fluid properties [104].
The discretisation of the governing equations remains formally fully conservative, as a set of governing equations in conservative form is solved for each mesh cell. A conservative discretisation of the governing equations is regarded as a prerequisite for the accurate prediction of both the speed and strength of shock waves and for satisfying the Rankine-Hugoniot relations [8,9]. Issues associated with a discontinuous change of fluid properties at the interface are circumvented with ACID, while retaining the information associated with compression and expansion waves, thus allowing acoustic waves, shock waves and rarefaction fans to interact with the interface without artificial (numerical) obstruction. The assumed mechanical and thermal equilibrium at the interface can lead to an unphysical flux in the interface region when both bulk phases have different temperature and different pressure. However, in practice, this spurious flux is typically small and appears to be no problem for dynamic interfacial flows, as observed in the results presented below. Applying a TVD scheme, the discretisation of velocity at the interface remains high-order, which is important when capturing instabilities at the interface (e.g. the Richtmyer-Meshkov instability).
The proposed method shares conceptual similarities with the GFM [40] and with flux-splitting algorithms, such as the one of Abgrall and Karni [15]. GFM and ACID are both one-fluid models, since some form of a single fluid is assumed when discretising the governing equations. GFM copies variables required for the discretisation to ghost cells on "the other side" of the interface and assumes only one of the phases occupies a given mesh cell. ACID, on the other hand, defines the fluid properties in the entire computational stencil based on the value of the interface indicator function taken from the discretised cell, without making explicit assumptions about which phase, or combination of phases, is present in the cell. Similar to flux-splitting algorithms, ACID results in asymmetric fluxes at the cell-faces, i.e. a different flux is applied to a given face dependent on which of its adjacent cells is discretised. To clarify, using ACID, the advecting velocity ϑ f is symmetric at face f , but the density ρ f and the specific isobaric heat capacity c p, f may be asymmetric in the interface region. While ACID evaluates the fluxes based on the advecting velocity ϑ f and the fluid properties only, previously presented flux-splitting algorithms, e.g. [15,45], rely on a Riemann solver to evaluate the fluxes. Furthermore, ACID also adjusts the values at previous time-levels to consistently account for the changes in fluid properties at the interface.
In preliminary studies, the proposed algorithm exhibited convergence issues in some cases where a perfect-gas fluid (e.g. air) and a stiffened-gas fluid (e.g. water) interact, if the Minmod scheme (or other TVD schemes) was applied. Instead of a monotonically reducing residual of the discretised system of equations, oscillations about a constant finite value larger than the predefined solver tolerance η were obtained. After careful inspection, this behaviour is attributed to the inherent nonlinearity of TVD schemes [83], because the definition of the flux limiter depends on the advected scalar field itself. As a result of the typically large differences in density ρ and specific total enthalpy h of perfect and stiffened gases, such as air and water, even small changes of the TVD flux limiter ξ f can have a significant impact on the solution and its convergence.
Since both density ρ and specific total enthalpy h are discontinuous at the interface, TVD schemes revert to the first-order upwind scheme, with the flux limiter ξ f → 0, to avoid an oscillatory solution. Enforcing explicitly the application of the first-order upwind scheme (ξ f = 0) at faces in the interface region, where |ψ P − ψ Q | > η, for the advection of density ρ and specific total enthalpy h is an effective remedy for the convergence issue associated with TVD schemes.

Solution procedure
The linear system, Aφ = b, of the governing equations for a three-dimensional flow discretised as described in Section 3, including the modifications at the interface proposed in Section 5, on a computational mesh with N cells is given as where · denotes the L 2 -norm and η is the predefined solution tolerance.
An inexact Newton method [108] is applied to account for the nonlinearity of the governing equations, with nonlinear iterations in which the linearised governing equations are solved and the deferred variables are subsequently updated. A flow chart illustrating this iterative solution procedure is shown in Fig. 5. Following the work of Xiao et al. [76], the linear system of equations (63) is solved in each nonlinear iteration by assuming the flow is barotropic, i.e. density is only a function of the pressure (but not the temperature). Note that the flow is not considered to be isothermal, but only the update of density neglects changes in temperature in the barotropic loop. When the L 2 -norm of the residual vector δ of the linear system of equations (63) satisfies are both satisfied simultaneously. This solution procedure was shown to be robust for single-phase flows at all speeds and without underrelaxation of the governing equations [76].  method is applied to satisfy the stringent time-step constraints of the CICSAM scheme [77,97] at an only moderate increase of computational time, with a different time-step t ψ applied to solve the VOF advection equation than the time-step t 1 applied to solve the equations governing the fluid flow. In order for the advection of the colour function to be consistent, the fluid time-step t 1 has to be equal to or an integer multiple of the VOF time-step t ψ . For all presented simulations, the Courant number associated with t ψ is Co ψ = t ψ |u|/ x ≤ 0.05.

Results
In order to verify and validate the proposed algorithm, the results for different representative test cases are presented: interface advection in Section 7.1; the conservation of momentum, mass and energy in Section 7.2; the propagation, reflection and transmission of acoustic waves in Section 7.3; the propagation and interface interaction of shock waves in Section 7.4; shock tube problems in Section 7.5; the interaction of shock waves with bubbles and drops in Sections 7.6-7.8. The presented results focus on the propagation of acoustic waves, shock waves and rarefaction fans in interfacial flows. Results for single-phase flows are only briefly discussed to demonstrate the accurate prediction of acoustic phenomena and shock waves by the proposed algorithm, since a broad variety of results obtained with the applied numerical framework for single-phase flows at all speeds and on unstructured meshes has been published by Xiao et al. [76].   at the prescribed velocity u 0 . Although this might appear to be trivial, this has been notoriously difficult to achieve with conservative methods [79].

Interface advection with constant velocity
The interface is represented by a stepwise change in colour function ψ , initially located at x = 0.1 m. Fig. 6 shows the profiles of velocity u, pressure difference p = |∂ p/∂ x| x and density ρ at time t = 0.7 s, with and without ACID. The interface and the associated discontinuity of the fluid properties do not affect the flow if ACID is applied, and the interface is advected with the correct speed of u = u 0 . The equilibrium of pressure p and velocity u as well as the stepwise change of density ρ at the interface are retained, contrary to previously published algorithms, e.g. [31]. Note that the discontinuity in colour function and density at the interface does not numerically diffuse, which is a common issue of diffuse interface methods methods for compressible flows [109][110][111]. However, if ACID is not applied, the flow field is clearly affected in an unphysical way, with large changes in velocity and pressure, and the interface is not advected with the correct speed, as evident in Fig. 6.
As a second test-case, a smooth interface with a linear variation of ψ from ψ = 1 to ψ = 0 over an interface thickness of 0.1 m (50 cells) is simulated, with the interface initially located between x = 0.1 m and x = 0.2 m. Such cases, with cells partially occupied by both fluids, present a particular difficulty for density-based algorithms with approximate Riemann solvers. The profile of the velocity difference u = |∂u/∂x| x, pressure difference p = |∂ p/∂ x| x and density ρ at t = 0.7 s are shown in Fig. 7. The differences of pressure and velocity are negligible in the entire domain, while the density profile and, consequently, the position of the interface are predicted accurately, demonstrating that a smooth interface does not present a problem for the presented algorithm. As in the case of the sharp interface discussed above, the interface (i.e. the discontinuity in colour function and fluid properties) does not artificially diffuse.

Conservation of momentum, mass and energy
In order to avoid spurious oscillations at the interface, the conservation of momentum, mass and energy at the interface is often relaxed or sacrificed [2,5,14,23]. As discussed in Section 5, the proposed discretisation method ACID applies a correction to the discretisation at the interface to account for the change in fluid properties. The fully-conservative formulation and discretisation of the governing flow equations is, nevertheless, retained in the finite-volume discretisation stencil of each mesh cell. Furthermore, the compressive VOF method applied in the proposed algorithm was previously shown to conserve the volume of each phase in incompressible flows [77].
The advection of a circular bubble with constant velocity is simulated. The fluid properties are identical to the onedimensional interface advection in Section 7.1. The computational domain is represented by an equidistant Cartesian mesh where superscript (0) denotes values at reference time t = 0. Fig. 8 shows the relative conservation errors ε of momentum, mass and energy for the entire domain and for the bubble, as a function of time. The global (domain) conservation errors are negligible and the governing conservation laws are satisfied accurately. With respect to the bubble, the conservation errors are still insignificant but up to one order of magnitude higher than the corresponding global errors. This can be attributed to the volume conservation error of the applied VOF method, shown in Fig. 8d, which is the leading error of the conservation of the bubble.

Acoustic waves
The propagation of acoustic waves in a one-dimensional domain is simulated to study the capabilities of the proposed algorithm to predict acoustic effects. The accurate prediction of the formation and propagation of acoustic waves is an important feature for the simulation of compressible flows and was previously found to be a particularly challenging problem [76,112], because even small inconsistencies in the discretisation or a lack of convergence can lead to a visible change in the amplitude and speed of the waves. In the presented simulations, the acoustic waves are generated at the domain-inlet by a sinusoidal velocity perturbation with (small) amplitude u 0 . For small perturbations to the flow, with density amplitude ρ ρ 0 and velocity amplitude u a 0 , the resulting wave is a sound wave. The propagation speed of such a small-amplitude acoustic wave is the speed of sound a 0 and, according to linear acoustic theory, the resulting amplitudes of the density waves ρ 0 and the pressure waves p 0 are [102] Five different materials are considered, for which the fluid properties are given in Table 1. In all cases, the computational domain is initialised with uniform velocity u 0 = 1.0 m s −1 and pressure p 0 = 10 5 Pa. The central differencing scheme is applied for the presented single-phase flows and the Minmod scheme is applied for the interfacial flows, while the secondorder backward Euler scheme is applied to discretise the transient terms in both single-phase and interfacial flows. At the domain-inlet, pressure and temperature are extrapolated from the closest cell centre, and at the domain-outlet, a zerogradient condition is specified for all variables. The propagation of acoustic waves in single-fluid flows, either by considering only a single fluid or a constant mixture of two fluids (emulating the interface region between two fluids), is discussed in Table 1 Fluid properties for the propagation of acoustic waves. The heat capacity ratio γ and pressure constant of water and copper are taken from [6].

Propagation of acoustic waves
Acoustic waves in air and water, as well as in the interface region of air-helium, air-water and water-copper flows are simulated to demonstrate the accurate prediction of the speed of sound and pressure waves in the bulk phases and the interface region. The computational domain has a length of 1 m, represented by a Cartesian mesh with mesh spacing x = 2 × 10 −3 m, and the applied time-step corresponds to a Courant number of Co = a 0 t/ x = 0.44 − 0.52. The velocity at the domain-inlet is u in = u 0 + u 0 sin(2π f t), with frequency f and amplitude u 0 = 0.01 u 0 .  In summary, the propagation of acoustic waves in gases, liquids and solids is predicted accurately, especially considering that the density wave typically has an amplitude that is five orders of magnitude smaller than the ambient density. The pressure and density waves exhibit a constant amplitude as they propagate downstream, which is expected for an inviscid flow simulated with second-order discretisation schemes, as discussed in Section 3.6. In fact, at the end of a domain with a length of 50 m ( x = 2 × 10 −3 m, Co = 0.5), the amplitude of velocity, pressure and density of an acoustic wave in air ( f = 500 s −1 ) are all > 99.75% of their initial value at the domain inlet. The definition of the speed of sound given in Eq. (57) is found to be accurate in the bulk phases and throughout the interface region (i.e. where a mixture of two bulk phases is present), demonstrating that ACID yields a thermodynamically-consistent discretisation.

Reflection and transmission at fluid interfaces
The propagation of a single acoustic wave in a helium-air flow, an argon-air flow and an air-water flow is simulated in a one-dimensional domain with mesh spacing x = 2 × 10 −3 m and with a time-step that corresponds to Co = a 0 t/ x = 0.48. The single acoustic wave is initiated by the inlet-velocity with u 0 = 0.02 u 0 . The frequency of the acoustic wave is f = 5000 s −1 in the helium-air flow and the air-water flow, and f = 2000 s −1 in the argon-air flow.
Assuming the incident acoustic wave travels from left to right, a part of the pressure wave is transmitted to the right phase when the incident wave reaches the fluid interface, while the remaining part of the incident wave is reflected in the left phase. Based on linear acoustic theory, as the incident wave with its pressure amplitude p incid.
with p trans.
L,0 , and Z L = ρ L a L and Z R = ρ R a R are the acoustic impedance of the left and right phase, respectively. The pressure amplitude of the reflected wave, thus, follows as p refl.
In Fig. 12, the pressure profiles of the acoustic wave in the helium-air flow and the argon-air flow are shown, before and after the interaction of the incident wave with the fluid interface. In the helium-air case, Fig. 12a, the pressure amplitude p incid. He = 3.306 Pa of the incident wave is predicted accurately compared to linear acoustic theory ( p incid. He,0 = 3.307 Pa).   Air,0 = 4.686 Pa). In the argon-air case, Fig. 12b, the pressure amplitudes of the incident, reflected and transmitted acoustic waves are also in very good agreement with linear acoustic theory. The ratio of pressure amplitudes is predicted as p trans.
Air / p refl. Ar = −5.919, compared to the theoretical value of p trans. Air,0 / p refl. Ar,0 = −5.871. Since the presented flows of two interacting ideal gases with acoustic waves is isentropic, the change in specific entropy is s = 0, where The computed change in entropy as the acoustic wave propagates in the bulk phases of the helium-air flow is clearly negligible, as seen in the inset of Fig. 13. An increase of | s| can be observed in Fig. 13 when the acoustic wave interacts with the helium-air interface (at t/t ≈ 1, with t the time required for the peak of the acoustic wave to reach the interface), a change in specific entropy that can, nevertheless, be regarded as inconsequential given its insignificant magnitude. Fig. 14 shows the pressure profile for the acoustic wave in the air-water flow, alongside the theoretical pressure amplitude and wavelength given by linear acoustic theory. Contrary to the interaction of two perfect gases discussed above, water is described by the stiffened-gas model and, thus, large differences in acoustic impedance and density between the two interacting bulk phases ensue; the acoustic impedance and the density of water are approximately three orders of magnitude larger than for air. The amplitudes of the reflected and transmitted pressure waves are in very good agreement with linear acoustic theory, with the ratio of pressure amplitudes obtained from the simulation being p trans.
Water / p refl. Air = 1.995, compared to the theoretical value of p trans.
The initial velocity in both cases is u 0 = 0.30886 m s −1 .
The result of the first case ( Z = 423.588 Pa s m −1 ) is shown in Fig. 15a for a sinusoidal velocity perturbation with f = 2000 s −1 and u 0 = 0.01u 0 . The observed amplitude of the pressure wave is in excellent agreement with linear acoustic theory in both phases, while the change in wavelength as a result of the different speeds of sound is predicted accurately and is clearly visible. Note that a reflected wave in the left phase would lead to a strong interference with the oncoming incident waves, which would be visible in the pressure profile.
The result of the second case ( Z = 500 Pa s m −1 ) is shown in Fig. 15b for a single sinusoidal wave with f = 5000 s −1 and u 0 = 0.02u 0 . A small reflected wave can be identified in the pressure profile of the left phase at t = 0.9 × 10 −3 s. Nevertheless, this reflected wave has only a minor effect on the transmitted wave in the right phase, which has a pressure amplitude of p trans.

Shock waves
The simulation of the propagation of shock waves poses a particular difficulty for numerical algorithms, because shock waves are discontinuous and because valid solutions of the governing conservation laws (presented in Section 2) are not guaranteed to satisfy the second law of thermodynamics across shock waves [9]. This raises the question whether the proposed algorithm reliably converges to the physically-correct weak solution of the governing conservation laws, which is a prerequisite for the accurate prediction of both the speed and strength of strong shock waves [8,9]. The simulation of strong shock waves in air and water is examined in Section 7.4.1 and the interaction of shock waves in a single-phase flow is studied in Section 7.4.2. The interaction of shock waves with an air-helium interface and an air-water interface are presented in Sections 7.4.3 and 7.4.4, respectively. The interaction of a shock wave with an interface separating to fluids with the same shock impedance is discussed in Section 7.4.5.

Shock wave propagation
The propagation of a shock wave travelling to the right with Mach number M s = u s /a II in air and water, with γ and given in Table 1, in a computational domain with a length of 1 m is simulated. The pre-shock region (II) has pressure  in the post-shock region is B Water,I = ρa 2 = 2.9 × 10 13 Pa (for comparison, B Air,II = 1.4 × 10 5 Pa), which demonstrates the robustness of the proposed numerical algorithm even for strong shock waves (large pressure ratio) and marginally compressible fluids. The smearing of the shock wave reduces with mesh refinement, shown in Fig. 16, and the computed shock wave reproduces the solution of the corresponding Riemann problem with increasing precision. The common intersection point of the pressure profiles for different mesh resolutions observed in Fig. 16, which coincides with the intersection point with the corresponding Riemann solution, demonstrates that the speed of the shock wave is predicted accurately irrespective of the mesh resolution.
In Section 5.3, the Rankine-Hugoniot relations for the interface region are presented. These relations are tested using a constant air-water mixture, with the properties of the individual (pure) phases being the same as above. The domain is initialised based on the Rankine-Hugoniot relations given in Section 5.3. As for the single-phase cases discussed in the previous paragraph, the applied time-step corresponds to Co = 0.5, the shock wave is initially located at x s,0 = 0.1 m and the simulations are concluded at t end = 0.7 m/u s . Fig. 18 shows the pressure profile of a shock wave with M s = 10 in an air-water flow with ψ ∈ {0.25, 0.50, 0.75}, where ψ = 0 represents air and ψ = 1 represents water. The pressure profiles are in excellent agreement with the theoretical Riemann solution and converge with mesh refinement.
Given the pressure-based formulation of the proposed algorithm, the explicit evaluation of the density may be considered the "Achilles' heel" of the numerical framework. The L 1 -norm of the relative error associated with the density field, given for a computational mesh with N cells as where ρ (comp.) P is the density computed with the presented numerical algorithm at cell P and ρ (exact) P is the corresponding exact density given by the Riemann solution at the centre of cell P , on meshes with different resolutions for the shock wave with M s = 10 is shown in Fig. 19 for air and water single-phase flows, as well as the air-water mixture with ψ = 0.75. For  all considered cases, the density error is of similar magnitude and converges with first order, as imposed by the applied monotone discretisation schemes and as expected for an oscillation-free numerical simulation of a shock wave [94].
The error convergence under mesh refinement, supported by the improved shock resolution observed qualitatively in Figs. 16 and 18, is expected from a discretely conservative algorithm and, by virtue of the Lax-Wendroff theorem [113], suggests convergence to the physically-correct weak solution of the governing conservation laws [8]. Note that a better shock resolution can, of course, also be obtained by applying a TVD scheme or by reducing the applied time-step, which is not shown here in the interest of conciseness, but was previously demonstrated by Xiao et al. [76] using the single-phase framework the proposed algorithm is based on.

Shock wave interaction
The interaction of two shock waves in a closed, one-dimensional domain, as previously studied by Woodward and Colella [114], is simulated. The one-dimensional domain is 1 m in length and occupied by a gas with γ = 1.  [114]. Note that the results of Woodward and Colella [114] were obtained on a mesh with 3096 cells, which was adaptively refined around flow discontinuities by factor 8, yielding an effective resolution equivalent to an equidistant mesh with 24768 cells. The results obtained on the finer mesh are in very good agreement with the reference results and a significant improvement of the solution, in particular for the density profile, is observed when the mesh spacing is reduced.

Air-helium interface
The interaction of a planar shock wave with a planar air-helium interface is simulated. The shock wave is initially located at x s,0 = 0.05 m and travels in the left phase, which is assumed to be air, with speed u s towards the air-helium interface that is initially situated at x ,0 = 0.15 m. The shock wave separates the post-shock region I and the pre-shock region II, which are initialised with u I = 125.65 m s −1 , p I = 1.59060 × 10 5 Pa, T I = 402.67 K,   The heat capacity ratio and the specific heat capacity at constant volume of air are γ Air = 1.4 and c v,Air = 720 J kg −1 K −1 , respectively. Hence, the density of the air in the pre-shock region is ρ II,Air = 1 kg m −3 , resulting in a speed of sound of a II,Air = 376.65 m s −1 . Following the experiments of Haas and Sturtevant [115] and the numerical study of Quirk and Karni [116] of the interaction between a shock and a helium bubble, the helium phase is assumed to be contaminated with 28% air by mass, with a heat capacity ratio of γ He = 1.648 and a specific heat capacity at constant volume of c v,He = According to the corresponding Riemann solution, the shock wave is observed to travel with a speed of u s = 459.50 m s −1 , which corresponds to M s = u s /a II,air = 1.22. Fig. 22 shows the profiles of density, pressure and Mach number, at t = 2 × 10 −4 s after the shock has interacted with the interface. A rarefaction fan is reflected in the air phase and a shock wave transmitted in the helium phase. All computed variables compare very well with the corresponding Riemann solution, including the strength and speed of the shock wave transmitted to the helium phase, which is propagating at u s = 1079.90 m s −1 (M s = 1.13). The position of the contact discontinuity (coinciding with the fluid interface) and the associated change in density and Mach number are also predicted accurately, as observed in Fig. 22.

Air-water interface
The interaction of a shock wave with M s = 10 with an air-water interface is simulated. Using the heat capacity ratio γ and pressure constant of air and water given in Table 1, the post-shock region I and pre-shock region II are initialised with u I = 2869.3 m s −1 , p I = 1.165 × 10 7 Pa, ρ air,I = 6.614 kg m −3 , u II = 0 m s −1 , p II = 10 5 Pa, ρ air,II = 1.157 kg m −3 .
The applied one-dimensional domain is 1 m in length and represented with a mesh of 1000 equidistant cells ( x = 10 −3 m) and a time-step of t = 10 −8 s. The shock wave and the air-water interface are initially situated at x s,0 = 0.25 m and x ,0 = 0.50 m, respectively. Fig. 23 shows the profiles of density, pressure and velocity, at t = 2.78 × 10 −4 s after the shock wave has interacted with the interface. The density profile computed with the proposed algorithm is in excellent agreement with the exact Riemann solution, despite the large density ratio of air and water. The pressure and velocity profiles are also in good agreement with the theoretical solution. However, small discrepancies can be identified in the pressure profile in comparison to the theoretical solution, with a wave forming upstream of the shock. This pressure wave actually originates during the initial contact of the shock wave with the interface, and propagates with the flow. Nevertheless, it does not have a lasting effect on the fidelity of the pressure distribution.

Gas-gas flow with shock impedance matching
The interaction of a shock wave with a fluid interface when both bulk phases have the same shock impedance is considered. As a result of the shock impedance matching, no shock wave or rarefaction fan is reflected at the interface. Cases with shock impedance matching have previously been found to be problematic with GFM-based methods [46,47], where spurious shock waves or rarefaction fans are reflected at the interface. For a non-reflective shock-interface interaction, the fluid properties and the pressure ratio across the shock wave have to satisfy [46] (γ L − 1)ρ II,L + (γ L + 1)ρ II,L Assuming the same setup as in Section 7.4.3 for the shock wave with M s = 1.22 in an air-helium flow, γ L = γ air , ρ II,L = ρ II,air and the right phase has γ R = 1.648, the specific isochoric heat capacity of the right phase follows from Eq. (74) as c v,R = 512.41 J kg −1 K −1 . Fig. 24 shows the profiles of density, pressure and Mach number, at t = 2 × 10 −4 s after the shock-interface interaction. All variables are predicted accurately compared to the corresponding Riemann problem, without spurious reflections in the left phase or spurious oscillations at the interface. As observed for the previously discussed cases involving shock waves, the speed and position of the shock wave are computed with very high accuracy.

Shock tubes
Shock tubes are extensively used to test and scrutinise numerical frameworks for compressible flows, because an exact reference solution based on the associated Riemann problem is available and because they typically feature all three primary wave types (shock wave, rarefaction fan and contact discontinuity).

Gas-gas shock tube
A subsonic and a transonic shock tube with two-phase flow are simulated to assess the handling of discontinuous data and the prediction of shock waves, rarefaction fans and contact discontinuities by the presented algorithm. In both cases the discontinuity of the fluid states is initially located in the middle of the one-dimensional domain with a length of 1 m, which is represented with 400 equidistant mesh cells ( x = 2.5 × 10 −3 s), with a time-step that corresponds to Co = a R t/ x ≤ 0. 27 The profiles of density, pressure, velocity, Mach number and acoustic impedance are shown in Fig. 25 for the subsonic shock tube at t = 8 × 10 −4 s and in Fig. 26 for the transonic shock tube at t = 6 × 10 −4 s, alongside the corresponding theoretical Riemann solution. All shown quantities are in excellent agreement with the theoretical Riemann solution in both cases, irrespective of the Mach number regime. The accurate prediction of the acoustic impedance at the fluid interface further supports the notion that ACID conserves the acoustic properties of the two-phase flow, while the accurate prediction of the Mach number suggests that the hydrodynamic-thermodynamic coupling of the proposed algorithm is correct. The correct position of the shock wave, the rarefaction fan and the contact discontinuity in both cases demonstrates that the conservative discretisation of the governing flow equations leads to physically sound numerical results satisfying the Rankine-Hugoniot relations and entropy condition. Furthermore, the computational result converges to the theoretical Riemann solution under mesh refinement, as demonstrated by the density profile of the subsonic shock tube in Figs. 27a and 27b.  The contact discontinuity occurring in the solution of the shock tubes, which coincides with the fluid interface, is a linearly degenerate wave and, thus, represents the primary difficulty and main source of error with respect to the convergence of the applied finite-volume method under mesh refinement [11,104]. As shown rigorously by Banks et al. [104], the convergence of the L 1 -norm of the solution error associated with a linearly degenerate wave for a q-th order advection scheme without compressive limiting is of order q/(q + 1). Hence, for the applied second-order Minmod scheme the order of convergence of the L 1 -norm of density ρ at the contact discontinuity is 2/3 [104]. The L 1 -norm of the density profile, given by Eq. (73), for the subsonic gas-gas shock tube is shown in Fig. 27c. The error of the density profile converges with order 0.657 (based on a regression fit of the data points, with R 2 = 0.9987), which is in very good agreement with the theoretical value for the convergence rate of 2/3. This suggests that the ACID method does not adversely affect the convergence of the numerical framework, considering that the contact discontinuity coincides with the material interface.

Gas-liquid shock tube
A gas-liquid shock tube is simulated, with air as the left phase and water as the right phase.   with the nonlinear change of the speed of sound in the interface region (see Fig. 11b), which has, however, no lasting influence on the overall result of the shock tube.

Two-dimensional shock-bubble interaction
The interaction of a shock wave with M s = 1.22 in air with a circular bubble of helium and of R22 is simulated. These cases have previously been studied experimentally by Haas and Sturtevant [115] and considered in a number of numerical studies to investigate the governing physical mechanisms [116,117] and to test numerical algorithms, e.g. [5,31,41]. The computational setup is schematically illustrated in Fig. 29. Since the problem is symmetric about the centreline of the domain, as indicated in Fig. 29  computations and experiments is also observed at τ = 1.77 with respect to the primary shock wave, the reflected waves immediately following the primary shock wave and the Mach stem. Also note that the density gradient and, consequently, the interface representation stay sharp throughout the simulation, without artificial smearing of the interface. Furthermore, the presented results are also found to be in very good agreement with previous numerical studies, notably [31,41,116,118]. The growth and magnitude of the instabilities forming at the interface, e.g. on the downstream side of the interface in Fig. 30f, exhibit significant differences between various numerical studies (cf. [5,118]) and the "correct" development of the instabilities, although clearly present, cannot be precisely deduced from the experimental images of Haas and Sturtevant [115]. Johnsen and Colonius [118] showed that spurious pressure oscillations at the interface increase these instabilities artificially, leading to more pronounced interface instabilities predicted by methods that do not completely eliminate spurious pressure oscillations.

R22 bubble
The properties of the R22 bubble are taken from the experimental study of Haas and Sturtevant [115], with heat capacity ratio γ R22 = 1.249 and specific isochoric heat capacity c v,R22 = 365 J kg −1 K −1 . The contours of the density gradients and the Mach number are shown for different time instances in Fig. 31. Similar to the results presented for the shock interaction with a helium bubble in the previous section, the simulation results obtained with the proposed pressure-based algorithm are in very good agreement with the experimental results reported by Haas and Sturtevant [115], cf. Fig. 11 in [115] where   wave has passed the bubble, see τ = 0.68 in Fig. 31c and τ = 0.89 in Fig. 31d, it collapses back on itself, leading to a number of reflected waves captured by the simulations, which are also observed in the experiments of Haas and Sturtevant [115]. The presented result also compare well to previous simulation results [110,117,119]. As mentioned in Section 3, the proposed algorithm is applicable on unstructured meshes. Fig. 32 shows the R22 bubble during and after the impact of the shock wave on a triangular mesh with an equivalent mesh spacing (i.e. the mesh spacing of an equidistant Cartesian mesh with the same number of cells) of x = d 0 /293, alongside the results obtained on an equidistant Cartesian mesh with x = d 0 /300. The qualitative features pertaining to the shock wave and its reflections, as well as the interface shape, are in very good agreement. In general, more flow features are resolved in the bulk phases on the Cartesian mesh, while instabilities developing along the interface are more pronounced on the triangular mesh. Despite these small qualitative differences, the density, pressure and velocity profiles obtained on both meshes on a line along the x-axis at y = 0.002 m, shown in Fig. 33 for τ = 0.89, are in excellent agreement.  For the same Cartesian and triangular meshes, the conservation of mass in the entire computational domain (accounting for the additional mass that enters through the domain inlet), shown in Fig. 37, is conserved accurately, with an error < 0.01% on the Cartesian meshes and < 0.03% on the triangular meshes. Similarly, the conservation of mass of the bubble is very good, see Fig. 38, being within 0.1% on the Cartesian meshes and within 1.0% on the triangular meshes.

Influence of the advection differencing scheme
The applied numerical method has a direct influence on the development and evolution of the compression and expansion waves, as well as instabilities at the interface, as briefly mentioned in Section 7.6.1. For instance, Johnsen and      Colonius [118] observed a substantial influence of spurious pressure oscillations at the interface on the onset and evolution of interfacial instabilities. Although spurious pressure oscillations are absent in the results presented above, a clear influence of the applied differencing schemes can be observed for the shock interaction with a helium bubble and a R22 bubble.
Figs. 39 and 40 show the contours of the density gradient and the velocity magnitude for the helium bubble at τ = 4.50 and τ = 7.70, respectively, obtained using the Minmod scheme and the Superbee scheme [84], on a Cartesian mesh with x = d 0 /500. Because of the compressive limiting applied with the Superbee scheme, the compression and expansion waves are resolved sharper. Moreover, interface instabilities form more rapidly with the more compressive Superbee scheme.
A similar behaviour can be observed for the R22 bubble at τ = 0.89 using the first-order upwind scheme, the Minmod scheme and the Superbee scheme, for which the resulting density gradient and velocity magnitude are shown in Fig. 41. Apart from the expected sharper resolution of compression and expansion waves with increasing compression of the applied differencing scheme, the Mach number at the interface differs substantially; while the flow is entirely subsonic with the upwind scheme, the flow inside the bubble becomes transonic with the Minmod scheme and supersonic with the Superbee scheme. Although the temporal evolution of the volume of the R22 bubble is virtually the same irrespective of the applied differencing scheme, as seen in Fig. 42a, the conservation of the bubble mass is strongly affected by the applied differencing scheme, as observed in Fig. 42b. Based on these results, the Minmod scheme appears to provide the most feasible result, although this issue requires a further, more comprehensive investigation to draw firm conclusions about the most appropriate differencing scheme.

Three-dimensional shock-bubble interaction
The interaction of a shock wave with M s = 1.68 in air with a three-dimensional helium bubble is simulated, as previously studied in [117]. Since the problem is axisymmetric with respect to the x-axis, only one quarter of the domain in the y − z plane is simulated. The domain has the dimensions 6r 0 × 4r 0 × 4r 0 , where r 0 = 0.0254 m is the initial radius of the bubble, and is represented by a Cartesian mesh with x = 5 × 10 −4 m (50.8 cells per r 0 ). Following Niederhaus et al. [117], the heat capacity ratio and specific isobaric heat capacity of air and helium are γ Air = 1.399 and c p,Air = 1006.36 J kg −1 K −1 , and    The contours of density and Mach number are shown in Fig. 43 for τ s = t u s,He /r 0 ∈ {1.4, 2.6, 5.2}, where u s,He is the speed of the shock wave in helium, which are the same dimensionless time instances τ s as considered by Niederhaus et al. [117].
The density field, position of the Mach stem, and the position and shape of the primary shock wave and the reflected rarefaction fan computed by the proposed algorithm agree very well with the results reported in [117] for all considered  time instances. The density field at τ = 2.2 is shown in Fig. 44a together with the isocontour ρ = 2.5 kg m −3 , which separates the pre-shock and the post-shock regions next to the bubble, and which follows the Mach stem and the rarefaction fan in the immediate vicinity and behind the bubble. The same result is shown in Fig. 44b but with the isocontour ρ = 1.5 kg m −3 that represents the shock front. Both isocontours in Fig. 44 are clearly axisymmetric, demonstrating that the proposed algorithm captures the three-dimensionality of the solution reliably.

Shock-drop interaction
The interaction of a shock wave with a liquid drop is simulated, with different Mach numbers of the shock wave. Contrary to the cases presented in Sections 7.6 and 7.7 where two perfect-gases interact, in the cases presented below a liquid drop described by a stiffened gas ( p) is surrounded by a perfect gas ( = 0).

Shock-drop interaction with M s = 1.47
The interaction of a shock wave with Mach number M s = 1.47 in air with a circular water drop, as previously studied in [30,[120][121][122], is simulated. The flow is assumed to be symmetric and the computational setup is schematically illustrated in Fig. 45.   Following the properties applied by Meng and Colonius [121], air has the properties γ Air = 1.4, Air = 0 and R Air = 287.08 J kg −1 K −1 , while water is described by γ Water = 6.12, Water = 3.43 × 10 8 Pa and R Water = 7170.23 J kg −1 K −1 .
The computational domain is represented by an equidistant Cartesian mesh with x = d 0 /200 and the applied time-step corresponds to a Courant number of Co = a Water,II t/ x = 0.15. Fig. 46 shows the density gradient and the velocity magnitude of the flow at different times t ∈ {7.74 μs, 16.18 μs, 23.00 μs}, with t = 7.74 μs corresponding to the dimensionless time t * = 0.017 in [121]. The results are qualitatively in very good overall agreement with previously reported numerical results for this particular shock-drop interaction [121,123]. The error in mass conservation ε ρ , shown in Fig. 47, is, however, about one order of magnitude larger than for the shock-bubble cases presented in Section 7.6.2. This increase in mass error extends to both the drop as well as the entire computational domain. Since the proposed framework is nominally mass-conservative, as demonstrated in Section 7.6.2, this error is likely the result of the finite interface thickness introduced by the applied VOF method and the associated misrepresentation of the fluid properties in this interface region, as further studied below.

Shock-drop interaction with M s = 6
The interaction of a shock wave with Mach number M s = 6 in gas with a circular liquid drop, similar to previous work of Shukla et al. [109], is simulated. The computational setup is schematically illustrated in Fig. 48. The post-shock region I and the pre-shock region II are initialised with u I = 5.789 m s −1 , p I = 42.388 Pa, T I = 2.794 × 10 −2 K, u II = 0 m s −1 , p II = 1.013 Pa, T II = 3.518 × 10 −3 K.    Fig. 11 in [109]. The results are overall in very good agreement with the results reported by Shukla et al. [109]. Comparing the results on different meshes at t = 0.34 s, see Figs. 49c and 50, small but clearly visible differences in the shape and structure of the compression and expansion waves are observed. This indicates mesh-dependent differences in the interaction of the waves with the gas-liquid interface, since the mesh-dependence of the speed and position of these types of waves is in general negligible, as demonstrated in Sections 7.4 and 7.5. Because of the large differences in fluid properties, notably the density and heat capacity, in gas-liquid flows, even small differences of the interface position, likely caused by small errors in the interface advection (see Section 4.2), or mesh-dependent thickness of the interface, can have a significant effect on the wave-interface interactions. These differences also manifest in mass-conservation errors with respect to the computational domain and the drop, see Fig. 51. Yet, as expected from errors associated with the applied interface capturing scheme, the mass-conservation error reduces with mesh refinement. Even though qualitative differences are visible in the contour plots, and the mass-conservation errors are not negligible for the results obtained on different meshes, the position and velocity of the centre of mass of the drop and the volume of the drop, shown in Fig. 52, are in good agreement on the different meshes, and converge with mesh refinement.

Conclusions
A pressure-based algorithm for compressible interfacial flows has been presented, based on a fully-coupled finite-volume framework with a collocated variable arrangement that is applicable to unstructured meshes. The governing flow equations, which are discretely conservative in mass, momentum and energy, are solved for velocity, pressure and specific total enthalpy. Contrary to all previously published algorithms for compressible interfacial flows, density is not a solution variable but is computed based on pressure, temperature and the local fluid properties via an equation of state. To this end, the stiffened-gas model has been considered in this study to describe gases, liquids and solids. The proposed algorithm includes an acoustically-conservative interface discretisation (ACID), which performs a local manipulation to the discrete values of density and total enthalpy based on the interface position, and which does not require the solution of a Riemann problem or a priori knowledge of the wave structure to evaluate fluxes.
Results for a variety of representative compressible gas-gas and gas-liquid interfacial flows have been presented, including the propagation of acoustic waves, shock tubes and the interaction of shock waves with one-, two-and threedimensional interfaces. The transmission and reflection of acoustic waves at the interface has been shown to be captured accurately by the presented algorithm, including cases with acoustic impedance matching, a capability not previously reported in the literature. Similarly, the interaction of shock waves with the interface in one-, two-and three-dimensional flows, as well as gas-gas and gas-liquid shock tube problems, have been favourably compared against exact Riemann solutions and previously reported experimental and computational results, especially the temporal evolution of the shock position, of reflected shock waves and rarefaction fans, of the wave structure and the interface shape. In all considered cases, the presented algorithm reliably retains the acoustic properties of the fluids, while the accurate prediction of the Mach number suggests a correct hydrodynamic-thermodynamic coupling. Notably, the speed, position and strength of shock waves are predicted accurately and the Rankine-Hugoniot relations are satisfied even for strong shock waves, in the bulk phases and in the interface region, demonstrating that the proposed algorithm converges to the correct weak solution of the governing conservation laws [8][9][10]. This enforces the notion that, irrespective of small conservation errors associated with the applied interface capturing scheme, the governing equations remain discretely conservative and implies that the proposed algorithm implicitly satisfies the second law of thermodynamics.
Despite the demonstrated success of the presented algorithm, several open questions remain. The EOS of the perfect-gas model or the stiffened-gas model is used in this study to evaluate the density based on pressure and temperature. However, most EOS are formulated in terms of internal energy rather than density. While this is not a problem for the considered gas models, and some real-gas models can readily be formulated in terms of density, such as the Peng-Robinson EOS [124], more complex EOS may not be easily inverted to a density formulation. Furthermore, the proposed interface discretisation method ACID implies a mechanical and thermal equilibrium. Mechanical equilibrium is a reasonable assumption for interfacial flows, since mechanical relaxation is associated with acoustic effects [21,26] and, thus, occurs at a time-scale comparable with the typically applied computational time-steps. The implicit assumption of thermal equilibrium, however, may present a limiting factor when more complex EOS are considered [4] and/or thermal relaxation is governed by diffusion [59], for instance if reactive flows are to be simulated. Neither of these questions have been addressed in the context of pressure-based methods and warrant further study.
In summary, the proposed pressure-based algorithm has been shown to be a promising alternative to traditional densitybased algorithms for the simulation of compressible interfacial flows. The pressure-based formulation of the governing equations facilitates the definition of consistent mixture rules at the interface, including the Rankine-Hugoniot relations, and applies naturally to flows in all Mach number regimes.