An All-Mach Number HLLC-Based Scheme for Multi-Phase Flow with Surface Tension

: This paper presents an all-Mach method for two-phase inviscid ﬂow in the presence of surface tension. A modiﬁed version of the Hartens–Lax–van Leer Contact (HLLC) solver is developed and combined for the ﬁrst time with a widely used volume-of-ﬂuid (VoF) method: the compressive interface capturing scheme for arbitrary meshes (CICSAM). This novel combination yields a scheme with both HLLC shock capturing as well as accurate liquid–gas interface tracking characteristics. It is achieved by reconstructing non-conservative (primitive) variables in a consistent manner to yield both robustness and accuracy. Liquid–gas interface curvature is computed via height functions and the convolution method. We emphasize the use of VoF in the interest of interface accuracy when modelling surface tension effects. The method is validated using a range of test-cases available in the literature. The results show ﬂow features that are in sensible agreement with previous experimental and numerical work. In particular, the use of the HLLC-VoF combination leads to a sharp volume fraction and energy ﬁeld with improved accuracy. plat-form. The developed solver is rigorously assessed via application to a range of benchmark test-cases from 1-D to 2-D.


Introduction
High-speed multi-phase compressible flow induced by blast or shock waves is of interest to both basic science and engineering [1][2][3][4][5][6][7][8]. For example, when a sample of a solid metal is subjected to a high-power laser beam, the large negative pressures created in the metal lead to its instant melting followed by micro-spalling [9]. In an underwater explosion, the detonation of an unconfined charge leads to the growth of a stable gas bubble [10]. In the design of liquid blast mitigants, the detonation of a confined explosive surrounded by a liquid layer leads to the jetting of finger-like structures or instabilities [11]. In these examples, the compressible effects, the liquid-gas interface motion and surface tension are key fluid physics phenomena [12,13].
Over the past three decades, efficient schemes for multi-phase compressible flow have been put forward. These may be broadly divided into three families of multi-phase modelling. First, the "two-phase" models where each phase is treated explicitly by its own set of equations. In this regard, the Baer and Nunziato [14] and Abgrall and Saurel [15] seven-equation models are popular [16][17][18]. The second approach reduces the two energy equations to a single one, which is equally valid. Examples include the reduced fiveequation model of Kapila et al. [19] used in [1,[20][21][22]. The third and final approach is a homogeneous flow model known as the "one-fluid" formulation, where an equilibrium is assumed to exist between the liquid and the gas phases. This one-fluid formulation has been extensively used in incompressible flow [23][24][25][26], and recently in compressible shock

Governing Equations
Inviscid immiscible compressible two-phase flow is considered where the governing equations are based on that by Abgrall and Saurel [15]. Here, a one-fluid formulation is derived from their work by defining a cell-average density ρ, pressure p and a uniform velocity u. In the absence of gravitational, viscous effects and heat conduction or phase change, these read in compact form as where U = ρ, ρu, ρE, α T is the set of conserved variables, F and S are the flux and source vector written, respectively, as In the above, ρ, u and p are, respectively, the fluid density, velocity and pressure. Further, I is the identity matrix, and the volume fraction, α, is computed as a cell average of the volume occupied by the tracked phase i to the volume of the cell (α = V i V l ). The volume fraction transport equation is expressed in its semi-conservative form as per Johnsen et al. [44] and by Weymouth et al. [36]: where the right-hand-side is referred to as the dilatation term (compression and expansion of the tracked phase). Moreover, surface tension is included as a volumetric force as per the continuum surface force (CSF) method [4] as with δ s the Dirac surface function, n the interface-normal, σ the surface tension coefficient, and κ the interface curvature. The density ρ is volumetrically weighted by with subscripts denoting the liquid and gas, respectively. The total energy ρE is expressed as the sum of the kinetic and internal energy ρe: where ρe is given by ρe = αρ 1 e 1 + (1 − α)ρ 2 e 2 , and the internal energy is computed using the stiffened gas equation of state (EOS).
Finally, normal to the interface, this system is hyperbolic where there are three unique and real eigenvalues (see Appendix B), λ 1 = u · n − c, λ 2 = u · n and λ 3 = u · n + c with c denoting the acoustic velocity. The next subsection details the closure (EOS) for the above system.

Equation of State (EOS)
The governing equations are closed by the caloric relationship p = p(ρ, e) established via the stiffened gas EOS. This EOS has been extensively used in multi-phase compressible flow [12,44,45] and is written in Mie-Gruniessen form as follows: where γ i denotes the adiabatic expansion coefficient and p ∞ i is the empirical pressure for the phase i. This pressure arises due to the intermolecular forces such as van der Waals forces, which occur in a real gas or liquid and reduce to the ideal gas law when p ∞ = 0. The speed of sound c i for each phase is defined as where the nomenclature has previously been defined.

Mixture Rules
We note, from previous work [12,20,44,45], that a particular form of the material property, ( 1 γ−1 , γp ∞ γ−1 ) must be employed to guarantee an oscillation-free velocity and pressure field. This follows from an analysis of the internal energy equation, notably, when considering the propagation of a free surface in a steady uniform velocity and constant pressure field. As per Shyue [45], the material properties are volumetrically weighted as where the above expressions ensure consistency between the VoF and energy equation. Finally, the mixture acoustic velocity c is expressed as with the nomenclature as previously defined.

Thermal EOS
In most articles [18,20,45], the caloric form of the stiffened gas EOS is provided p = p(e, ρ). However, little attention is given to its thermal formulation i.e., T = T(e, ρ). Here, we follow the procedure outlined by physicists in [46] for the Noble-Abel stiffened gas EOS to obtain the following expression for the temperature T as (see Appendix A for derivation) where c v is the specific heat at constant volume. In the work by Metayer et al. [46], it was noted that allowing for A to be non-zero leads to isothermal curves that are non-monotonic where the acoustic velocity may become imaginary. To ensure that physical temperatures are computed, it is recommended in [46] to impose A = 0. Hence, Finally, for a two-phase mixture, T is computed using where The next section details the numerical procedure used to solve the set of conservation equations.

Finite Volume Dual-Cell Mesh
In this work, the chosen discretisation method is the finite volume vertex-centred median dual-cell mesh variant since it is applicable to general meshes. As illustrated in Figure 1, a computational median dual cell, Ω l , is constructed around node l by connecting the mid-points and centroids of the adjacent edges and elements, respectively. The cell, of volume V l , is bounded by facets of the set ∂Ω l with normals n l,m i and boundary faces of the set A b with normals n b, f i . l m n l,m,1 In ELEMENTAL ® , in the interest of computational efficiency, the spatial discretisation is done in an edge-wise manner. For Ω l , the governing equations read in semi-discrete form as follows: where for an edge connecting nodes l to m, denoted t l,m , the face coefficient is given by the sum of the product of the area and the normals of the shared face segments that straddle the edge. Mathematically, this reads for an internal edge as An| ∂Ω l,m = ∑ i A ∂Ω l,m i n l,m i and for a boundary face, For the purpose of computing the face-flux values, HLLC is employed. The involved left and right states are reconstructed via a second-order accurate method as detailed in the next section.

Second-Order Spatial Reconstruction
The Monotone Upstream-centred Scheme for Conservation Laws [41] (MUSCL) is implemented in this work to obtain spatial second-order accuracy on smooth fields for the left and the right state. For a variable φ, across the face of the internal edge t l,m , the left φ L and right φ R state are computed via a limiter, which is required due to the presence of sonic shocks as where with ∆φ = φ m − φ l and where t l,m = x l − x m . Here, ψ(r) is the flux limiter with r = ∆φ k ∆φ and k denotes an arbitrary side k ∈ {L, R}. Finally, the van Albada [47] flux limiter is implemented in this work as follows: where the nomenclature has previously been defined.

Temporal Integration and Stability
In this work, the time integration scheme is the explicit second-order Runge-Kutta (mid-point) method defined as with n the current time-step n + 1 2 a Runge-Kutta sub-iteration, n + 1 is the updated time-step and the global stable time step-size, ∆t, is selected as For inviscid compressible flow modelling on isotropic meshes, with N denoting the total number of nodes in the domain, c l is computed as per Equation (12), ∆x l denotes the effective mesh spacing for l, CFL denotes the Courant-Frederichs-Lewy number where CFL ∈ (0, 1) for a stable scheme. Finally, for surface tension effects, stability requires as per Brackbill et al. [4].

HLLC Solver
In this work, the HLLC solver [30] is employed to obtain the inviscid edge flux. This solver is preferred here for two reasons. First, it ensures positivity of variables and secondly, a volumetric computed flux based on characteristic wave-speeds is readily available.
As illustrated in Figure Rarefaction Equation (22) gives the generic form for all existing HLLC solvers in the literature. This defines the starting point for the inclusion of CICSAM and surface tension effects in the solver. In particular, the improvements that will be proposed for the computation of the intermediate flux will draw from existing work on HLLC [1] and will be further expanded on in the next sections.

Consistency Conditions
To obtain the intermediate star-state flux, consistency conditions on the velocity and pressure fields must first be set. For the face normal and tangential components of the velocity field at the contact discontinuity, the following must hold: u * ,L · n = u * ,R · n = s * and u * k − (u * k · n)n = u k − (u k · n)n. (23) Traditionally, in the absence of surface tension, the following pressure condition is enforced: However, to account for the pressure jump induced by surface tension, a new consistency condition is derived. Using the generalised Riemann invariants and eigenvector k (4) = 0, 0, 0, σκ, 1 T , the Young-Laplace pressure jump condition is obtained: where this consistency pressure condition was first proposed in [1].

Numerical Energy Consistency Criteria
The incorporation of CICSAM into HLLC begins by considering a key difference between the manner in which MUSCL and CICSAM operate. Notably, MUSCL is a monotone up-winding scheme that interpolates left and right states on either side of a face. These are then blended by HLLC based on characteristic wave-speeds to arrive at a discretised face-flux. In contrast, CICSAM employs a blend of up-wind differencing and down-winding to reconstruct a flux based on availability (amount of fluid available in the cell). As a result, CICSAM yields a more accurate discretisation of the α flux in the VoF equation. As HLLC is more desirable where compressible flow characteristics are key (ρ, u, p etc.), a consistent blend with CICSAM is sought. To enable this, face-fluxes are discretised in terms of primitive variables (ρ, u, p) in this work. This is illustrated.
Consider the special case of an interface problem only involving constant pressure and velocity. Here, the momentum and energy equations reduce in 1-D to The VoF equation is ∂α ∂t Since the sharp propagation of the interface is desired, the semi-discrete form of the VoF equation reads where α f ,C denotes the face value associated with CICSAM and u f f is the face-flux given by Equation (33). If the energy field is reconstructed in terms of conserved variables using the HLLC-MUSCL approach, the thermal equation is discretised as where subscript H denotes the internal energy flux obtained from HLLC associated with the MUSCL reconstruction. Introducing the gas stiffened EOS, the above expression then reads ∂ ∂t The above equation can be split into two transport sub-equations as ∂ ∂t Since the pressure is constant and non-zero, Equation (27a) reduces to which can be further expanded in terms of the volume fraction field as ∂ ∂t where α f ,H represents the face volume fraction field that is consistent with the HLLC-MUSCL reconstruction of the energy field. Equation (28) can be summarised for a phase i ∈ {1, 2} as ∂ ∂t where the same conclusion can be drawn from Equation (27b). Further simplification of the above expression leads to where α i ∈ {α, 1 − α} for a two-phase flow. It is clear that Equation (29) is inconsistent with Equation (26) as α f ,C = α f ,H . This will lead to the appearance of spurious currents in the pressure as well as the velocity field. The above demonstrates that a consistent discretisation method must be employed on both the energy and VoF flux in order to maintain an oscillation-free field. This relationship is set up by the caloric gas stiffened EOS and will be referred to as the energy consistency criteria (ECC). Hence, with the introduction of the VoF method, the ECC can only be achieved if the energy field is discretised via the EOS by applying a blend of HLLC and CICSAM as follows: Expanding the above as per Equations (27) and (28) will lead to the desired VoF discretisation, i.e., Equation (26) is recovered. This ensures consistency between the VoF CICSAM flux and the energy term computed by HLLC. Furthermore, since at the face of an edge, sharp interface properties are sought, it is proposed in this work that the EOS properties are discretised using CICSAM as follows: This ensures an accurate discretisation of the energy field as will be shown in Section 4.

Intermediate Star-State
In this sub-section, the following notation will be employed to differentiate between different scheme variants in discretising a variable φ at a face: φ k,M will refer to MUSCL, φ f ,C to CICSAM and φ k on its own will define a combination of the two aforementioned operators for an arbitrary side k ∈ {L, R}. The HLLC normal flux may be re-written from Equation (22) as where the flux for an arbitrary side k is computed via and the intermediate star-state is derived as shown in Appendix C as The sonic characteristic, χ k , is defined as where s * denotes the contact wave computed as per [1,48], The first term results from the 'traditional' HLLC convective wave-speeds. The second term is the result of capillary effects over the two-phase interface where α k,σ denotes the face value of alpha that is consistent with the discretisation of the surface tension source term. This is further expanded on in Section 3.5.4.
The energy field is discretised as The wave-speeds are computed as where s L and s R denote the characteristic left and right wave-speed, respectively. As per Einfeldt et al. [49], the wave-speeds are estimated using the Roe average eigenvalues as Further, H k is the specific stagnation enthalpy computed via Finally,φ ∈ {ũ,H} is computed using where the nomenclature has previously been defined.

Face-Flux
Given the equation for the intermediate flux (Equation (31)), a mathematical equivalent expression for the face-flux, u f f , which is consistent with the proposed CICSAM discretisation, can be derived. Using an adaptation to the HLLC solver proposed by Johnsen et al. [44], it is assumed that αu f f = α f ,C u f f . The expression for the alpha flux hence follows as If s * > 0, it follows that , and similarly, for s * < 0, , Therefore, for the edge t l,m , the face-flux is In particular, the above definition allows for easy implementation of the VoF and HLLC solver in a manner that guarantees the contact-preserving properties.
Finally, though the VoF CICSAM interface reconstruction is used in this work, for the purpose of comparison, the alpha face value, α f ,H , is also discretised using HLLC-MUSCL where the different states are given as and the alpha flux is then computed as per Equation (31).

Surface Tension Term
HLLC solvers rarely take into account surface tension effects for multi-phase compressible flow. The first mention of an HLLC-based scheme with surface tension was that of Garrick et al. [1], where a compression technique was used to remove any numerical diffusion present in the VoF field. In contrast, in this work, the conservative VoF-CICSAM method is combined with the HLLC solver to allow for curvature to be computed accurately. Similar to [1], this requires a well-balanced implementation with the surface tension source term. A well-balanced implementation implies that the surface tension source term is balanced by the computed pressure gradient in a steady flow field with zero velocity and zero gravity. This may be mathematically expressed by considering a static bubble configuration [4,12,50,51] where the Euler equation in volumetric form with a CSF surface tension discretisation reduces to The above is discretised for a node l as Here, ε p and ε σ are the discretisation errors associated with pressure and alpha, respectively.
For the case where s + and s − are non-zero, the pressure face as per Equation (31) is given as with s L = min(−c, −c L ) and s R = max(c, c R ). The contact wave-speed simplifies to In the case of a static bubble, as the pressure field converges to the analytical solution, the consistency pressure condition (Equation (24)) is satisfied and hence the contact wavespeed tends to zero (or machine precision, ǫ). Equation (36) may therefore be expressed as where s * → ǫ.
As explained by Popinet et al. [51], to recover the discrete equilibrium of Equation (35), it is required for ε p to cancel ε σ to zero. To achieve this, the same operator used for the pressure face must also be employed for the volume fraction field in the surface tension term. Hence, α f ,σ is interpolated for the left or right state using MUSCL and up-winded using HLLC as follows: Similarly, the contact wave-speed now reads for the general case, where κ f is averaged as per [1], i.e., κ f = κ l +κ m 2 . This approach for the computation of κ f is also widely used for solving the pressure Poiseuille equation in semi-implicit projection solvers [50,51]. The approximation is valid for a sufficiently fine mesh and for cases where the curvature is uniform and constant.
Finally, the energy surface tension term is re-expressed in its conservative form and discretised as follows: This guarantees a well-balanced method as demonstrated numerically via the static bubble test-case in Section 4.2.4.

Curvature Reconstruction
In this work, the curvature is re-constructed using two methods. Firstly, using the convolution method as proposed by Brackbill et al. [4], where the divergence of the free surface normal of a filtered (smoothed) alpha field α * [25] is used: The above is in the interest of applicability to arbitrary meshes. Secondly, using height functions [2,3], which are available in ELEMENTAL ® and are used to obtain second-order spatial accuracy on Cartesian meshes.

VoF Equation
Finally, the VoF equation is discretised using CICSAM as follows: where α f ,C is computed as per [34] and the blending is done using a smooth filtered field as per [24,25].

Summary of Algorithm
On the first time step, α L/R is reconstructed using Godunov while the primitive variables are reconstructed using MUSCL. Using this initial guess for α L/R , HLLC computes a face-flux u f f using Equation (33). This face-flux is used in CICSAM to recompute the face value for the volume fraction field. The curvature is then computed. The reconstructed variables are inputs to the HLLC solver where the conservative F f and face-flux are computed and stored. It is important to note that in the computation of s * , CICSAM is used for the computation of { 1 γ−1 , γp ∞ γ−1 } for both the left and right state. Hence the convective flux term is now consistent with the VoF flux (the ECC are satisfied). However, for the surface tension term appearing in the contact wave-speed (Equation (32)), MUSCL is used for the left and the right state for alpha. This ensures consistency between the pressure gradient operator and the surface tension source term (the Young-Laplace condition is satisfied). The benefits of this algorithm is that the use of CICSAM ensures that the volume fraction field is always sharp. This allows curvature to be computed accurately. The HLLC-CICSAM algorithm is summarised in Algorithm 1.

Numerical Test Cases
The novelty in this article affects all terms in the governing equation. As such, a range of test-cases were conducted increasing in complexity from a physics perspective. Unless otherwise specified, in the interest of VoF accuracy, a CFL = 0.1 is used.

1-D Problems
The single and two-phase compressible attributes of the proposed CICSAM-HLLC solver are evaluated in this sub-section using 1-D test-cases. For all described cases, the domain is an equispaced structured mesh where outflow conditions are prescribed at each end.

Gas-Liquid Riemann Problem
The first test-case considered is that of a 1-D two-phase Riemann problem [45]. The objective is to assess the accuracy due to CICSAM as compared to HLLC in discretising the VoF equation. The test-case consists of a gas-liquid unit domain where the following initial conditions are applied:   To quantify this improved accuracy, a mesh independence study was conducted using the L 1 norm, ε 1 , of the absolute error, where φ(x l , t) denotes the value for a variable at position x l and time t. Further, φ nu and φ an denote the numerical and analytical solutions, respectively.
From Figure 4 right, for the volume fraction field, it is seen that CICSAM achieves an L 1 norm of at least an order of magnitude lower than that of HLLC. Similarly, the sharper interface leads to an improvement in the accuracy of the internal energy Figure 4 left. Localised second-order spatial accuracy is also observed with the proposed CICSAM method when refining from a 400 to a 1600 node mesh.

2-D Problems
In this sub-section, the ability of the scheme to handle multi-dimensional flow is demonstrated. Six test-cases are presented. The first is concerned with evaluating the accuracy of the VoF scheme while the second and third involve two-phase compressible flow problems (without surface tension effects). The final three test-cases are concerned with surface tension modelling accuracy. For all test-cases, VoF initialisation is done using the Arbitrary Grid Initialiser (AGI) by Jones et al. [52]. Unless otherwise stated, slip conditions are prescribed at boundaries.

Advecting Bubble in an Oblique Velocity Field
In this test-case, two key components of the interface capturing scheme are evaluated. First, the accuracy of the VoF method. Secondly, the compatibility of CICSAM with HLLC, i.e., its ability to maintain pressure and velocity equilibrium for uniform flow.
Here, a circular bubble of radius R 0 = 0.16 is initialised at the bottom left corner, x = (0.25, 0.25), of a unit domain. The bubble is advected in an oblique uniform and steady velocity field, u = 1.0ê x + 1.0ê y . The initial pressure is set to 1. The material properties are given as per Shyue et al. A mesh convergence study is conducted by repeating the calculation of a set of increasingly finer meshes from 64 2 to 512 2 . The VoF accuracy is assessed using the shape error. The L 1 norm of the shape error, ε g , is computed via (44) Figure 5 demonstrates a better qualitative description of the volume fraction and the energy field. Further, as seen in Figure 6, the shape error with CICSAM is at least an order of magnitude below that of HLLC. CICSAM is clearly superior for advecting the liquid-gas interface.

Underwater Explosion
Here, we implement the popular underwater explosion test-case as described in [20,27]  A rectangular structured mesh with 400 × 250 nodes is used for this simulation. The total simulated time is 2.5 ms. The VoF-HLLC (denoted HLLC) is again compared to the proposed scheme (denoted HLLC-CICSAM). Figure 7 depicts the Schlieren-type images of the interface (density field) at t = 0.4 ms, 0.8 ms, 1.2 ms and 2.5 ms for the proposed HLLC-CICSAM method. Overall, we note that our results show flow features that are qualitatively in agreement with those obtained numerically by (Figure 10) in Shyue et al. [27]. In particular, it is seen in Figure 8 that the HLLC-CICSAM formulation allows for a qualitative improvement on the energy field (less diffuse). Further, to show comparative computed pressures, the cross-section along the line x = 0 is plotted as shown in Figure 9. Quantitatively, an overall satisfactory agreement is seen as compared to the interface capturing method proposed by Shyue et al. [27]. A higher pressure of up to a factor of 2 is seen in certain regions with HLLC-CICSAM as compared to HLLC. These differences are attributed to a sharper interface.

Shock-Bubble Interaction
Next, we consider the well-known shock-bubble interaction (see Figure 10), where a planar Mach 1.22 shock wave moving from right to left collides with an R22 bubble in air. This test-case effectively deals with the two species, R22 and air, modelled as ideal gases (p ∞ = 0). Different variations of this test-case exist in the literature [20,27,44].  Here, we implement the test-case as described by Shyue et al. [27]. An R22 gas bubble with a radius of R 0 = 25 mm is centred at x c = (225, 44 The simulation is run on a 3560 × 356 node mesh at a CFL of 0.9 for a total time of 1020 µs. Figures 11 and 12 illustrate the evolution of total energy as a function of time and again HLLC is compared with CICSAM for the VoF equation. For the same adiabatic coefficient, the speed of sound inside the R22 gas bubble is slower compared to the air around. Hence, the wave generated inside the bubble is a rarefaction wave while the incoming wave is a shock wave. The interaction of these two waves produces two outgoing transmitted waves that get reflected at the boundary and free surfaces. These subsequent waves result in the denser fluid (R22) being decelerated by the lighter fluid (air). These lead to the appearance of Rayleigh-Taylor instabilities as illustrated at t = 690 − 1020 µs. For the purpose of comparison with experimental data, the Schlieren-type images are shown in Figure 13. Total Energy, ρE, kJ Qualitatively, the plots show flow features that are consistent with numerical results illustrated in (Figure 5) in [27] with similar improvements as previously discussed with respect to the VoF and energy field seen when using CICSAM.

Spurious Currents in a Static Bubble
A classic test-case [1,34,51,53] to validate surface tension schemes is that of a static bubble. It is employed here to validate the accuracy of both HLLC and CICSAM VoF methods in conjunction with height function and convolution curvature calculation. An ideal gas bubble with radius R 0 = 0.4 is initialised in a slightly compressible liquid at the centre of a unit domain. The material properties are set as per Fuster et al. [12]: (ρ 1 , γ 1 , p ∞ 1 ) = (1.0, 1.4, 0.0)and(ρ 2 , γ 2 , p ∞ 2 ) = (1.0, 7.14, 300), with σ = 1 and where the grid is such that R 0 ∆x = 12.4 (equispaced structured). A uniform pressure field is initialised. The simulation is run for 15 s (≈179,229 steps at ∆t ≈ 8e −5 ) to allow the system to reach steady state. The infinity norm of the maximum pressure jump, ε ∆p ∞ , is computed using where ∆p is the difference between the maximum and minimum pressure in the domain, with p nu and p an the numerical and analytical pressure, respectively. Further, the velocity is non-dimensionalised with respect to the characteristic inviscid velocity U σ , which is computed as per Popinet et al. [50]: The quantitative results are summarized in Table 1 showing the clear superiority of CICSAM with height functions if used on a Cartesian mesh. Figure 14 illustrates the evolution of the velocity currents over time. The maximum velocity recorded with CICSAM and height functions is of the order of 10 −8 while that with HLLC is divergent. This shows that the proposed scheme is well-balanced as well as the importance of accurate curvature calculation. The errors when using convolution to compute a pressure jump are significantly higher due to a larger error in curvature similar to [54].

Oscillating Bubble
A more challenging test-case [1,29,53,55,56] to validate the implementation of surface tension effects in the solver is to simulate the periodic deformation of a liquid droplet in air in the absence of gravity. The equation for the initially deformed droplet, as per Torres et al. [55], is given in polar coordinates as where n is the integer mode of oscillation and R 0 is the initial radius. Here for n = 2, the above equation is written in Cartesian coordinates as with ǫ r = 0.01 and R 0 = 0.8.  where the densities of the gas and liquid are ρ 1 = 0.01 and ρ 2 = 1 with σ = 1. The theoretical small-amplitude inviscid oscillation frequency is then The domain is a square of size [−2, 2] m and the velocity field is initialised to zero everywhere. The simulation is run for a total time of two periods, using a CFL of 0.4. Figure 15 illustrates the evolution of the overall kinetic energy where the damping is due to the presence of numerical dissipation at low mesh resolution as seen in the work by Garrick et al. [1]. Here, the maximum % error, as compared to the analytical period on the coarsest mesh, is 1.1%, while on the finest, it is 0.15%. The results are comparable to those obtained by others [29,55,56]. As shown in Figure 16, second-order convergence is recovered while refining from a 200 2 to 400 2 mesh.

Rayleigh-Plesset Collapse Problem
A popular test-case [12,40,57] to evaluate the effect of non-linear terms in multiphase compressible flow is the so called Rayleigh-Plesset bubble expansion-contraction problem [58]. Here, we consider a cylindrical bubble (2-D) in an incompressible liquid. A gas bubble of radius R 0 at initial pressure p b,0 expands and contracts in an incompressible liquid due to a sudden change in the surrounding pressure. The evolution of the radius is governed by the Rayleigh-Plesset model [59]. The initial pressure field is where R,Ṙ,R are the bubble radius, interface velocity and acceleration of the bubble, respectively. Further, p S,0 is the initial pressure at a finite distance S. For the purpose of simulating the problem, a bubble of radius R 0 = 2 µm is first initialised at the centre of a square domain of size [−5, 5] µm and a Cartesian mesh is employed. The driving pressure function at the outer boundary is expressed as where ω c = 10208967.75 s −1 as per [59,60]. The material properties are as per Fuster et al. [12]: with σ = 0.1. At each iteration, a Dirichlet pressure boundary condition (pulse) is set at the outer boundary as per (46). The simulation is run at a CFL number of 0.5 for 0.5 s. Figure 17 illustrates a time lag on the numerical solution with an overshoot on rebound, which is invariant with respect to mesh refinement. This "incorrect" solution is due to the use of a finite Cartesian domain where such effects were first seen in [40]. To demonstrate this and also to show the applicability of the scheme developed in this paper to nonorthogonal grids, the calculation is repeated using a curvilinear mesh of radius 5 µm where elements are grown in size outside of the bubble. Since height functions are limited to structured grids, the convolution method is used for the curvature calculation.  Figure 18 illustrates that the numerical solution now asymptotes to the analytical, correcting the time lag and overshoot in Figure 17. In Table 2, the relative error on the minimum radius of collapse, the time period and the amplitude of rebound are recorded. Here, a maximum error of 4.9% is noted on the minimum radius on the coarsest mesh and 1.2% on the finest mesh. Satisfactory convergence with the analytical solution is obtained on the finest mesh ( R ∆x = 80). Figure 19 illustrates the evolution of the bubble at the start, mid-way and end time of the simulation. As shown, the CICSAM algorithm allows for sharp capturing of the interface retaining the geometric integrity of the bubble. This test-case demonstrates that the scheme can accurately predict the non-linear collapse of a poly-tropic gas bubble.

Conclusions
In this article, we have presented an all-Mach HLLC-based scheme capable of accurately modelling a homogeneous compressible two-phase flow in the presence of surface tension. The solver for the first time combines HLLC (with primitive variable reconstruction) and the conservative VoF CICSAM method for the sharp propagation of the interface. Here, the reconstruction of primitive variables allowed for the consistent integration of the VoF and energy flux where the ECC was satisfied. The addition of CICSAM to HLLC allowed for the accurate calculation of interface curvature and therefore accurate capture of surface tension effects.
Maxwell's rule relates the internal energy to the temperature, pressure and specific volume as follows: By substituting Equation (A3) into Equation (A4), the following expression is obtained: Further, substituting for the specific heat capacity at constant volume, Integrating with respect to T yields Next, differentiating the above with respect to the volume at constant temperature leads to Substituting Equations (A6) and (A7) into (A5) leads to the following expression: Substituting Equation (A9) into Equation (A8) yields Using an integrating factor such as e γ v dv , an expression for β(v) is obtained as follows: where A denotes the constant of integration. Then, substituting Equation (A11) into Equation (A7), an expression for p(v, T) is derived as Using a reference state (p 0 , T 0 , v 0 ), an expression for A is obtained as Re-expressing the above equation in terms of the reference density ρ 0 yields Finally, substituting Equation (A14) into Equation (A7) and re-expressing in terms of density, the equation for temperature reads where the nomenclature has previously been defined.

Appendix B. Eigenvalues and Vectors
The governing equation set presented in this work may expressed in strong compact form in 2-D as ∂U ∂t where subscripts x and y denote the flux in the respective Cartesian directions. Further,  The characteristic eigenvector corresponding to the Jacobian B matrix associated with the remaining dimension can easily be derived in a similar way and is not included here. Finally, for a system of equations of arbitrary dimension d, normal to the interface, there are three unique and real eigenvalues: λ 1 = u · n − c, λ 2 = u · n and λ 3 = u · n + c (multiplicity of spatial dimensions + one for each transport equation), where the interface normal is n.

Appendix C. Derivation of Star-State
Consider the general formulation for the HLLC solver: By substituting for known quantities U k and F k , an expression for the star-state region is obtained. First, for the continuity equation, this reads ρ * k u * k · n = ρ k u k · n + s k ρ * k − ρ k , ρ * k s * = ρ k u k · n + s k ρ * k − ρ k . Hence, and for the momentum equation, this is expanded as follows: ρ * k u * k u * k · n + p * k n = ρ k u k u k · n + p k n + s k ρ * k u * k − ρ k u k , p * k n = ρ k u k u k · n + p k n + (s k − s * )ρ * k u * k − s k ρ k u k , From consistency Equation (23), the following expression must hold: where u * n,k and u t are the normal and tangential velocities, respectively. Hence, u * k = u * n,k + u t , u * k = s * n + u k − (u k · n)n, u * k = u k + (s * − u k · n)n.
Substituting Equations (A18) and (A20) into Equation (A19), an expression for the intermediate pressure, p * k , is derived as follows, p * k = p k + ρ k (s k − u k · n)(s * − u k · n) (A21) Substituting Equation (A21) into Equation (A19), an expression for the intermediate momentum is obtained: Similarly, for the energy equation, the HLLC flux is written as Substituting the consistency Equations (23) and (A21) into Equation (A23) yields the expression for the intermediate star-state energy term as Finally, using consistency Equation (24), the following expression for the contact wave speed is obtained: where s k is computed as per Einfeldt et al. [49].