A unified algorithm for interfacial flows with incompressible and compressible fluids

Density-based algorithms, where the sought primary variables are the density, momentum and energy of the flow, are widely applied for compressible flows. With respect to the modelling of compressible interfacial flows, notable examples are the two-fluid Baer-Nunziato (BN) model [3], where each fluid is represented by its own set of governing equations, and the five-equation models [4, 5], with a separate continuity equation for each phase, and shared momentum and energy equations. In these algorithms, an exact or approximate Riemann solver is employed to evaluate the fluxes at cell-faces of the computational mesh [6, 7, 8]. The ghost-fluid method (GFM) [9] has established itself as an alternative to solving a Riemann problem, partly due to its conceptual simplicity. Recent developments have also seen a combination of GFM with Riemann solvers [10, 11] to improve the stability of simulations with strong shock-interface interactions and compressible gas-liquid flows.

Since pressure is not directly solved for with a density-based algorithm, the pressure field has to be reconstructed based on the applied thermodynamic model via an equation of state, which may be difficult in interfacial cells where two bulk phases coexist [4, 5, 12], and which requires hydrodynamically and thermodynamically plausible fluid properties. While density-based algorithms are naturally suited for compressible flows, they are poorly suited for low-Mach number flows [13, 14], where the coupling of pressure and density vanishes. Although density-based algorithms have been applied to low-Mach number flows with some success, this requires computationally expensive pre-conditioning techniques [15].

A different approach for single-phase flows at all speeds was pursued by van der Heul et al. [16], in which the energy equation is reformulated as an equation for pressure, while the continuity equation still serves as an equation for density. Similar methods were subsequently presented by Park and Munz [17] and Cordier et al. [18]. These algorithms formally converge to the incompressible flow equations in the limit of zero Mach number ($M\rightarrow 0$) [16, 17], with the reformulated energy equation enforcing a divergence-free velocity field. Boger et al. [19] extended the work of Park and Munz [17] to interfacial flows.

Pressure-based algorithms, in which the continuity equation is cast into a discretised equation for pressure (with pressure acting as the Lagrange multiplier that enforces $\mathbf{\nabla}\cdot\mathbf{u}\rightarrow 0$ for $M\rightarrow 0$), while density is either constant (incompressible fluids) or evaluated using a suitable equation of state (compressible fluids), are preferably applied to simulate the flow of incompressible fluids and weakly compressible flows, and may yield significant advantages for low-Mach number compressible flows, since the acoustic degeneration at low Mach numbers does not pose a problem [20]. The success of pressure-based algorithms is facilitated by the unique role of pressure in all Mach number regimes [20], with the pressure-velocity coupling dominant at low Mach numbers and the pressure-density coupling dominant at high Mach numbers [21, 22], and the convenient fact that the fully-conservative formulation of the governing conservation laws can still be satisfied accurately [1, 21].

The majority of pressure-based algorithms for incompressible interfacial flows are founded on pressure-correction methods, such as projection methods [23, 24], the SIMPLE method [25] and its subsequent derivatives, or the PISO method [26], or fully-coupled algorithms [27]. Pressure-based algorithms for compressible flows are much less prominent in the literature than their density-based counterparts, partly because deriving stable and efficient numerical schemes for the transonic regime [28] as well as formulating consistent shock-capturing schemes have proven difficult for pressure-based algorithms [14, 29]. Nonetheless, because pressure plays an important role in all Mach number regimes, pressure-based algorithms have a distinct potential for applications in all Mach number regimes or in which the Mach number varies strongly. Starting with the work of Harlow and Amsden [30], a variety of pressure-based algorithms for compressible single-phase flows has been developed [13, 21, 28, 31, 32]. Denner et al. [1] presented a fully-coupled pressure-based algorithm able to predict single-phase flows of incompressible, ideal-gas and real-gas fluids at all speeds ($0\leq M\leq 239$) accurately and robustly, using the same fully-coupled pressure-based algorithm and conservative finite-volume discretisation.

Despite these developments for single-phase flows, it was only recently that the first pressure-based algorithm in conservative form for compressible interfacial flows at all speeds was proposed [2]. This algorithm was proposed in conjunction with a new interface discretisation method that retains the acoustic features of the flow, without the need to employ a Riemann solver, and has been shown to be a promising alternative to traditional density-based algorithms, with the pressure-based formulation of the governing equations facilitating the definition of consistent mixture rules at the interface that apply naturally to flows in all Mach number regimes. This algorithm was subsequently extended to polytropic interfacial flows [33], i.e. treating the flow as a polytropic process.

Caiden et al. [34] were the first to propose a numerical algorithm dedicated to the simulation of general compressible-incompressible flows, solving different governing equations for the compressible fluid and the incompressible fluid, coupling the bulk phases at the interface using the GFM. A similar method to describe the behaviour of compressible bubbles in an incompressible fluid was subsequently presented by Aanjaneya et al. [35]. Wadhwa et al. [36] proposed a method for computing incompressible liquid drops in a compressible gas, representing the interface as a matching moving mesh and solving the flow of the incompressible fluid using the artificial compressibility method of Chorin [23]. Others proposed algorithms for compressible-incompressible flows applying the same governing equations in non-conservative form in both fluids [37, 38, 39] and these algorithms are, in general, only applicable to incompressible and low Mach number flows. However, the application of the governing equations in conservative form together with a conservative discretisation is a prerequisite for an accurate prediction of flows in all Mach number regimes [40], in particular shock waves, rarefaction fans and contact discontinuities. The density-based algorithm of Boger et al. [19] is readily applicable to interfacial flows in all Mach number regimes, including the incompressible limit ($M\rightarrow 0$), yet incompressible fluids ($M=0$) have not been considered in the context of such an algorithm.

In addition to the difficulties encountered when developing numerical methods for interfacial flows, such as a precise discrete force balance between surface tension and the pressure gradient, Caiden et al. [34] and Billaud et al. [37] identified the transmission of waves at the interface, as well as physically realistic and compatible discrete formulations of the governing equations describing the incompressible fluid and the compressible fluid, as the main additional difficulties associated with an accurate prediction of compressible-incompressible interfacial flows. The presented algorithm addresses and resolves these difficulties, as demonstrated by the chosen validation cases.

Considering, for example, the advective-diffusive transport of a general fluid variable, $\phi$, is given as

$$\frac{\partial(\rho\phi)}{\partial t}+\mathbf{\nabla}\cdot(\rho\mathbf{u}\phi)=\mathbf{\nabla}\cdot(\Gamma_{\phi}\,\mathbf{\nabla}\phi),$$

(16)

where $t$ is the time, $\rho$ is the fluid density, $\mathbf{u}$ is the fluid velocity and $\Gamma_{\phi}$ is the diffusion coefficient of $\phi$. Reformulating Eq. (16) in its integral form with respect to the control volume $V$, given as

$$\iiint_{V}\frac{\partial(\rho\phi)}{\partial t}\,\textup{d}V+\iiint_{V}\mathbf{\nabla}\cdot(\rho\mathbf{u}\phi)\,\textup{d}V=\iiint_{V}\mathbf{\nabla}\cdot(\Gamma_{\phi}\,\mathbf{\nabla}\phi)\,\textup{d}V,$$

(17)

allows a rather straightforward finite-volume discretisation of the advective-diffusive transport of $\phi$. In the interest of simplicity and brevity, the mesh is henceforth assumed to be Cartesian with a local mesh spacing $\Delta x$. Extending the presented methods to unstructured meshes is easily achieved by introducing corrections for mesh skewness and non-orthogonality, as for instance described in [1].

The transient term is discretised in the following using the Second-Order Backward Euler scheme, which is given for cell $P$ as

$$\iiint_{V}\frac{\partial\Phi}{\partial t}\,\textup{d}V\approx\frac{3\Phi_{P}^{(n+1)}-4\Phi_{P}^{(t-\Delta t)}+\Phi_{P}^{(t-2\Delta t)}}{2\Delta t}\,V_{P}+\mathcal{O}(\Delta t^{2}),$$

(18)

where $\Phi=\rho\phi$, $\Delta t$ is the time-step, superscript $(n+1)$ denotes the implicitly sought solution, superscript $(t-\Delta t)$ denotes the previous time-level and superscript $(t-2\Delta t)$ denotes the previous previous time-level.

The discrete form of the advection term is obtained by applying the divergence theorem,

$\displaystyle\iiint_{V}\mathbf{\nabla}\cdot(\rho\mathbf{u}\phi)\,\textup{d}V$

$\displaystyle=\oiint_{\partial V}\rho\mathbf{u}\phi\,\textup{d}\mathbf{\Sigma},$

(19)

where $\mathbf{\Sigma}$ is the outward-pointing surface vector on the surface $\partial V$ of control volume $V$. Since the surface of a discrete control volume is constituted by a finite number of flat faces $f$ of area $A_{f}$, as illustrated in Fig. 1, and applying the midpoint rule [44, 45], the advection term follows as

$\displaystyle\oiint_{\partial V}\rho\mathbf{u}\phi\,\textup{d}\mathbf{\Sigma}$

$\displaystyle\approx\sum_{f}\rho_{f}\vartheta_{f}\phi_{f}A_{f}=\sum_{f}\dot{m}_{f}\phi_{f}.$

(20)

The face value $\phi_{f}$ is typically not readily available and has to be interpolated from the cell-centred values, which in a finite-volume sense represent the values of $\phi$ averaged over the respective cell volume, using an appropriate interpolation scheme. The advecting velocity $\vartheta_{f}=\mathbf{u}_{f}\cdot\mathbf{n}_{f}$, where $\mathbf{n}_{f}$ is the normal vector of face $f$ (pointing outwards with respect to cell $P$, see Fig. 1), represents the velocity normal to face $f$ and forms, together with the face density $\rho_{f}$ and the face area $A_{f}$, the mass flux $\dot{m}_{f}=\rho_{f}\vartheta_{f}A_{f}$ through face $f$. In the context of the presented algorithm, the advecting velocity $\vartheta_{f}$ is determined using a momentum-weighted interpolation from the cell-centred values of velocity and pressure, as further discussed in Section 5.2. Since the values of density and velocity reside at the cell centres, a suitable interpolation must be applied, such as upwind differencing, central differencing or a TVD scheme.

Similar to the advection term, applying the divergence theorem in conjunction with the midpoint rule, the diffusion term of the transport equation (17) is discretised as

$$\iiint_{V}\mathbf{\nabla}\cdot(\Gamma_{\phi}\,\mathbf{\nabla}\phi)\,\textup{d}V\approx\sum_{f}\Gamma_{\phi,f}\left(\mathbf{\nabla}\phi|_{f}\cdot\mathbf{n}_{f}\right)A_{f},$$

(21)

where the diffusion coefficient $\Gamma_{\phi}$ at face $f$ is defined by a harmonic average of the cell-centred values as [46]

$$\Gamma_{\phi,f}=\frac{2}{\Gamma_{\phi,P}^{-1}+\Gamma_{\phi,Q}^{-1}},$$

(22)

with cell $Q$ the neighbour cell of cell $P$, as illustrated in Fig. 1. The face-centred gradient of $\phi$ along the normal vector $\mathbf{n}_{f}$ is approximated with second-order accuracy as

$$\mathbf{\nabla}\phi|_{f}\cdot\mathbf{n}_{f}\approx\frac{\phi_{Q}-\phi_{P}}{\Delta x},$$

(23)

where $\Delta x$ is the mesh spacing.

Applying the divergence theorem, the spatial gradient of $\phi$ averaged over the cell volume $V_{P}$ is readily computed as

$$\mathbf{\nabla}\phi_{P}\approx\frac{1}{V_{P}}\sum_{f}\phi_{f}\,\mathbf{n}_{f}\,A_{f}.$$

(24)

The momentum-weighted interpolation (MWI) is applied to define an advecting velocity $\vartheta=\mathbf{u}_{f}\cdot\mathbf{n}_{f}$ at cell faces, which is used in the discretised advection terms of the governing equations. MWI provides a robust pressure-velocity coupling for flows with low Mach numbers and flows of incompressible fluids by introducing a cell-to-cell pressure coupling and by applying a low-pass filter acting on the third derivative of pressure [47, 48], thus avoiding pressure-velocity decoupling due to the collocated variable arrangement.

Following the unified formulation of the MWI proposed by Bartholomew et al. [48], the definition of the advecting velocity includes modifications to the original formulation of the MWI, as introduced by Rhie and Chow [49], to account for large density ratios and source terms occurring in interfacial flows, and the transient nature of the considered problems. As demonstrated by Bartholomew et al. [48], only the driving pressure gradient $\mathbf{\nabla}P=\mathbf{\nabla}p-\mathbf{S}$, which is the pressure gradient associated with the flow field, should be coupled to the velocity field, whereas source terms and other external contributions should not be coupled to the velocity field. Taking this into account, the advecting velocity $\vartheta_{f}$ at face $f$ is given as

$$\begin{split}\vartheta_{f}=\overline{\mathbf{u}}_{f}\cdot\mathbf{n}_{f}&-\hat{{d}}_{f}\left[\mathbf{\nabla}P_{f}\cdot\mathbf{n}_{f}-\frac{\rho_{f}^{\ast}}{2}\,\left(\frac{\mathbf{\nabla}P_{P}}{\rho_{P}}+\frac{\mathbf{\nabla}P_{Q}}{\rho_{Q}}\right)\cdot\mathbf{n}_{f}\right]\\ &+\hat{{d}}_{f}\,\frac{\rho^{\ast(t-\Delta t)}_{f}}{\Delta t}\left(\vartheta^{(t-\Delta t)}_{f}-\overline{\mathbf{u}}_{f}^{(t-\Delta t)}\cdot\mathbf{n}_{f}\right),\end{split}$$

(25)

where the face velocities $\overline{\mathbf{u}}_{f}$ and $\overline{\mathbf{u}}_{f}^{(t-\Delta t)}$ are obtained by linear interpolation, and

$$\mathbf{\nabla}P_{f}\cdot\mathbf{n}_{f}\approx\frac{P_{Q}-P_{P}}{\Delta x}.$$

(26)

The face density $\rho_{f}^{\ast}$ 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 [48]. The coefficient $\hat{{d}}_{f}$ is defined as [27, 48]

$$\hat{{d}}_{f}=\dfrac{\left(\dfrac{V_{P}}{{e}_{P}}+\dfrac{V_{Q}}{{e}_{Q}}\right)}{2+\dfrac{\rho^{\ast}_{f}}{\Delta t}\,\left(\dfrac{V_{P}}{{e}_{P}}+\dfrac{V_{Q}}{{e}_{Q}}\right)},$$

(27)

where ${e}_{P}$ and ${e}_{Q}$ are the sum of the coefficients of the primary variable $\mathbf{u}$ arising from the advection and shear stress terms of the momentum equations [1].

The MWI formulation given in Eq. (25) is independent of the applied time-step and the error in kinetic energy introduced by the MWI is proportional to $\Delta x^{3}$ [48]. Hence, the convergence of the second-order accurate finite-volume method remains unaffected. The density weighting applied to the cell-centred pressure gradient in Eq. (25) has been shown to yield robust results for flows with large and abrupt changes in density [48], demonstrated for incompressible interfacial flows with a density ratio of up to $10^{24}$ [27, 50]. Including the transient term has previously been shown to be important for a correct temporal evolution of pressure waves [32, 48], which is particularly pertinent for acoustic effects in compressible flows.

Since the discretisation of the governing conservation laws is identical for single-phase and interfacial flows, the following presentation of the discretisation of the governing conservation laws focuses on single-phase flows. The extension of this discretisation to interfacial flows by an appropriate definition of the fluid properties and using the ACID method is described in Section 6.

Applying the numerical schemes described in the previous sections, the discretised continuity equation (1) for cell $P$ is given as

$$\left.\frac{\partial\rho}{\partial t}\right|_{P}V_{P}+\sum_{f}(\rho_{f}\,\vartheta_{f})^{(n+1)}\,A_{f}=0,$$

(28)

where superscript $(n+1)$ denotes implicitly sought solutions. A Newton linearisation [28, 32] is applied to the advection term of the discretised continuity equation, given as

$$(\rho_{f}\,\vartheta_{f})^{(n+1)}\approx\rho^{(n)}_{f}\,\vartheta_{f}^{(n+1)}+\rho^{(n+1)}_{f}\,\vartheta_{f}^{(n)}-\rho^{(n)}_{f}\,\vartheta_{f}^{(n)},$$

(29)

where superscript $(n)$ denotes the latest available solution. The Newton linearisation facilitates a smooth transition from low to high Mach number regions [21, 28] and provides an implicit pressure-velocity coupling for flows at low Mach numbers and flows of incompressible fluids [1, 32].

The discretised momentum equations (2) of cell $P$ are given as

$$\begin{split}&\left.\frac{\partial\rho u_{j}}{\partial t}\right|_{P}V_{P}+\sum_{f}(\rho_{f}\,\vartheta_{f}\,u_{j,f})^{(n+1)}\,A_{f}=-\sum_{f}\overline{p}_{f}^{(n+1)}\,{n}_{j,f}\,A_{f}\\ &+\sum_{f}\mu_{f}^{\ast}\left[\frac{u_{j,Q}^{(n+1)}-u_{j,P}^{(n+1)}}{\Delta x}+\left.\overline{\frac{\partial u_{i}}{\partial x_{j}}}\right|_{f}^{(n)}\,{n}_{i,f}-\frac{2}{3}\left.\overline{\frac{\partial u_{k}}{\partial x_{k}}}\right|_{f}^{(n)}\,{n}_{i,f}\right]A_{f}+S_{j,P}\,V_{P},\end{split}$$

(30)

where $\overline{\square}_{f}$ denotes linear interpolation of the values at the adjacent cell centres and

$$\begin{split}(\rho_{f}\,\vartheta_{f}\,u_{j,f})^{(n+1)}&\approx\rho^{(n)}_{f}\,\vartheta_{f}^{(n)}\,u_{j,f}^{(n+1)}+\rho^{(n)}_{f}\,\vartheta_{f}^{(n+1)}\,u_{j,f}^{(n)}\\ &+\rho^{(n+1)}_{f}\,\vartheta_{f}^{(n)}\,u_{j,f}^{(n)}-2\,\rho^{(n)}_{f}\,\vartheta_{f}^{(n)}\,u_{j,f}^{(n)}.\end{split}$$

(31)

The discretised energy equation (3) of cell $P$ is given as

$$\begin{split}&\left.\frac{\partial\rho h}{\partial t}\right|_{P}V_{P}+\sum_{f}(\rho_{f}\,\vartheta_{f}\,h_{f})^{(n+1)}\,A_{f}=\left.\frac{\partial p}{\partial t}\right|_{P}V_{P}+\sum_{f}k_{f}^{\ast}\frac{T_{Q}^{(n+1)}-T_{P}^{(n+1)}}{\Delta x}\,A_{f}\\ &+\sum_{f}\overline{u}_{i,f}^{(n+1)}\,\mu_{f}^{\ast}\left[\left(\left.\overline{\frac{\partial u_{j}}{\partial x_{i}}}\right|_{f}^{(n)}+\left.\overline{\frac{\partial u_{i}}{\partial x_{j}}}\right|_{f}^{(n)}\right)-\frac{2}{3}\left.\overline{\frac{\partial u_{k}}{\partial x_{k}}}\right|_{f}^{(n)}\right]{n}_{j,f}\,A_{f}\\ &+\,u_{i,P}^{(n+1)}\,S_{i,P}\,V_{P},\end{split}$$

(32)

with

$$\begin{split}(\rho_{f}\,\vartheta_{f}\,h_{f})^{(n+1)}&\approx\rho^{(n)}_{f}\,\vartheta_{f}^{(n)}\,h_{f}^{(n+1)}+\rho^{(n)}_{f}\,\vartheta_{f}^{(n+1)}\,h_{f}^{(n)}\\ &+\rho^{(n+1)}_{f}\,\vartheta_{f}^{(n)}\,h_{f}^{(n)}-2\,\rho^{(n)}_{f}\,\vartheta_{f}^{(n)}\,h_{f}^{(n)},\end{split}$$

(33)

and

$$h^{(n+1)}=c_{p}\,T^{(n+1)}+\frac{\mathbf{u}^{(n),2}}{2}.$$

(34)

The advecting velocity $\vartheta_{f}$ is given by Eq. (25) and is the same in all equations, ensuring a consistent transport of the conserved quantities. Following the work of Ferziger [46], the viscosity $\mu_{f}^{\ast}$ and the thermal conductivity $k_{f}^{\ast}$ at face $f$ are given by a harmonic average of the values at the adjacent cell centres.

Previous studies [2, 28] have demonstrated substantial improvements with respect to the performance and stability of the numerical solution algorithm associated with a Newton linearisation of all governing equations. In particular, such a linearisation enables a smooth transition from low to high Mach number regions [13, 28, 51]. To this end, the advection terms of the governing equations as well as the transient terms of the momentum and energy equations are discretised using a Newton linearisation, see Eqs. (31) and (33), as described in more detail in [28].

For flows in the incompressible limit, with $M\rightarrow 0$, the density is constant along the fluid particle trajectories [52], with

$$\frac{\mathrm{D}\rho}{\mathrm{D}t}=\frac{\partial\rho}{\partial t}+\mathbf{u}\cdot\mathbf{\nabla}\rho=0.$$

(35)

Consequently, the continuity equation is no longer effective as a transport equation for density, but becomes a constraint on the velocity field with $\mathbf{\nabla}\cdot\mathbf{u}\rightarrow 0$ [52], enforced by pressure which acts as a Lagrange multiplier [53]. Inserting Eq. (35) into Eqs. (1)–(3), the governing conservation laws for flows in the incompressible limit reduce to

	$\displaystyle\mathbf{\nabla}\cdot\mathbf{u}$	$\displaystyle=0$		(36)
	$\displaystyle\rho\left(\frac{\partial\mathbf{u}}{\partial t}+\mathbf{\nabla}\cdot(\mathbf{u}\otimes\mathbf{u})\right)$	$\displaystyle=-\mathbf{\nabla}p+\mathbf{\nabla}\cdot\boldsymbol{\tau}+\mathbf{S}$		(37)
	$\displaystyle\rho\left(\frac{\partial h}{\partial t}+\mathbf{\nabla}\cdot\left(\mathbf{u}h\right)\right)$	$\displaystyle=\frac{\partial p}{\partial t}-\mathbf{\nabla}\cdot\mathbf{q}+\mathbf{\nabla}\cdot\left(\boldsymbol{\tau}\cdot\mathbf{u}\right)+\mathbf{S}\cdot\mathbf{u}.$		(38)

Inserting $\rho=\rho_{0}$, see Eq. (10), into the discretised governing conservation laws presented above, the discretised continuity equation (28) reduces to

$$\sum_{f}\vartheta_{f}\,A_{f}=0\ ,$$

(39)

the discretised momentum equations (30) become

$$\begin{split}&\rho_{P}\left(\left.\frac{\partial u_{j}}{\partial t}\right|_{P}V_{P}+\sum_{f}\vartheta_{f}\,u_{j,f}\,A_{f}\right)=-\sum_{f}\overline{p}_{f}\,{n}_{j,f}\,A_{f}\\ &+\sum_{f}\mu_{f}\left[\left(\left.{\frac{\partial u_{j}}{\partial x_{i}}}\right|_{f}+\left.{\frac{\partial u_{i}}{\partial x_{j}}}\right|_{f}\right)\,-\frac{2}{3}\left.\frac{\partial u_{k}}{\partial x_{k}}\right|_{f}\right]{n}_{i,f}\,A_{f}+S_{j,P}\,V_{P},\end{split}$$

(40)

and the discretised energy equation (32) becomes

$$\begin{split}&\rho_{P}\left(\left.\frac{\partial h}{\partial t}\right|_{P}V_{P}+\sum_{f}\vartheta_{f}\,h_{f}\,A_{f}\right)=\left.\frac{\partial p}{\partial t}\right|_{P}V_{P}+\sum_{f}k_{f}\left.\frac{\partial T}{\partial x_{i}}\right|_{f}{n}_{i,f}\,A_{f}\\ &+\sum_{f}\overline{u}_{i,f}\mu_{f}\left[\left(\left.{\frac{\partial u_{j}}{\partial x_{i}}}\right|_{f}+\left.{\frac{\partial u_{i}}{\partial x_{j}}}\right|_{f}\right)-\frac{2}{3}\left.\frac{\partial u_{k}}{\partial x_{k}}\right|_{f}\right]{n}_{j,f}\,A_{f}+\,u_{i,P}\,S_{i,P}\,V_{P}.\end{split}$$

(41)

These are precisely the governing conservation laws of the incompressible limit, Eqs. (36)–(38), discretised with the schemes as discussed above. In addition, Eqs. (39) and (40) are identical to the discretised continuity and momentum equations of the fully-coupled pressure-based algorithm for incompressible interfacial flows of Denner et al. [27].

The discretised governing conservation laws are solved simultaneously in a single linear system of equations, $\mathbf{A}\cdot\boldsymbol{\phi}=\mathbf{b}$. Following Denner et al. [1], the system of equations is solved for the primary solution variables $\chi$, which are the pressure $p$, the velocity vector $\mathbf{u}$ and the temperature $T$. For a three-dimensional computational mesh with $N$ cells, the linear system of governing equations is given as

$$\begin{pmatrix}\mathbf{A}_{p}&\mathbf{A}_{u}&\mathbf{A}_{v}&\mathbf{A}_{w}&\mathbf{0}\\ \mathbf{B}_{p}&\mathbf{B}_{u}&\mathbf{B}_{v}&\mathbf{B}_{w}&\mathbf{0}\\ \mathbf{C}_{p}&\mathbf{C}_{u}&\mathbf{C}_{v}&\mathbf{C}_{w}&\mathbf{0}\\ \mathbf{D}_{p}&\mathbf{D}_{u}&\mathbf{D}_{v}&\mathbf{D}_{w}&\mathbf{0}\\ \mathbf{E}_{p}&\mathbf{E}_{u}&\mathbf{E}_{v}&\mathbf{E}_{w}&\mathbf{E}_{T}\\ \end{pmatrix}\cdot\begin{pmatrix}\boldsymbol{\phi}_{p}\\ \boldsymbol{\phi}_{u}\\ \boldsymbol{\phi}_{v}\\ \boldsymbol{\phi}_{w}\\ \boldsymbol{\phi}_{T}\end{pmatrix}=\mathbf{b}\ ,$$

(42)

where the coefficient submatrix $\mathbf{A}_{\chi}$ of size $N\times N$ holds the coefficients of the primary variable $\chi$ associated with the continuity equation (28), submatrices $\mathbf{B}_{\chi}$, $\mathbf{C}_{\chi}$, $\mathbf{D}_{\chi}$, all of size $N\times N$, hold the coefficients of the primary variable $\chi$ associated with the momentum equations (30), for the $x$-, $y$- and $z$-axes of the Cartesian coordinate system, and submatrix $\mathbf{E}_{\chi}$, also of size $N\times N$, contains the coefficients of the primary variable $\chi$ associated with the energy equation (32). The subvector $\boldsymbol{\phi}_{\chi}$ of length $N$ holds the solution of the primary variable $\chi$. All contributions from previous non-linear iterations or previous time-levels are gathered in the right-hand side vector $\mathbf{b}$ of length $5N$.

The solution procedure performs nonlinear iterations in which the linear system of governing equations (42) is solved using the Block-Jacobi preconditioner and the BiCGSTAB solver of the software library PETSc [54], as described in detail by Denner et al. [2]. No underrelaxation of the governing equations is required.

To advect the fluid interface between two fluids, Eq. (6) is applied to the colour function, Eq. (43), and reformulated as

$$\frac{\partial\psi}{\partial t}+\mathbf{\nabla}\cdot(\mathbf{u}\psi)-\psi\,\mathbf{\nabla}\cdot\mathbf{u}=0.$$

(44)

The advection of the fluid interface is then described by Eq. (44) using an appropriate discretisation method. Two different VOF methods, an algebraic VOF method [2] and a piecewise-linear interface calculation (PLIC) method with Lagrangian advection of the interface [55], are considered for the advection of the fluid interface.

In the algebraic VOF method [2], the advection equation of the colour function, Eq. (44), is discretised using the Crank-Nicolson scheme for the discretisation of the transient term and the CICSAM scheme [56] for the spatial interpolation of the colour function $\psi$. In the applied VOF-PLIC method, the interface is reconstructed based on the local colour function $\psi$ and the normal vector of the interface [57]. The interface advection equation (44) is then rewritten in integral form and the reconstructed interface is advected using the Lagrangian split advection scheme of van Wachem and Schouten [55]. In both considered interface advection methods, the advection is based on the same advecting velocity $\vartheta_{f}$ as for all advection terms of the governing equations, thus ensuring an accurate prediction of volume changes [2].

Nevertheless, the finite-volume discretisation and pressure-based algorithm presented in Section 5 are not limited to the employed VOF methods and other methods to represent the bulk phases and advect the interface, including level-set and front-tracking methods, may equally be applied.

The definitions of the fluid properties largely follow the principles outlined by Denner et al. [2]. The density of the fluid is defined based on the colour function $\psi$ as

$$\rho=(1-\psi)\,\rho_{\mathrm{a}}+\psi\,\rho_{\mathrm{b}},$$

(45)

where the partial densities $\rho_{\mathrm{a}}$ and $\rho_{\mathrm{b}}$ of the bulk phases are given by Eq. (13). This linear interpolation of the density is required in order to satisfy the discrete conservation of mass, momentum and energy. In the context of compressible flows, it is equivalent to an isobaric closure assumption [4, 58], while for incompressible fluids this formulation reduces to the typically used definition of density [59]. The heat capacity ratio also follows from the isobaric closure assumption as

$$\frac{1}{\gamma-1}=\mathcal{C}_{\mathrm{a}}\,\frac{1-\psi}{\gamma_{0,\mathrm{a}}-1}+\mathcal{C}_{\mathrm{b}}\,\frac{\psi}{\gamma_{0,\mathrm{b}}-1},$$

(46)

where $\mathcal{C}_{\mathrm{a}}$ and $\mathcal{C}_{\mathrm{b}}$ are the binary compressibility operators defined in Eq. (12) associated with fluid $\mathrm{a}$ and fluid $\mathrm{b}$, respectively. The specific isobaric heat capacity is defined by a mass-weighted interpolation [2], which is essential for the conservation of the total energy, given as

$$\rho\,c_{p}=(1-\psi)\,\rho_{\mathrm{a}}\,c_{p,\mathrm{a}}+\psi\,\rho_{\mathrm{b}}\,c_{p,\mathrm{b}},$$

(47)

where the partial densities $\rho_{\mathrm{a}}$ and $\rho_{\mathrm{b}}$ are given by Eq. (13), density $\rho$ is given by Eq. (45), and the partial specific isobaric heat capacities $c_{p,\mathrm{a}}$ and $c_{p,\mathrm{b}}$ are given by Eq. (14). The viscosity $\mu$ and thermal conductivity $k$ are defined as

	$\displaystyle\mu$	$\displaystyle=(1-\psi)\,\mu_{\mathrm{a}}+\psi\,\mu_{\mathrm{b}},$		(48)
	$\displaystyle k$	$\displaystyle=(1-\psi)\,k_{\mathrm{a}}+\psi\,k_{\mathrm{b}}.$		(49)

The Continuum Surface Force (CSF) model of Brackbill et al. [59] is applied to model the force due to surface tension, represented by the source term

$$\mathbf{S}=\sigma\,\kappa\,\mathbf{\nabla}\psi,$$

(50)

where $\kappa$ is the interface curvature, which is computed using a second-order height-function method [60], and $\sigma$ is the surface tension coefficient. In the interest of conciseness, but without loss of generality, the surface tension coefficient is assumed to be constant.

The ACID method [2] assumes that, for the purpose of discretising the governing conservation laws for a given cell, all cells in its finite-volume stencil are assigned the same colour function value, i.e. the colour function is assumed to be constant in the entire finite-volume stencil. The relevant thermodynamic properties that are discontinuous at the interface, i.e. density and enthalpy, are then evaluated based on this locally constant colour function field, as described by Denner et al. [2] in detail. This recovers the contact discontinuity associated with the interface [61, 62] and enables the application of the conservative discretisation described in Section 5, identical to the one applied for single-phase flows. Denner et al. [2] reported robust and accurate results for acoustic and shock waves in interfacial flows, supporting the notion that the interface discretisation indeed conserves the acoustic features of the flow and retains the conservative discretisation of the governing equations.

A circular bubble of a compressible gas surrounded by an incompressible fluid is simulated to study the parasitic currents associated with the numerical treatment of the surface tension contribution. Similar to the test-cases used in previous studies [68, 69], the considered two-dimensional bubble has a diameter of $d_{\mathrm{B}}=0.8\,\mathrm{m}$ and a surface tension coefficient of $\sigma=1.0\,\mathrm{N}\,\mathrm{m}^{-1}$. The ambient pressure is $p_{0}=1.0\,\mathrm{Pa}$, the incompressible fluid has a density of $\rho_{0,\mathrm{f}}=1.0\,\mathrm{kg}\,\mathrm{m}^{-3}$ and the compressible gas has a density of $\rho_{\mathrm{g}}=1.0\,\mathrm{kg}\,\mathrm{m}^{-3}$ at $p_{0}$. The compressible gas has a heat capacity ratio of $\gamma_{0}=1.4$ and both fluids have a specific isobaric heat capacity of $c_{p,0}=1008\,\mathrm{J}\,\mathrm{kg}^{-1}\,\mathrm{K}^{-1}$. The viscosity $\mu$, which is the same for both fluids, follows from the considered Laplace number, $\mathrm{La}=\rho_{0,\mathrm{f}}d_{\mathrm{B}}\sigma/\mu^{2}=120$. The bubble is placed at the centre of a domain with edge length $2\,\mathrm{m}\times 2\,\mathrm{m}$, represented by an equidistant Cartesian mesh with $64\times 64$ mesh cells, and the applied time-step is $\Delta t=10^{-3}\,\mathrm{s}$, which satisfies the capillary time-step constraint [50]. The fluid interface is advected using the VOF-PLIC method.

Figure 2 shows the total kinetic energy $E_{\mathrm{kin}}=\sum_{P=1}^{N}\rho_{P}\mathbf{u}_{P}^{2}/2$, with $N$ the total number of mesh cells, in the computational domain as a function of the dimensionless time $\tau=t/t_{\mu}$, with $t_{\mu}=\rho_{0,\mathrm{f}}d_{\mathrm{B}}^{2}/\mu$ the viscous timescale. Since the bubble is in equilibrium, the observed kinetic energy is the result of parasitic currents only. The imbalance caused by errors in the numerical evaluation of the interface curvature leads to an initial production of parasitic currents and the associated kinetic energy. As the interface topology relaxes towards a numerical equilibrium, the production of parasitic currents diminishes and the existing parasitic currents dissipate as a result of viscous stresses. Consequently, the kinetic energy in the domain decays rapidly and the exact balance is recovered on the discrete level.

To demonstrate the accurate prediction of viscous flows and surface tension by the proposed algorithm, the oscillation of a standing capillary wave between two viscous fluids is simulated. This test case is particularly challenging because the wave amplitude is small and the temporal evolution of the wave amplitude, governed by the dispersion (due to surface tension) and attenuation (due to viscous stresses) of the capillary wave, is very sensitive to the implementation of the viscous stress terms, the surface tension model, numerical diffusion and spurious oscillations [65, 68]. The analytical solution for the initial value problem of a freely-oscillating capillary wave with infinitesimal amplitude, derived by Prosperetti [70], which is valid for two-phase systems in which both fluids have the same kinematic viscosity $\nu=\mu/\rho$ or in which one fluid is neglected, serves as a reference. The considered capillary wave has a wavelength of $\lambda=0.1\,\mathrm{m}$, both fluids have a kinematic viscosity of $\nu=1.6394\times 10^{-4}\,\mathrm{m}^{2}\,\mathrm{s}^{-1}$ and uniform initial temperature $T_{0}=300\,\mathrm{K}$, and the surface tension coefficient is $\sigma=0.25\,\pi^{-3}\,\mathrm{N}\,\mathrm{m}^{-1}$. The density of the two fluids are $\rho_{\mathrm{a}}=1\,\mathrm{kg}\,\mathrm{m}^{-3}$ and $\rho_{\mathrm{b}}=16\,\mathrm{kg}\,\mathrm{m}^{-3}$. The computational domain has the dimensions $\lambda\times 3\lambda$ and is represented by an equidistant Cartesian mesh with $\Delta x=\lambda/32$. A compressible-compressible flow, with $\gamma_{0}=1.4$, as well as a compressible-incompressible flow (fluid ’b’ is considered to be incompressible) are considered. The fluid interface is advected using the algebraic VOF method.

Figure 3 shows the dimensionless wave amplitude $A/A_{0}$ as a function of the dimensionless time $t\omega_{0}$, where $\omega_{0}=\sqrt{\sigma k^{3}/(\rho_{\mathrm{a}}+\rho_{\mathrm{b}})}$ is the undamped angular frequency of the capillary wave and $k=2\pi/\lambda$ is the wavenumber, with initial wave amplitude $A_{0}=0.01\lambda$, obtained for both flow types. In both cases the computed result is in excellent agreement with the analytical solution of Prosperetti [70], with respect to the amplitude as well as the frequency.

The propagation of a single acoustic wave in an air-water flow, where air is treated as a compressible fluid and water is taken to be either a compressible fluid described by the stiffened-gas model or an incompressible fluid with constant density, is simulated in a one-dimensional domain with mesh spacing $\Delta x=2\times 10^{-3}\,\mathrm{m}$, with initial velocity $u_{0}=1\,\mathrm{m}\,\mathrm{s}^{-1}$, initial pressure $p_{0}=10^{5}\,\mathrm{Pa}$ and initial temperature $T_{0}=300\,\mathrm{K}$. Air has the fluid properties $\gamma_{\mathrm{0,Air}}=1.4$ and $\Pi_{\mathrm{0,Air}}=0\,\mathrm{Pa}$, and the density at initial conditions is $\rho_{\mathrm{Air}}=1.16\,\mathrm{kg}\,\mathrm{m}^{-3}$. Compressible water has the fluid properties $\gamma_{\mathrm{0,Water}}=4.1$, $\Pi_{\mathrm{0,Water}}=4.4\times 10^{8}\,\mathrm{Pa}$ with a density at initial conditions of $\rho_{\mathrm{Water}}=1000\,\mathrm{kg}\,\mathrm{m}^{-3}$. Incompressible water has a constant density of $\rho_{\mathrm{0,Water}}=1000\,\mathrm{kg}\,\mathrm{m}^{-3}$. The acoustic wave is initiated by the inlet-velocity

$$u_{\mathrm{in}}=\begin{cases}u_{0}+\Delta u_{0}\,\sin\left(2\pi ft+\dfrac{3}{2}\pi\right)&\mbox{if $t<f^{-1}$}\\ u_{0}-\Delta u_{0}&\mbox{if $t\geq f^{-1}$},\end{cases}$$

(62)

with $\Delta u_{0}$ the amplitude and $f=2500\,\mathrm{s}^{-1}$ the frequency of the velocity perturbation. Two perturbation amplitudes are considered, $\Delta u_{0}\in\{0.01\,u_{0},100\,u_{0}\}$, to study the linear and nonlinear acoustic regimes.

In the linear acoustic regime ($\Delta u\ll a_{0}$ and $\Delta\rho\ll\rho_{0}$), which applies to the considered velocity perturbation amplitude of $\Delta u_{0}=0.01\,u_{0}$, the perturbation leads to a sound wave. Based on linear acoustic theory [61], the acoustic wave reflected at the air-water interface should have a pressure amplitude of

$$\Delta p_{\mathrm{Air},0}^{\mathrm{reflected}}=\dfrac{\Delta p_{\mathrm{Air},0}^{\mathrm{incident}}}{\dfrac{2Z_{\mathrm{Water}}}{Z_{\mathrm{Water}}-Z_{\mathrm{Air}}}-1},$$

(63)

where $Z=\rho a$ is the characteristic acoustic impedance and $\Delta p_{\mathrm{Air},0}^{\mathrm{incident}}=Z_{\mathrm{Air}}\,\Delta u_{0}$ is the pressure amplitude of the incident wave. While the theoretical pressure amplitude of the reflected wave is $\Delta p_{\mathrm{Air},0}^{\mathrm{reflected}}=\Delta p_{\mathrm{Air},0}^{\mathrm{incident}}$ if water is considered to be an incompressible fluid ($\rho_{\mathrm{Water}}=\mathrm{const.}$, $a_{\mathrm{Water}}=\infty$), the reflected pressure amplitude is $\Delta p_{\mathrm{Air},0}^{\mathrm{reflected}}\simeq 0.999\,\Delta p_{\mathrm{Air},0}^{\mathrm{incident}}$ if water is considered to be a compressible fluid ($\rho_{\mathrm{Water}}=f(p,T)$, $a_{\mathrm{Water}}=1343\,\mathrm{m}\,\mathrm{s}^{-1}$). Hence, the differences in the results between the incompressible and compressible treatments of water should be negligible, at least in the linear acoustic regime, apart from the transmitted wave that should only be present for the compressible treatment of water.

Figure 4 shows the pressure profile of the acoustic waves, after they have interacted with the interface, for both considered perturbation amplitudes using both the compressible treatment and the incompressible treatment of the liquid phase. For the small perturbation amplitude, $\Delta u_{0}=0.01\,u_{0}$, the amplitude of the reflected pressure wave is in excellent agreement with linear acoustic theory with both treatments of the liquid phase. When the liquid phase is treated as a compressible fluid, a pressure wave is transmitted through the interface, while no wave is transmitted when the liquid phase is treated as an incompressible fluid, as expected. For the large perturbation amplitude, $\Delta u_{0}=100\,u_{0}$, the pressure profile of the acoustic wave departs significantly from its originally sinusoidal shape as a result of nonlinear wave steepening. Nevertheless, the reflected waves predicted by the compressible treatment and the incompressible treatment of the liquid phase are in excellent agreement, suggesting that the compressible-incompressible treatment of gas-liquid flows can be applied even when large amplitude acoustic waves are present in the gas phase.

The pressure-driven growth and collapse of gas bubbles, in particular cavitation, is a widely observed phenomenon, and the ability to focus large amounts of energy by a collapsing bubble is being utilised in an increasing number of applications. The simplest analytical model for the pressure-driven collapse of a bubble in a liquid is the Rayleigh-Plesset equation [71] without considering viscous dissipation and surface tension, given as

$$R\ddot{R}+\frac{3}{2}\dot{R}^{2}=\frac{p_{\mathrm{g}}-p_{\infty}}{\rho_{0,\ell}},$$

(64)

where $R$ is the bubble radius, $\rho_{0,\ell}$ is the (constant) density of the liquid, $p_{\mathrm{g}}$ is the (uniform) gas pressure inside the bubble and $p_{\infty}$ is the ambient liquid pressure at infinite distance from the bubble. Interestingly, Eq. (64) is based on the assumption of an incompressible liquid, with $\rho_{0,\ell}=\mathrm{const.}$, and a compressible gas bubble.

The collapse of a spherical gas bubble due to an overpressure in the liquid, the so-called Rayleigh collapse [72], is considered to validate the prediction of pressure-driven flows of a compressible gas in contact with an incompressible liquid by the proposed algorithm. Following the setup of Denner et al. [33], a bubble with initial radius $R_{0}$ and an initial gas pressure of $p_{\mathrm{g}}=4000\,\mathrm{kPa}$ is situated in a liquid with an initial pressure of $p_{\ell}(r)=p_{\infty}+(p_{\mathrm{g}}-p_{\infty})\,R_{0}/r$, with $p_{\infty}=10^{5}\,\mathrm{Pa}$ and $r$ the radial coordinate. The gas has a heat capacity ratio of $\gamma_{0}=1.4$ and the liquid has a constant density of $\rho_{0,\ell}=1000\,\mathrm{kg}\,\mathrm{m}^{-3}$. Viscosity, thermal conduction and surface tension are neglected. The fluid interface is advected using the algebraic VOF method. Because the liquid is incompressible and viscous dissipation and thermal conduction are neglected, the oscillations of the bubble radius should continue indefinitely with unchanged frequency and amplitude.

Figure 5 shows the dimensionless radius $R/R_{0}$ as a function of the dimensionless time $t/t_{\text{R}}$, where $t_{\text{R}}=0.915\,R_{0}\sqrt{\rho_{0,\ell}/p_{\infty}}$ is the Rayleigh collapse time, obtained with different spatial and temporal resolutions using the proposed algorithm. The solution of the Rayleigh-Plesset equation (64) is also shown as a reference and regarded here as the exact solution of the dynamic bubble behaviour under consideration. The results predicted by the proposed algorithm and by the Rayleigh-Plesset equation are in excellent agreement, as shown in Fig. 5, for a sufficiently large spatial and temporal resolution.