An extension of the all-Mach number pressure-based solution framework for numerical modelling of two-phase flows with interface

In this paper, we present the extension of the pressure-based solver designed for the simulation of compressible and/or incompressible two-phase flows of viscous fluids. The core of the numerical scheme is based on the hybrid Kurganov — Noele — Petrova/PIMPLE algo-rithm. The governing equations are discretized in the conservative form and solved for velocity and pressure, with the density evaluated by an equation of state. The acoustic-con-servative interface discretization technique helps to prevent the unphysical instabilities on the interface. The solver was validated on various cases in wide range of Mach number, both for single-phase and two-phase flows. The numerical algorithm was implemented on the basis of the well-known open-source Computational Fluid Dynamics library OpenFOAM in the solver called interTwoPhaseCentralFoam. The source


I. Introduction
Compressibility of two-phase flows plays a key role in many industrial processes. For example, interaction between liquid metal jets and gas shocks has a significant impact on sizes of resulting metal granules [1,2]. Small gaseous inclusions in liquid metals (such as micro-bubbles) yield variation of acoustic and resonant properties affecting metal casting production [3]. Another curious effect is the possibility of shock appearance in a two-phase mixture when the velocity of 1 Corresponding author: Email: os-cfd@yandex.ru the flow seems to be subsonic for every pure mixture component. It is well known that speed of sound changes non-linearly and nonmonotonically with volume fraction of mixture components, which might produce shocks (see [4]). Also, phase separation processes, that occurs during discharge of supercritical fluid into low pressure volume, must be taken into account [5,6] in a fluid motion model.
Numerical simulation of such processes involves approximation of complex mathematical models, including compressible sub-and supersonic viscous flow equations coupled with real-gas equations of state. However, numerical methods differ significantly depending on flow regime and a range of dimensionless numbers (such as !", #$, %", &' and others). This diversity and inconsistency of numerical methods for different flow regimes create an obstacle to multiphysics simulations, especially for problems without dominance of a single phenomenon (viscosity, inertia, gravity, etc). For example, the ability of a numerical approach to recover behavior of the incompressible flow is important for resolving hydrodynamic instabilities and transient effects, which appear at small Mach and Reynolds number values (for example, in fluidstructure interaction problems, where a standard PIMPLE (Pressure Implicit with Splitting Operators and Semi-Implicit Pressure-Linked Equation) method have shown already its robustness [7]). On the other hand, the combination of the high-speed flow and viscous effects is a fairly common case in many problems. Therefore, a numerical scheme must be able to reproduce corresponding terms properly [8].
Obviously, the presence of two or several phases amplifies the complexity of a numerical method due to extension of a range of characteristic dimensionless numbers.
For the single phase modelling including multicomponent gases, a pressure-based numerical tool employing the hybrid approximation of convective fluxes [9,10] was developed earlier. The solver has been validated against different physical conditions [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] and for the wide range of Mach numbers. The next step for the framework development was the adaptation for simulation of high-speed multicomponent real-gas flows [27,28] in context of consideration of the phase separation phenomenon. The hybrid numerical framework was used in research by M. Pfitzner's group over last several years [5,6,27,[29][30][31]. The works were consolidated in PhD thesis "Real-Gas Effects and Single-Phase Instabilities during Injection, Mixing and Combustion under High-Pressure Condition'' [25]. It was shown that the numerical approach is able to predict an appearance of new phases from single fluid, when thermodynamic parameters pass below critical point. However, the problem of interfacial flow numerical modelling (including surface tension forces, abrupt change of properties, etc) within single model remains open.
This paper presents the next step towards generalization of the pressure-based hybrid framework using the KNP/PIMPLE or Kurganov -Noelle -Petrova/ PIMPLE [9] approach. The general idea to apply this approach for the two-phase flows was inspired by works [32,33], where it was implemented within the fractional step method for the prediction of cavitation. To preserve the consistency for thermodynamic variables and to ensure the stability of the interface modelling, the hybrid KNP/PIMPLE approximation of two-phase flow has been augmented by the following ideas: • the monotonicity of primitive variables (the velocity, the pressure and the temperature) on the interface is preserved by Acoustically Conservative Interface Discretization (ACID) technique, developed by F. Denner [34]; • the convective fluxes for each phase are approximated by original hybrid approximation method [10]; • the pressure is recovered after each iteration from the pressure equation formulated for the gas-liquid mix ture; • the motion of the phase interface is described by the volume fraction advection equation for liquid phase. The article is organized as follows. The first section recalls the governing equations of the two-phase mixture motion. The second section contains the key points of numerical approximation, based on the hybrid KNP/PIMPLE approach and ACID technique. The third section describes the validation of the described numerical algorithm on the several test problems. Main conclusions of the work are gathered in the "closure'' section.

II. Governing equations
The viscous two-phase fluid is presented as a pure mixture. The motion of the mixture is described by the system of equations, based on the Kapila's reduced model [35,36]. Kapila's model is used as the hyperbolic part of viscous flow equations, therefore, all assumptions for original model of Kapila are applicable to the following governing equations.
The system of equations comprises the mixture continuity equation: (1) the mixture momentum equation: the mixture energy equation: the phase mass conservation equation: and the liquid volume fraction transport equation: Here ( is the mixture density field, ) * *⃗ is the mixture velocity vector field, , is the mixture pressure field, σ . is the mixture viscous stress tensor field, ℎ ! = ℎ + 1 2 ⁄ 5) * *⃗ 5 " is the mixture specific total enthalpy field, ℎ is the mixture thermodynamic specific enthalpy field, 6 ⃗ is the mixture diffusive heat flux, 7 # is the 8 −th phase volume fraction field, ( # is the 8 −th phase thermodynamic density field, : $" is the interface compression coefficient field, ; ⃗ is the gravity acceleration, 8 = 1,2 is the phase index, 1 corresponds to the liquid phase and 2 corresponds to the gas phase. The system (1)-(5) is closed with the following relations: • the perfect gas equation of state for the gas phase: , = ( " ! # " ⁄ =, • the perfect fluid equation of state for the liquid phase: • the relation for the mixture compression coefficient: • the Fourier law for the heat flux: • the Newton and Stokes assumptions for the viscous stress tensor: K is the deformation rate tensor, D = 7 # D # + 7 $ D $ is the mixture dynamic viscosity, > $ is the liquid phase dynamic viscosity, > " is the gas phase dynamic viscosity, I ! is the identity tensor, B = 7 # B # + 7 $ B $ is the mixture heat conductivity coefficient, B " is the heat conductivity coefficient for k-th phase, < # and < $ are the molar weights of liquid and gas, respectively, is the acoustic impedance of k-th phase, L " = MN " ; < " ⁄ > is the sonic speed of k-th phase, N " = O (," O )," ⁄ is the heat capacity ratio of k-th phase, B %,# is the constant isobaric heat capacity coefficient of k-th phase, B ',# is the constant isochoric heat capacity coefficient of k-th phase, = is the temperature of mixture, ( (,$ is the liquid phase initial density, ! is the universal gas constant.
Mixture enthalpy ℎ is calculated as the weighted sum of phase enthalpies: ⁄ is the mass fraction of k-th phase, ℎ # = B %,# = is the enthalpy of k-th phase.

III. Computational method
The system of governing equations (1)(2)(3)(4)(5) is discretized by the Finite Volume Method (FVM) with co-located variables storage on unstructured polyhedral grids [25,37,38]. The choice of the discretization method was dictated by its flexible implementation in the OpenFOAM framework, which gives the opportunity to extend the numerical model to account for more sophisticated simulation problems.
The procedure of approximation involves five steps: the application of ACID technique to the governing equations, the finite volume approximation of balance mass-weighted equations (1,4), the application of the hybrid KNP/PIMPLE procedure to formulate convective fluxes for the approximated system, the formulation of pressure equation to close the system thermodynamically, and the approximation of the liquid volume fraction transport equation (5).

A. The ACID approach in governing equations
The ACID technique helps to suppress unphysical oscillations near the interface. The approach is the substantially analogous to the Ghost Fluid Method [39,40]: the interface between two phases might be considered as the moving internal boundary [39], which requires appropriate conditions for the flow fields (temperature, velocity, pressure). Authors of the original approach [34] have proposed to treat cells intersected by interface in a special way: • a spatial distribution of the liquid volume fraction is assumed to be uniform within the stencil of the considered cell (i. e. the cell is surrounded by imaginary cells, where thermodynamic properties are saved to be original, but a mixture composition is equal to the cell one); • the temporal change of the mixture composition is neglected. Therefore, convective fluxes become asymmetric in interfacial cells. If we consider two adjacent cells, where volume fractions are different, the material flux from the first cell to the second would differ from the reversed flux value. Application of this rule to discretized equations produces rather complex algebraic expressions.
However, the ACID idea can be applied directly to balance equations, written in the differential form. If we consider the discrete conditions imposed by ACID onto mixture composition, they can be summarized as the following finite difference relations for the liquid volume fraction 7 $ : ,! is the temporal numerical derivative.
If we substitute these expressions into the FVM approximation for the material derivative of the volume fraction, we come to the conclusion, that ACID technique set these terms to zero. Therefore, ACID technique introduces into original balances an approximation error, which is equal to subtraction of terms with the multiplier ∂ t α 1 +U **⃗ ⋅∇α 1 .

B. FVM discretization of governing equations
The further FVM approximation of the (6)-(8) system under ACID assumptions is built by summation of phase-wise discretization of balances weighted with the corresponding volume fractions. Such approach has a clear benefit: it is possible to resolve wave propagation in each phase separately, instead of considering the Riemann problem for the whole mixture. Therefore, the discrete model of two-phase mixture motion comprises the FVM approximation for modified balance equations (6)-(8) and original liquid volume fraction equation (5) (equations (9)-(12)).
In equations (9) multiplied by its normal vector K *⃗ . , M is the volume of cell, NO is the time step. If the time level (superscript) is not specified, the current available instance is used (updated iteratively). The explicit approximation for the discrete mass conservation equation of 8 -th phase reads (13).

C. Hybrid approximation of fluxes
Convective fluxes I #,. are approximated for each phase individually using the previously proposed hybrid approach [10,27].
For diffusive fluxes, the normal component of the gradient is approximated using linear interpolation [37,38]. For example, for the heat flux 6 ⃗ could be written: where P P ⁄ K *⃗ . is a numerical approximation of a face normal derivative.
The total contribution to overall mixture balance is calculated by summing divergences of flux from each phase.
An extension of the all-Mach number pressure-based solution framework for numerical modelling of two-phase flows with interface

D. Pressure equation
Usage of hybrid flux approximation together with acoustically conservative interface discretization has allowed us to employ the projection-based method PIMPLE [25] to construct an iterative algorithm for solution of discrete system (6)- (12). The final step of all such algorithms involves solving of the pressure equation, that is derived from continuity, momentum, and other equations (energy and state) to close the system and "nudge" it to mass conservation.
The pressure equation is derived by standard procedure [10,25,27], therefore, only key points are outlined here. First, solution to velocity from momentum equation (10) is represented as the sum of operator Q * *⃗ and pressure gradient R, divided by the diagonal matrix S: Then expressions for velocity (15) and phase densities (i. e. equation of state) are substituted to the equation for mixture density (9), where fluxes I #,. are already approximated using hybrid KNP [9,41] approach.
The resulting equation reads: Phase mass fluxes are recovered from the pressure equation components. After that, the mixing procedure [10] is applied to them in order to switch numerical scheme between compressible and incompressible formulations. It is necessary to stress that face interpolation of all terms . in (16) is constructed in accordance with hybrid KNP/PIMPLE approximation procedure [10,27,41].

E. Liquid phase volume fraction transport equation
The discrete transport equation (12) for liquid phase volume fraction is solved using the explicit MULES approach [42,43]. This technique guarantees boundedness and monotonicity of the solution while preserving 2 /1 order of approximation for smooth fields.
The flux I -,. = ) * *⃗ ⋅ J ⃗ . is calculated to obey mixture volumetric continuity equation. This equation is derived by summing phase mass equations (4), normalized by corresponding density ( # : Fluxes are calculated using a procedure similar to PISO/SIMPLE methods: 1 3. correction I -,. ′ is assumed to be proportional to the gradient of some correction pressure ,′: Poisson equation is formulated for ,′: where ,5 ( ,! is the explicit approximation of material derivative for ( # calculated from (16).
The final value of I -,. is recovered from the solution of equation (18) F. The overall numerical algorithm The overall solution algorithm contains the following steps. 1. Initialize variables. 2. Compute the next time value:

Store values of variables and fields from
the previous time step. 4. Predict the density of every phase individually by solving the continuity equation for each phase (13) using values of the mass fluxes from the previous time step. 5. Start the PIMPLE loop: a. Update fluid properties, which depend on temperature and density (i.e. the compressibility coefficient X # , speed of sound B # , acoustic impedance @ # , and the stiffness coefficient : $" ). b. Solve the volume fraction transport equation (12), (18). c. Assemble the matrix for the momentum equation (10)  Stability of the algorithm depends primarily on convective terms due to implicit approximation for the diffusive one. For the convective flow, two stability criteria are used: • the mixture flux stability criterion: • the front tracking stability criterion: where ] 7 is the indicator function, which determines a position of the interface: for example, ] 7 can be equal to 1 in the region where α 1 ∈(0.01,0.99) and set to zero elsewhere. Usually, the characteristic velocity-based CFL (Courant -Friedrichs -Levy) criterion is used for compressible flows. In this case, the fully implicit PIMPLE approach is employed in small Mach number regions. For large Mach number regions, the characteristic CFL number is satisfied automatically if B_ 8 < 1 is preserved. Therefore, using of stability criteria (19)-(20) is sufficient.

IV. Results of numerical tests
Numerical tests were run for the following problems to demonstrate the ability of the proposed computational method to work in different regimes (compressible/incompressible, viscous/inviscid, single-phase/twophase): • single-phase inviscid 1D Riemann problems: -Sod's problem; propagation of pressure wave through liquid; • two-phase inviscid 1D Riemann problems: movement of contact discontinuity separating two phases; interaction of high-pressured gas with normal-pressured liquid; All physical parameters and quantities were nondimensionalized where possible.

A. 1D Riemann problems
First two 1D single-phase tests help to assess the ability of the computational method to resolve propagation of waves in high-speed gas and low-speed liquid flows separately. Spatial and temporal grid convergence were studied. In all subsequent 1D problems, the case setup is similar (Fig. 1): a channel with constant cross-section area is divided by impermeable membrane initially. An initial state in the left part of the channel is determined by volume fraction of liquid 7 $,9 , temperature = 9 , pressure , 9 and velocity ) 9 ; an initial state in the right part is determined by volume fraction of liquid 7 $,: , temperature = : , pressure , : and velocity ) : . Right and left boundaries of the channel (a = 1 and a = 0 respectively) are placed far enough from the impermeable membrane to prevent interaction of propagated disturbances with boundary conditions before the end of computations. Boundary conditions for all fields are set to zero normal derivative. When a calculation starts, the membrane is removed, and initial disturbance propagates from the center to the left and to the right.

Gas shock tube
Propagation of discontinuities in 1D channel was studied. The channel is filled by perfect gas with heat capacity ratio γ2 = 1.4. Initial values for left and right states of gas are presented in Table I   Results are presented in Figs. 2-3. The solution converges and the numerical dissipation vanishes with decreasing of time step and with refinement of a spatial grid. The test demonstrates that single-phase gas dynamics is recovered in the present approximation; the behavior of hybrid KNP/PIMPLE scheme for pure single-phase solver was studied for this test in previous work [10].

Propagation of pressure wave in liquid
Propagation of acoustic wave in lowcompressible medium is demonstrated in this case. Initial conditions are presented in Table II. Compressibility of liquid X $ was set to 1/(R1T), where R1 = 6934 to get a particular speed of sound. Molar mass was chosen equal to M1 = 18×10 -3 , reference density was chosen equal to ρ0,1 = 1033, dynamic viscosity was set to zero. Results of simulations are presented in Figs. 4-5.  Qualitatively, the solution is reproduced in accordance with corresponding physical processes: shock wave and rarefaction waves travel away from initial discontinuity with the same speeds, because the temperature variation is negligible in low-compressible media. The region between two waves is occupied by constant pressure equal to the average of minimum and maximum values. Solution converges to an ideal case with spatial and temporal mesh refinement.
However, it can be seen in Fig. 5, that numerical diffusion for the case of subsonic flow is highly dependent on the Courant number. This observation demonstrates changes in the numerical scheme, that arise as it switches from KNP (Co=0.0002) to PIMPLE (Co=0.1) formulations. When the time step is large enough to allow acoustic waves travelling more than one cell per time step, an additional diffusion appears. This diffusion indicates the acoustic solution time averaging. On the other hand, numerical diffusion helps to filter out numerical oscillations.

Moving contact discontinuity
This test verifies the continuous behavior of pressure, velocity, and temperature near the constantly moving contact discontinuity, which imitates a phase interface. Two neighboring volumes of gas and liquid move with constant speed in space and time. Initial conditions for the case are presented in Table III. Liquid and gas properties were similar to previous cases IV A 1 and IV A 2 Simulation time Tend is 0.01.  Results of simulations are presented in Fig. 6. Notably, pressure, velocity, and temperature preserve continuity across the phase interfaces, even during movement.

Pressure discharge from gas into liquid
The problem studies an interaction of a gas compression wave with a liquid column [45,46]: the high-pressured gas pushes liquid. Initial conditions are presented in Table IV; liquid and gas properties are taken identical from the case IV A 3. All features of compressible flow in this scenario are presented: the rarefaction wave is going to the left into gas, the compression wave and the continuous high-speed velocity front are running to the right (Fig. 7). Mesh convergence is presented by computations on two meshes (500 and 5000 cells per channel length).

Pressure discharge from liquid into gas
The opposite to IV A 4 situation is studied: the high-pressured volume of liquid interacts with the gas, creating an expansion wave in the liquid and a compression wave in the gas [45,46]. Initial conditions are presented in Table V, liquid and gas properties are taken identical from case IV A 3. Due to relatively small velocity in the liquid region (Ma ≈ 0.002), the numerical solution in this case might demonstrate both compressible and incompressible behavior. If the time step is adjusted in accordance with flow velocity, then acoustic Co in the liquid region exceeds unity and flow is diffused additionally (Fig. 8). When the time step is adjusted to keep acoustic Co less than 0.5 in the whole domain, the sharp behavior of velocity, pressure, temperature, and other associated properties recovers (Fig. 9).
Table V -Initial state for 1D pressure discharge from liquid into gas (case IV A 5)

B. Laminar flow over backwardfacing step
Behavior of the proposed method in the laminar incompressible regime is verified using the problem of flow over backward facing step [47]. The computational domain consists of two straight horizontal channels with constant height. The short channel of the height h expands abruptly into the long one of the height H = 2h and the length L (Fig. 10). The length of the short part is equal to h and the length of the long part is L = 26h. The flow at the inlet of the channel is laminar and obeys Poiseuille's law with a parabolic velocity profile: U x (y)=U max Y1-y 2 h 2 ⁄ Z, where y is measured from the centreline of the small channel. The average velocity value is The length of the separation zone xl behind the backward-facing step is one of the main parameters characterizing the flow. When the flow is laminar and two-dimensional, the dependency of xl(Re) is almost linear. The average inlet velocity was kept constant for all calculations, while the value of Reynolds number was adjusted by variation of fluids dynamic viscosity.
The constant horizontal velocity profile is specified on the left vertical line (inlet). Fixed pressure level is prescribed on the right vertical line (outlet). Flow is assumed to be single-phased and isothermal. Horizontal walls and the rear edge of the step are treated as no-slip impermeable walls. Uniform spatial grids with 10, 20 and 40 cells per height h were used. For the last two grids, the difference in the separation zone length xl was negligible, which proves the grid convergence of the numerical solution.
The results of the newly developed solver were compared against simulations performed by the simpleFoam (standard Open-FOAM solver) in Table VI. According to the experimental data [47], the flow for Re h =U av h ν ⁄ up to 400 is two-dimensional and steady. The presented measurements clearly show the good agreement between experimental data and numerical results for both methods when Reynolds numbers up to 300. The flow streamlines are presented in Fig. 11.

C. Dam break problem
The ability of the developed solver to model complex gas-liquid interfacial flows at low Mach number conditions is validated using the dam break problem [48]. The case setup was obtained from the numerical study [49] (Fig. 12). Since the flow dynamics is driven mainly by gravity force, this example allows to demonstrate interaction of different momentum equations terms.
Snapshots of the liquid phase volume fraction are compared with experimental photos [50]. Three mesh resolutions were used to check mesh convergence: 2 cells per h value (cph), 4 cph and 8 cph. Mesh convergence is presented in Fig. 13. The numerical solution demonstrates the key features similar to experimental observations: deformation of the column in early stages, ejection of the liquid sheet, formation of the gas bubble beneath this sheet when it hits the right wall (Fig. 14).

D. Interaction between liquid droplet and blast wave
Interaction between planar blast wave and liquid two-dimensional column is studied. The case allows to validate numerical scheme for complex problem, where transient and spatial processes of high-speed gas-liquid interaction take place. A sketch of the computational domain is presented in Fig. 15. The rectangular domain is filled initially with quiescent medium at the following parameters: = " = 300: , , " = 580 × 10 ; f$ , 7 $ = 0 in the ignition zone; = $ = 300: , , " = 10 ; f$ , 7 $ = 1 in the droplet zone; = $ = 300: , , " = 10 ; f$ , 7 $ = 0 elsewhere. After the start of simulation, the disturbance near the ignition zone start to propagate towards liquid droplet at the speed, equal to 2.4 velocities of sound at ambient conditions. Pressure magnitude in the ignition zone is adjusted in such a way to match the prescribed velocity (≈ 823.2 m/s).
Experimental results and details of the problem statement are given in the original paper [51]. Comparison of numerical simulation and experimental observations measured by two sensors [51] are presented in Fig. 16; shadow photographs [51] and the pictures of numerical simulation are compared in Fig. 17. Computed data pressure-time series near sensors were averaged over time and space according to conditions of experiment [51]. It can be noted that first stages of shock wave interaction to droplet have very good agreement with experiment. However, the further time evolution of the process shows the divergence between model and observation, especially in the region, filled by fluid. These discrepancies might be related to the emergence of bubbles due to cavitation and the chosen equation of state for liquid phase (perfect fluid).

V. Summary
The approximation of the compressible two-phase flow model based on the reduced model of Kapila [35] is introduced. The numerical approximation of the system is built with the following methods: (a) hybrid Kurganov -Noelle -Petrova / PIMPLE method for convective fluxes, (b) pressure equation for ensuring mixture continuity, (c) implicit approximation for viscous fluxes, (d) explicit equation for the liquid volume fraction coupled with mixture mass fluxes through the correction equation and (e) acoustically-conservative interface discretization (ACID) technique for the mass, energy and momentum balances near the interface. The presented hybrid approximation of the two-phase system provides several advantages over other computational techniques, such as pure implicit pressurebased or explicit Godunov methods: it might be used for flows with surface tension and/or turbulence closure using conventional numerical tools, which are usually employed in 2 nd order finite volume programs; it uses monotonicity-preserving numerical schemes (KNP and ACID) for the numerical solution near discontinuities; it employs the pressure equation for the mixture mass balance, what usually impacts positively the overall robustness of a numerical scheme.
The described numerical algorithm was implemented as the OpenFOAM solver called interTwoPhaseCentralFoam. The solver was built as an extension of the previously developed hybrid algorithm for single-phase all Mach number flows. The source code is available on GitHub: https://github.com/unicfdlab/hybridCen-tralSolvers.
The solver was tested against several problems (1D Riemann problems for singlephase and two-phase flows, laminar incompressible flow over backward facing step, dam break, interaction of liquid droplet and gas blast wave). Source code of test cases is also available in the solver repository. Numerical tests have demonstrated the ability of the code to resolve both compressible and incompressible flows with phase interface adequately. The functionality of the numerical algorithm for single-phase flows was recovered in regions placed sufficiently far from the interface.