A pressure-based diffuse interface method for low-Mach multiphase flows with mass transfer

This study presents a novel pressure-based methodology for the efficient numerical solution of a four-equation two-phase diffuse interface model. The proposed methodology has the potential to simulate low-Mach flows with mass transfer. In contrast to the classical conservative four-equation model formulation, the adopted set of equations features volume fraction, temperature, velocity and pressure as the primary variables. The model includes the effects of viscosity, surface tension, thermal conductivity and gravity, and has the ability to incorporate complex equations of state. Additionally, a Gibbs free energy relaxation procedure is used to model mass transfer. A key characteristic of the proposed methodology is the use of high performance and scalable solvers for the solution of the Helmholtz equation for the pressure, which drastically reduces the computational cost compared to analogous density-based approaches. We demonstrate the capabilities of the methodology to simulate flows with large density and viscosity ratios through extended verification against a range of different test cases. Finally, the potential of the methodology to tackle challenging phase change flows is demonstrated with the simulation of three-dimensional nucleate boiling.


Introduction
Many flows of engineering interest exhibit compressibility effects, even if the flow velocities are relatively small and much smaller than the speed of sound in extended regions of the domain. Labelled as low-Mach flows, these types of flows share characteristics with both the fully compressible regime (pressure-density coupling dominates) and the fully incompressible regime (pressure-velocity coupling dominates), making their numerical simulation challenging [1]. To further complicate things, real-life applications usually involve two or more fluids that interact dynamically, with the possibility of mass transfer between different phases.
A typical example of this type of flow is boiling; the phenomenon occurring when a heated liquid reaches or exceeds its saturation temperature at a specific pressure. At that point, the vaporisation process is initiated and bubbles are formed either on the heated surface or in an adjacent liquid layer, due to the presence of microscopic surface imperfections or other impurities [2]. With continuous heat transfer, the bubbles grow and eventually detach from the heated surface, and rise inside the liquid. Even though the liquid phase can be considered approximately incompressible, the large heat and mass transfer rates induce compressibility effects in the gas phase that cannot be ignored. Additional characteristics such as the sensitivity on wetting and surface roughness render boiling as a very complex physical phenomenon, one whose physical interpretation is yet to be fully revealed. Nonetheless, boiling is utilised in steam generators [3], heat exchangers [4] and electronics cooling [5] amongst other applications. Other phase transition processes such as cavitation and evaporation can also exhibit weak or strong compressibility effects. Even though these phase transition processes are seemingly different to boiling, they are all governed by the same physical mechanism, the equilibrium of the local Gibbs free energy between the two phases [6].
One of the most important aspects of the numerical modeling of multiphase flows is the method used to describe the movement of the liquid-gas interface. This is usually done with the use of a marker function that acquires different values for each fluid, and helps identify the interface. These methods can be broadly categorised as [7] (i) surface tracking methods, where the marker function is reconstructed by marker points on the interface that are advected [8,9], and (ii) surface capturing methods, where the marker function is advected directly. Prominent examples of this second category are the level-set methods [10], the volume-of-fluid methods [11] and the diffuse interface (DI) methods [12,13]. Even though DI methods have the drawback of adopting an interface thickness that is significantly larger compared to the physical thickness, these methods have important advantages when used for multiphase compressible flows. First of all, the thermodynamic consistency is retained everywhere, even at the interface where the averaging of the properties of each phase takes place. Also the diffused shape of the interface allows the numerical resolution of the property gradients, which is very beneficial on the overall accuracy and stability of the solution methodology. Moreover, the dynamic creation and disappearance of interfaces emerges naturally, a feature that is of great importance in boiling simulations [14].
The most general two-phase DI model is the Baer-Nunziato model [15] (and the variant of Saurel-Abgrall [13]), consisting of seven equations: two equations for the mass conservation in each phase, two equations for the momentum conservation, two equations for the total energy conservation and one for the evolution of the volume fraction. This model is often characterised as a non-equilibrium model, meaning that in the regions where both phases coexist, there is no requirement for kinetic equilibrium (same velocity), mechanical equilibrium (same pressure), thermal equilibrium (same temperature) or chemical equilibrium (same Gibbs free energy). From this parent model, a hierarchy of models arises via relaxation processes that drive the system to specific equilibrium states [16], such as, • the six-equation model of Saurel et al. [17], a kinetic equilibrium model where stiff pressure relaxation is applied (see also [18,19]), • the five-equation model of Kapila et al. [20], with kinetic and mechanical equilibrium (see also [21,22,23,24,25]), • the four-equation model of Abgrall [26], with kinetic, mechanical and thermal equilibrium (see also [27,28,29,30,6]).
Amongst these models, the five-equation model provides a good compromise between physical complexity and performance, and it is the most widely used model for compressible two-phase simulations. Nonetheless, in the presence of conductive heat transfer, the additional simplification of stiff thermal relaxation is justified [30]. More specifically, when the thermal boundary layers on either side of the liquid-gas interface are properly resolved, the temperature at the interface should be continuous and a four-equation model becomes appropriate. Since the numerical methodology proposed in this study aims to simulate boiling flows where conduction heat transfer is prominent, a four-equation model will be adopted.
To numerically solve the adopted model in the low-Mach regime, two broad categories of solution strategies exist, namely (i) the density-based approach, originating directly from methodologies for compressible flows (e.g. [31,14,32]), and (ii) the pressure-based approach, originating from methodologies for incompressible flows (e.g. [33,34]). Even though density-based approaches have been shown to perform well in a range of different multiphase flows, they rely on preconditioning techniques to overcome the stiffness problem in the low-Mach limit [35]. Preconditioning techniques continue to develop and become more sophisticated, but there is still much room for improvement, especially for the simulation of unsteady flows [36] or three-dimensional cases where the computational cost of preconditioning becomes hardly feasible for practical applications. As commented in the Future Issues section of the recent review paper of Saurel and Pantano [37], preconditioning techniques for two-phase low-Mach simulations need to become more efficient. On the other hand, the pressure-based approach can achieve good performance in terms of computational cost. Moreover, it has the advantage of preventing pressure oscillations at interfaces, since the pressure is solved for and not retrieved from the energy. This approach was mainly employed in single-phase cases, with the first true all-Mach number flow solver presented in [38]. Numerical schemes specific to single-phase low-Mach number flows were also presented in [39,40,41], inspiring further development of novel pressure-based low-and all-Mach number methods [42,43,44]. Only a small number of studies were devoted to weakly compressible multiphase flows, such as [45,46] based on the Baer-Nunziato model, [47] based on the six-equation model, and other sharp interface methods [34,48,33,49]. These few studies showed promising results, both in terms of numerical efficiency and overall performance for several test cases such as bubble oscillations, explosion shocks, oscillating water column, etc. In addition, few studies demonstrated the possibility to add phase transition to pressure-based methods, using either a sharp interface approach [50,51,52,53] or the Cahn-Hilliard phase-field method [54,55].
This study presents the development of a novel pressure-based methodology for the solution of a four-equation DI model. The methodology is capable of simulating low-Mach flows, where weak compressibility effects cannot be ignored even though the flow 3 velocities are much smaller than the speed of sound. Mass transfer is also taken into account via a Gibbs free energy relaxation procedure that is activated whenever the proper thermodynamic conditions are locally met. The physical description is enriched with additional terms, which account for viscous stresses, surface tension, heat conduction, and gravity. Moreover, the method is able to incorporate complex equations of state, specific to each phase, and can handle large density and viscosity ratios. The rest of this paper is organised as follows: Section 2 presents the mathematical model, followed by the description of the proposed solution methodology in Section 3. The verification of the methodology against benchmark single-phase and multiphase cases, with and without mass transfer, is presented in Section 4, where the conservation of mass and total energy are also evidenced. Furthermore, the potential of this method to simulate the computationally demanding nucleate boiling flow is demonstrated in Section 5. Finally, Section 6 lists some possible extensions and improvements to be addressed and Section 7 concludes the study with a summary of the key findings.

Governing equations
As indicated in the Introduction, the diffuse interface model presented here uses a four-equation model describing a two-phase flow in kinetic, mechanical and thermal equilibrium [29,30,6]. This two-phase model results from the velocity, pressure and temperature relaxation of the full Baer-Nunziato model [15]. For completeness, the main steps of the derivation of the relaxed model are illustrated in Appendix A. As a convention, phasic quantities are identified with subscript k or explicitly with {1,2}, while mixture quantities bare no such identification. We will denote with a k the volume fraction of phase k, with ρ the mixture density, ρ = a 1 ρ 1 + a 2 ρ 2 , with u the velocity field, and with E the internal energy per unit volume. The pressure and temperature equilibrium four-equation two-phase model in the literature is commonly written in terms of the conserved variables a 1 ρ 1 , a 2 ρ 2 , ρ u and E = E + ρ | u| 2 2 (mixture total energy per unit volume). The formulation in terms of these variables, including here viscous stresses, surface tension, heat conduction, and gravity effects reads: Alternatively, one equation for one partial density (e.g. a 2 ρ 2 ) can be replaced by the equation for the mixture density ρ: The source terms on the right hand side of Eqs. (1), are defined as follows: where ν is the chemical relaxation parameter, and g is the Gibbs free energy.
where µ is the dynamic viscosity.
• Surface tension Σ where σ is the surface tension coefficient. This modeling follows the continuum surface force model (CSF) proposed by Brackbill et al. [56].
• Gravity G G = ρ g, where g is the acceleration of gravity.
where λ is the heat conduction coefficient and T is the temperature.
In the present work we adopt a different choice of primary variables for reasons related to the numerical discretization method which will be discussed in the next section. Our four-equation model formulation uses the governing equations for the volume fraction a 1 of one phase, the temperature T , the velocity u, and the pressure p: The momentum Eq. (10c) is recast in a different form from Eq. (1c) using the conservation of mass, Eq. (2). This is done to avoid the need of the updated density since, as shown in the algorithmic part of this study (see Algorithm 1), the density is updated after the numerical solution of Eq. (10c). In the above equations c is the speed of sound, given by [57], 5 The quantities S T , S (1) p and S (2) p are defined as: where, In the above expressions, Γ = ∂p ∂E ρ is the Grüneisen coefficient, κ p is the specific heat capacity at constant pressure and h is the specific enthalpy. The extensive heat capacity at a constant pressure for each phase is given by C pk = a k ρ k κ pk . The derivation of the form of the source terms appearing in the model formulation (10) from the conservative formulation (1) is presented in Appendix B. Let us note that in the equations (10) the source terms S T (D E +K) and S (2) p (D E +K) model viscous dissipation and heat conduction associated to the evolution of the same variables, a 1 , T and p. Note also that the contribution a ∇ · u in the volume fraction equation for liquid-gas mixtures accounts for mechanical cavitation, leading to an increase of the gaseous volume fraction in expansion regions and its decrease in compression regions. The term S (3) T ∇· u in the temperature equation leads to an increase of the temperature in compression regions and a decrease in expansion regions.
Let us remark that even though the model is fully compressible, discretizations based on the adopted non-conservative formulation (10) are not appropriate for the numerical solution of flows with strong compressibility effects. Since the total energy and mass are not solved for, and momentum is solved in a non-conservative form, the conservation of these quantities is not rigorously guaranteed at the discrete level. As a consequence, highly compressible flows such as those involving shock waves cannot be captured accurately. Nonetheless, as shown in the following sections, the model is appropriate for the simulation of weakly compressible flows with phase change, which are the applications targeted by the present study.

Equation of state
To close the system of equations above we need to specify an equation of state (EOS) for each phase, for instance by providing a pressure law p k (E k , ρ k ) and a temperature law T k (p k , ρ k ). Given the equation of state of each phase, the equation of state for the mixture is determined by the pressure and temperature equilibrium conditions p 1 = p 2 = p, T 1 = T 2 = T , by the mixture density relation ρ = α 1 ρ 1 + α 2 ρ 2 and by the mixture The thermodynamic closure used in this study is the Noble-Abel stiffened-gas (NASG) EOS [58]. The use of this specific EOS is not in any way mandatory, since the proposed method can be coupled with any complex EOS. The NASG EOS is expressed by, Since the pressure and temperature fields are solved for using Eqs. (10d) and (10b), Eq. (14a) is used to calculate the phasic densities ρ k , by setting p k = p and T k = T (pressure and temperature relaxation). Following the above relations, the phasic specific entropy s k , specific enthalpy h k , Gibbs free energy g k , speed of sound c k and Grüneisen coefficient Γ k become: where κ vk (phasic specific heat capacity at constant pressure), γ k , η k ,η k , p ∞ and b k are the phasic parameters of the EOS. It is noted that κ pk = γ k κ vk . The set of adopted values will be presented separately for each test case. Given the equations of state for the liquid and vapor phases of a species, the theoretical pressure-temperature saturation curve is determined by the Gibbs free energy equilibrium condition g 1 = g 2 . For liquid and vapor phases governed by the NASG EOS [58] this gives the following equation defining the p-T saturation curve: where, The parameters of the NASG EOS of the two phases are defined so that the theoretical saturation curves fit the experimental ones of the chosen material in a certain temperature range [58].

Numerical methodology
This section details the solution of Eqs. (10a)-(10d) on a staggered (marker and cell) Cartesian grid, where all scalar fields are defined on cell centers while the velocity is defined on cell faces. To handle numerically the mass transfer source term M we use an operator splitting technique, which is commonly employed to treat relaxation terms in Baer-Nunziato type models (e.g. [13,59]): the equations are first solved without the mass transfer term M , and then a relaxation procedure is applied to integrate this source term accounting for phase transition.

Pressure-based solution method for the system without mass transfer term
The equations are integrated in time using an explicit, 3 rd order Runge-Kutta (RK3) method [60]. Within the context of the RK3 method, each time-step n is split into three sub-steps m = 1, 2, 3. With this notation, the following convention is followed to represent any quantity q at the beginning and the end of each time-step, q n,m=1 = q n and q n,m=4 = q n+1 .

Temperature
Similar to the volume fraction, the solution of Eq. (10b) (without the mass transfer term) is advanced in time as, For consistency, the same spatial discretisation schemes used for Eq. (18) are also applied to Eq. (20).

Predicted velocity
To decouple velocity and pressure, a fractional-step approach is adopted (in the spirit of [63]), where a predicted velocity field u n,m * is first calculated without considering the pressure gradient term. In this form, the predicted velocity can be obtained as, Following the discretisation used for the volume fraction and temperature equations, the convection term ∇ · ( u ⊗ u) is discretised using the van Leer flux limiter, while all other terms are discretised with central differences.
Once the updated pressure field p n,m+1 becomes available, the corrected velocity field u n,m can be obtained as,

Pressure solution
Eq. (10d) (without the mass transfer term) is discretised in time as, Term ∇ · u n,m+1 is replaced by the divergence of Eq. (22), yielding the following Helmholtz equation for the pressure, In the present study, the Helmholtz equation is solved using the parallel semicoarsening multigrid (PFMG) solver combined with a Red-Black (RB) preconditioner, both available in the HYPRE library [64]. Once p n+1 is obtained, the corrected velocity field u n+1 is updated using Eq. (22).

Phase transition solver
As mentioned above the mass transfer is treated via a (first-order) operator splitting method: we first solve the system without the source term M and then we solve a system of ordinary differential equations accounting for the mass transfer. Hence we consider: p M ] T is the vector of the source terms in Eq. (10) containing the mass transfer term M = ν(g 2 − g 1 ). As one can easily see from this system of ODEs, during this chemical relaxation process the mixture density, mixture energy and velocity remain constant. To determine the state after the mass transfer step we need to determine three independent variables, for instance a 1 , T and p. We have implemented two different relaxation techniques proposed in the literature to determine the updated values of a 1 , T and p after mass transfer.
The first one [30,19] assumes instantaneous chemical relaxation, ν → ∞, so that thermodynamic equilibrium is instantaneously attained. In this case we do not need to solve a system of ODEs, but instead impose directly the equilibrium condition g 1 (p, T ) = g 2 (p, T ). This gives an algebraic system of equations to be solved for the thermodynamic equilibrium state with the mixture relations ρ = a 1 ρ 1 (p, T ) + a 2 ρ 2 (p, T ) and The second relaxation procedure, based on [65,66,67], allows the modelling of chemical relaxation of arbitrary rate, finite-rate (for instance with a given function to define ν) or instantaneous. It is based on the idea of approximating the relaxation process toward the equilibrium g 1 = g 2 by an exponential behavior. Within this approximation, a semi-exact exponential solution of the system of ODEs ∂ t [a 1 , T, u, p] T = Φ M can be found. This approximate solution is used to define the solution after the mass transfer step. This second approach is simpler since we update the variables using explicit formulas, whereas in the first approach we need the solution of a non-linear system of algebraic equations. We typically use this second relaxation technique, nonetheless the first method was also tested within the context of the present study (where we assume instantaneous mass transfer in all the tests) with very similar results.
Let us note that the temperature is equal to its saturation value T sat (p) at the equilibrium g 1 (p, T ) = g 2 (p, T ). Denoting here with subscript 1 the liquid and with subscript 2 the vapor, if g 1 > g 2 then T > T sat (P ) and liquid-to-vapor transition occurs (evaporation), whereas if g 1 < g 2 then T < T sat (P ) and vapor-to-liquid transition occurs (condensation). This mass transfer processes modelled by the chemical relaxation term may be activated and deactivated in the numerical model depending on the desired criteria. In the present study, mass transfer is activated only in the presence of a two phase mixture, i.e. when both a 1 > and a 2 > , with = 10 −8 . Moreover, in some numerical tests we may activate chemical relaxation only if the condition for evaporation T > T sat (P ) is met, as done for various tests for instance in [68,69,19]. This criterion is applied in the test taken from these references in Section 4.4. In such case, for the NASG EOS we solve the equation (16) using p n,m+1 as an independent variable, which provides T n,m+1 sat .

Algorithm overview and additional remarks
To make the proposed methodology as clear as possible, a step by step description of the overall solution procedure is presented in Algorithm 1. For the purposes of the present study, the algorithm was implemented using the framework already available in [70], with the substitution of the fast Fourier transform library with the Hypre library.
Algorithm 1 Overall solution procedure of the proposed methodology.

13:
Helmholtz Eq. (25) is solved and p n,m+1 is obtained. 14: u n,m+1 is calculated from equation (22). 15: T are updated. 16: if phase change conditions then 17: a n,m+1 1 , T n,m+1 and p n,m+1 are locally modified following the relaxation procedure of [65,66,67]. 18: a and S Since the adopted four-equation model and the proposed solution methodology are significantly different from what was previously used in multiphase diffuse interface studies, a few key points on various aspects of the methodology are discussed below: • Set of equations: Common formulations of the four-equation model in the literature adopt either (a 1 ρ 1 , a 2 ρ 2 , ρu, E) [26,29] or (a 1 ρ 1 , ρ, ρu, E) [30] as the set of primary variables (see the set of Eqs. (1), (2)). In the present study, the choice to follow a pressure-based methodology was made to avoid pressure oscillations at interfaces when the pressure is retrieved from the energy. To construct a pressure-based methodology using either of these sets of equations, the total energy equation can be transformed to a pressure equation, and the temperature can be calculated from the EOS. This approach was followed in [43] for single-phase problems, using a pressure equation similar to Eq. (10d). Using this approach, the authors of the present study noted that the calculation of the temperature field was not very accurate, especially at the interface. Even a small error in the calculation of the temperature field in phase transition simulations could cause the temperature to artificially exceed or fall below the saturation temperature, severely affecting the results. For this reason, we adopted (a 1 , T, ρu, p) as the primary variables, so to have a more accurate calculation of the temperature field. As mentioned earlier, this model is not appropriate for the numerical solution of fully compressible flows, because there are no discrete equations that guarantee the conservation of the total energy and mass, while the momentum is solved using a non-conservative form. Nevertheless, the following sections demonstrate the validity of this model in simulations of compressible flows with phase change at low speeds.
• Mass and energy conservation: Since the adopted four-equation model does not consider equations for the conservation of mass and energy, these quantities are not automatically conserved. Nonetheless, with proper spatial and temporal resolution, mass and energy are indeed conserved over long simulation times. This observation will be demonstrated and quantified in the verification cases presented in Section 4.
• Momentum discretisation: As described in Section 3.1.3, the discretised momentum equations adopt a non-conservative form of the advection term. The reason behind this is to avoid invoking the updated density field which is yet to be computed at this point. An alternative treatment (not used in the present study) would be to over-constrain the system of equations by solving an extra equation for the mass conservation before solving the momentum equation. In that case, mass will be conserved by definition and the updated density field would be available to be used in a fully conservative momentum equation. This of course would add the cost of having to solve an additional equation and encounter some loss of consistency between the various thermodynamic quantities because in that case ρ = a 1 ρ 1 (p, T ) + a 2 ρ 2 (p, T ). A similar treatment of over-constraining the system for algorithmic purposes was adopted in sharp-interface formulations [33,48].
• Time step restrictions: It is generally accepted that explicit pressure-based methods bare overwhelming time-step restrictions. In the proposed methodology, this is overcome by using the updated pressure field in the momentum equation, that results in a Helmholtz equation (Eq. (25)) for the pressure field. In addition, the adoption of the RK3 method for numerical integration improves the overall time step restrictions [60]. The time-step restriction employed in this study is [71], where ∆t c , ∆t σ , ∆t µ and ∆t λ are the maximum allowable time steps due to convection, surface tension, momentum and thermal energy diffusion. These are determined as suggested in [71]: where |u i,max | is an estimate of the maximum value of the ith component of the flow velocity, ρ k,min is the minimum density of phase k in the domain and ∆x, ∆y, ∆z are the grid spacings along the x, y, z directions. Since this study considers only weakly compressible flows, the acoustic time-step restrictions are not taken into account. Setting C ∆t = 0.25−1.0 was seen to be sufficient for a stable and accurate time integration. For a more accurate and fair comparison against reference results, a constant time step was adopted in some of the test cases. This is clearly specified in each test case.

Verification
The methodology presented in the previous sections will be verified in a number of different test cases, under incompressible and compressible conditions, with and without mass transfer. When mass transfer is activated it is always assumed as an instantaneous process (chemical relaxation parameter ν → +∞). In the following sections, wherever a two-phase mixture is present, subscript (1) refers to the liquid phase while subscript (2) refers to the gas phase.

Gresho vortex
The proposed numerical algorithm is first tested against the Gresho vortex benchmark [72], a rotating vortex which is a time-independent solution of the incompressible Euler equations. The aim of this test is to assess the accuracy of the method against an exact solution and its ability to preserve the vortex structure for different Mach numbers. To this purpose, we employ a variant of the original benchmark, where the analytical solution, reported in Eqs. (28) depends on a reference Mach number [73] and it is a continuous differentiable function [74]: where u φ is the angular velocity and u r a reference velocity. In this steady configuration, the pressure gradient is balanced by the centrifugal force and the density is uniform and equal to a reference value, i.e. ρ = ρ r . Therefore the radial momentum balance reads, By splitting the pressure p into a reference pressure p r and a second order pressure p (2) , i.e., p = p r + p (2) M a 2 r and using Eq. (29), p r results to be a uniform and constant field, while p (2) can be computed as, where M a r = ρ r /(p r γ r ) is the reference Mach number. The EOS parameters used for the calculation of the reference Mach number and all other necessary quantities are listed in Table 1, modelling an ideal gas. Using Eq. (30) finally provides an expression for p, The integrals in Eq. Note that in Eq. (31) the reference pressure field p r is given by p r = ρ r u 2 r /(γ r M a r ) (with ρ r = 1, u r = 1), r is the radial coordinate, given by r = (x − l r /2) 2 + (y − l r /2) 2 and l r a reference length. Since Eq. (28) is formulated in a polar reference frame, a coordinate transformation is performed to obtain the Cartesian velocities components, i.e. u(x, y) = u φ sin(θ) and v(x, y) = v φ cos(θ) with θ = arctan 2(y − l r /2, x − l r /2).  28) and (31), is prescribed as initial condition for three different Mach numbers, M a r = 10 −1 , 10 −2 and 10 −3 . Simulations are conducted up to tu r /l r = 2 (i.e., one complete revolution of the vortex) using a constant time-step ∆tu r /l r = 2.5 × 10 −3 . Note that this value represents the maximum allowable time-step to ensure a stable time integration for the highest grid resolutions cases (i.e., 128 × 128) and is employed for the coarser cases, irrespective of M a r . Fig. 1 shows the Mach number distribution for the various cases considered. It is clear that the proposed numerical methodology is able to preserve the vortex shape regardless of the employed M a r . This result is possible thanks to the implicit treatment of the acoustic part of the pressure field (using a prediction-correction approach [75]), which 15 γ ηη p ∞ b κ v 1.4 0 0 0 0 717.5 Table 1: EOS parameters adopted for (a) the "Gresho vortex" test case (Section 4.1) and (b) the thermally driven flow in a differentially heated cavity (Section 4.2). ensures a stable and bounded solution of the pressure equation, even for M a r → 0. The excellent ability of method in preserving the vortex shape is reflected also in the good conservation property of the kinetic energy E k = 1/2 V ρ u · udV , whose temporal evolution is reported in Fig. 2(a) for different grid resolutions. Note that E k is conserved at the highest resolution cases with M a r = 0.001 and a similar behaviour has been observed also for M a r = 0.01 and 0.1 (not shown).
We conclude the analysis with the accuracy assessment in terms of the L 1 -norm, evaluated as: where s represents either one of the Cartesian velocity components, u or v or the pressure, p, while s ex represents the corresponding exact fields computed with the analytical solution. If L 1,N is the L 1 -error using N x × N y grid point and L 1,2N is the L 1 -error evaluated with 2N x × 2N y grid points, the order of accuracy n L1 is computed as: Both L 1 -error and n L1 are evaluated at tu r /l r = 1 and the results are reported in Fig. 2(b) for M a r = 0.001. As expected, a second-order accurate solution for (u, v, p) is achieved for all cases (a similar trend has been observed also for M a r = 0.01 and 0.1), confirming the correct behaviour of the proposed method irrespective of M a r .

Thermally driven flow in a differentially heated cavity
In this section, the flow of air in a closed two-dimensional square cavity with heated and cooled side walls and adiabatic horizontal walls is considered. The ascending and descending buoyant currents next to the heated and cooled walls form a circulation current, while the central region of the cavity features an almost stagnant fluid with stratified temperature. In thermally driven flows, the compressibility effects are directly related to the temperature difference between the thermally active side walls. For temperature differences less than approximately 30 K the flow of air is considered incompressible [76], and the Oberbeck-Boussinesq approximation is typically adopted [77,78,79]. Outside the limits of the Oberbeck-Boussinesq approximation, weak compressibility effects appear, and different low-Mach methodologies were utilised to study these effects [80,81,82]. For the purposes of this study, this test case helps to assess the ability of the methodology to accurately incorporate the effects of viscosity and thermal conductivity.
The case simulated here follows the setup presented in [80,81] and involves a temperature difference ∆T = 720 K, around a reference temperature of T r = 600 K. The reference thermophysical quantities and the height of the cavity L are chosen such that the Rayleigh and Prandtl numbers are equal to, Ra = gρ 2 r κ p β∆T L 3 µλ = 10 6 , and Pr = µκ p λ = 0.71, where β is the thermal expansion coefficient. The ideal gas EOS is used as the thermodynamic closure (parameters listed in Table 1) and the reference density value ρ r is calculated using the reference temperature T r =600 K and the reference pressure p r =101325 Pa. All other thermophysical quantities are considered constant. Furthermore, no-slip boundary conditions are applied to the walls. Initially the air in the cavity is stagnant and isothermal with T (t = 0, x, y) = T r , and a hydrostatic pressure field with p(t = 0, x = L/2, y = L/2) = p r is applied. The time step is dynamically adjusted according to Eq. (26), with C ∆t = 0.5. Even though this is a simple configuration, the introduction of large temperature difference increases the resolution requirements due to the thinning of the thermal boundary 17 layer at the cooled wall [83]. In [81], a non-uniform grid of 256 × 256 was used, with clustering of grid points next to the walls, while in [80] a uniform grid of size 512 × 512 was chosen. In the present study, four different uniform grid sizes were used, namely 64 × 64, 128 × 128, 256 × 256 and 512 × 512 grid nodes. As shown in Fig. 3(a), the highest resolution considered is able to achieve conservation of the total mass with approximately 1% error. The heat transfer rate inside the cavity is expressed through the Nusselt number, defined as, where h is the heat transfer coefficient, ∇T w is the temperature gradient on any of the thermally active vertical walls andn w is the corresponding unit normal vector on the wall. The temporal evolution of the Nusselt number at the heated wall is plotted in Fig. 3, for different grid sizes. After a steep drop during the initial stages of the simulation, the Nusselt number increases gradually to a steady state value. All grid sizes considered capture the evolution of the Nusselt number fairly well, converging to the reference solution from [81] with increasing resolution. More specifically, the steady state solution of the 512 × 512 grid differs by 1.3% with respect to the reference solution.
Noting that the solution methodology followed in the reference study is based on the low-Mach number asymptotic expansion of the Navier-Stokes equations, and therefore is significantly different to the present method, the agreement between the two solutions is considered satisfactory.

Rising bubble
The well-established "Rising bubble" test case is employed here as a two-phase numerical benchmark [84]. This test case helps to assess the ability of the proposed numerical methodology to capture topological changes of a moving interface, in the presence of surface tension. More specifically, the evolution of the shape, position and velocity of the center of mass of a rising bubble in a two dimensional liquid column will be compared against the reference solution in [84].
By introducing a reference length l r , a velocity u r , gas and liquid densities ρ g,r and ρ l,r and gas and liquid dynamic viscosity µ g,r and µ l,r , we can define the five dimensionless groups which governs the problem: the Reynolds number Re = ρ g,r u r l r /µ g,r , the Weber number W e = ρ g,r u 2 r l r /σ with σ equal to the surface tension, the Froude number F r = u 2 r /(| g|d 0 ), the density ratios λ ρ = ρ l,r /ρ g,r and the viscosity ratio λ µ = µ l,r /µ g,r . The liquid column has a dimension l x = 2d 0 and l y = 4d 0 where d 0 is the bubble diameter, whose initial center of mass position is (x c,0 , y c,0 ) = (d 0 , d 0 ).
This section considers two different rising bubble test cases. For the first test case, simulations are conducted setting l r = d 0 , u r = | g|d 0 , Re = 35, W e = 1, F r = 1, λ ρ = 10 and λ µ = 10 in a domain discretized with N x ×N y = [32×64, 64×128, 128×256] grid points. Note that the EOS parameters reported in Table 2 are chosen to match the specific λ ρ . The top and bottom boundaries are no-slip non-moving walls, while periodic conditions are prescribed in the horizontal directions. The initial velocity field is zero and the pressure is uniform. A constant time step ∆t | g|/d 0 = 3 × 10 −4 is used. This value is the maximum allowable time-step to ensure a stable time integration for the highest grid resolutions cases (i.e., 128 × 256) and is employed for the coarser cases in order to ensure that the same time discretization error is introduced in all the cases.  First, Fig. 4 shows the position of the interface at different time instances, providing a qualitative assessment of the numerical solution for the three grid resolutions. An excellent agreement is observed between the numerical solution on the most refined grid and the reference solution in all the analysed time instances. Similarly, an excellent agreement is confirmed by comparing the temporal evolution of the bubble's center of mass position and vertical velocity, reported in Fig. 5. Note that the proposed benchmark is often employed for assessing the accuracy of incompressible two-phase codes without temperature variation in the bulk regions of the phases. In the present work, since the two phases exhibit compressible effects, local temperature variations cannot be avoided a-priori given the local variation of the pressure field. However, the maximum temperature variations remain below ∆T /T r ≈ 10 −5 ; confirming the excellent behaviour of the proposed numerical method in the incompressible limit.
The second test case considers larger property variations and involves extreme events such as small bubble breakups. More specifically, the parameters considered in this second benchmark are Re = 3.5, W e = 0.125, λ ρ = 1000 and λ µ = 100, while all other simulation parameters remain the same as the first test case. As before, the EOS  parameters reported in Table 2 are chosen to match the specific λ ρ considered here. Three uniform grid resolutions are used to simulate the flow: 128 × 256, 256 × 512 and 512 × 1024. First, the comparison of the interface position at dimensionless time t | g|/d 0 = 3 is presented in Fig. 6, for the three different grid resolutions considered. In this case, we observe the formation of two trailing bubble tails instead of two small circular bubbles that break away from the original structure in the reference solution. As a result, the bubble rise velocity is lower than the reference results, as depicted in Fig. 7. This deviation is reduced as the grid becomes finer, reaching a value of approximately 9% at dimensionless time t | g|/d 0 = 3, for the 512 × 1024 grid. Even though this method is able to accurately simulate flows with large property variations (as demonstrated in Sections 4.4 and 4.5) and capture large bubbles that break away from larger structures (demonstrated in Section 5), it requires a fine grid resolution to capture the break away of small bubbles. This drawback is solely attributed to the inherent diffusive nature of the diffuse interface method. Small bubbles that break away cannot get fully detached from the bigger structures without additional treatment of the interface thickness. Section 6 describes several approaches that can be incorporated in 20  the proposed methodology to manage the interface thickness.

Water liquid-vapor expansion tube
The one-dimensional water liquid-vapor expansion tube test case was first proposed in [68] and was later adopted in other studies such as [69,19]. This test case is presented here to verify the code in a fully compressible flow and to demonstrate the accuracy of the phase transition solver. A tube of L x = 1 m length is filled with a two-phase mixture of liquid water with a uniformly distributed small amount of vapor a 2 = 0.01. The tube is open at both ends and the water is subject to atmospheric pressure, p = 10 5 Pa, with a liquid density of ρ 1 = 1150 kg m −3 . Initially, a velocity discontinuity is present at the centre of the cavity x = 0.5 m, with u = −u 0 for x < 0.5 m and u = u 0 for x > 0.5 m, where u 0 = 2 m s −1 . Viscosity, thermal conductivity and surface tension effects are not considered. The EOS parameters adopted for this test case are listed in Table 3. Using this parameters, the temperature of the two-phase mixture is calculated T = 354.728 K and the vapor density ρ 2 = 0.63 kg m −3 .
Two cases were simulated: (i) phase transition is not activated and (ii) phase transition is activated when the condition T > T sat (p) is met. A uniform grid with 1024 points is used and a constant time step ∆t = 3.2 × 10 −6 s is adopted. To avoid numerical instabilities, the initial velocity discontinuity used in the present simulations is approximated via a hyperbolic tangent function in the form, where u 0 = 2 m s −1 and a value ofc = 10 −2 m is adopted for the sharpness parameter. The flow was allowed to develop for 1000 time steps, and the relevant fields at t = 0.032 s are used for comparison against the reference results from [19]. The comparison between the present and reference results is shown in Fig. 8. Starting from the case without mass transfer, two rarefaction waves appear, moving in opposite directions. The induced small mechanical expansion of the vapor phase causes an increase in vapor volume fraction at the center of the tube. Due to the presence of the rarefaction waves, the liquid phase becomes metastable and the inclusion of mass transfer influences the solution significantly. More specifically, in the presence of mass transfer, the expansion of liquid water causes the decrease of the pressure at the centre of the tube, reaching its saturation value. Consequently, a small amount of vapour is generated, as evidenced in the increase of the vapor volume and mass fractions. These characteristics provide the basis of the comparison depicted in Fig. 8, where the agreement between present and reference results is verified with only minor discrepancies. Overall, the proposed methodology accurately captures the solution of this compressible flow with and without mass transfer.
Since the equations for the conservation of mass and energy are not explicitly solved within the context of the proposed methodology, the correct calculation of these quantities needs further assessment. Assuming that the rarefactions do not reach the boundaries, the exact total mass M tot and total energy E tot in the domain as a function of time can be calculated as, The ratios between these quantities and their corresponding initial values are shown in Fig. 9, for the case where mass transfer is activated. An excellent agreement is observed between the computed and analytically derived quantities.

Water liquid-vapor filled tube with a superheated region
The test case presented in this section was specifically designed to verify the code in conditions that are more relevant to boiling flows. A 1.7m one-dimensional long tube is filled with a water mixture, with a liquid volume fraction of a 1 = 0.9 for x < 0.68 m and a 1 = 0.1 for x > 0.68 m. The tube is open at both ends and the mixture is subject to a pressure of p = 3.0104051 × 10 6 Pa. The gas-dominated region x > 0.68 m is set at saturation conditions, while the liquid-dominated region is overheated by 2.0K. The EOS parameters adopted for this test case are listed in Table 4. Using these parameters, the temperature for x > 0.68 m is set at T = T sat = 513.21628 K, while for x < 0.68 m is set at T = 515.21628 K. Viscosity, thermal conductivity and surface tension are not considered.  Table 4: EOS parameters adopted for the study of water liquid-vapor filled tube with a superheated region. To assess the accuracy of the results produced by the present methodology, the same test case was also simulated by using an established methodology for compressible twophase flows with phase transition, documented in [19,65,66,67]. The reference methodology (formally second-order accurate) is based on a six-equation two-phase diffuse interface model, and it uses a HLLC-type Riemann solver for the homogeneous equations together with mechanical, and thermo-chemical relaxation procedures for inter-phase processes. A grid of 1024 points and a constant time step of 6×10 −6 s was used to produce the present results, where the temperature and volume fraction discontinuities were approximated with a hyperbolic tangent function in the form of Eq. (36). The reference results were obtained with a grid of 5000 points and a varying time step with fixed Courant number equal to 0.5, without considering any smoothing of the initial discontinuities. For this numerical test, chemical relaxation is activated everywhere in the domain.
The comparison between the present and reference solution at t = 0.006 s is shown in Fig. 10. The pressure on the liquid side (x < 0.68 m) almost instantaneously jumps to the saturation value. Therefore a pressure discontinuity is generated, giving rise to two waves that propagate to induce a thermodynamic equilibrium between the two regions. Overall, a good agreement is observed between the results of the two numerical methods. The reference results display larger numerical diffusion, mainly due to the use of a dissipative HLLC-type Riemann solver compared to the Van-Leer flux limiter used in the present work. Some numerical diffusion is also added due to the fact that the six-equation numerical model of [19] needs two additional relaxation steps (mechanical and thermal) in the fractional step algorithm after the solution of the homogeneous equations. The present method, therefore, provides a more accurate representation of the sharp variations, despite the smoothing applied to the initial temperature and volume fraction fields.

Three-dimensional nucleate boiling in water
To demonstrate the potential of the proposed methodology to simulate challenging boiling flows, a three-dimensional nucleate boiling simulation was carried out. for this case was previously presented in [30], albeit in two-dimensions. A closed cuboid box of dimensions L x ×L y ×L z = 7×7×12 cm is filled with a water liquid-vapor mixture, with a 1 = 0.9999 for z < 6 cm and a 1 = 0.0001 for z > 6 cm. A schematic representation of the domain is shown in Fig. 11. All the domain boundaries are considered solid walls, therefore no-slip boundary conditions are applied for the velocity field. The vertical and the top walls are considered adiabatic ( ∇T | w ·n w =0), while a time-varying temperature is applied on the bottom wall, in the form, where T sat = 372.74 K is the saturation temperature at z = 0 and ∆T = 15 K. In this way, the bottom wall is slowly heated beyond the saturation temperature of the surrounding liquid, without causing too strong pressure waves. In addition to the threedimensional representation, a more accurate EOS is used here in accordance to the NASG framework, as shown in Table 5.
The acceleration of gravity is set to | g| = 9.81 m s −2 , the surface tension coefficient σ = 0.073 N m −1 while the thermal conductivity in each phase λ c,1 = 0.6788 W m −1 K −1 and λ c,2 = 0.0249 W m −1 K −1 . Similarly to [30], the effects of viscosity are not considered. Initially, the two-phase mixture is stagnant and the pressure in the cavity is 101325 Pa. The temperature is equal to the saturation value based on the local value of pressure. A uniform grid of N x × N y × N z = 256 × 256 × 512 is used (approximately 33.6 million grid points), while the time step is dynamically adjusted in accordance to Eq. (26), with C ∆t = 0.5.
At the very early stages of the flow development a vapor film starts forming at the bottom wall, which soon becomes unstable and breaks up. Instantaneous snapshots of a 1 = 0.5 iso-surfaces after the initial vapor film breakup are shown in Fig. 12. The film gives way to a torus-like structure centered around the vertical axis at the center of the cavity. In addition, bubbles are formed at the four bottom corners of the cavity and other smaller structures along the bottom edges of the cavity. As the flow develops further, the torus structure breaks up into four large bubbles along the x − y diagonals. Eventu-  ally, bubbles of different sizes reach the interface and release their vapor content to the vapor-filled top half of the cavity. Overlooking the complexity of the three-dimensional structures, this phenomenological description is similar to what was reported in the reference two-dimensional study of [30], where the initial vapor film breaks up in an elongated bubble at the center of the cavity and two smaller bubbles at the two bottom corners. Figure 12: Instantaneous snapshots of a 1 = 0.5 iso-surfaces, coloured using the temperature field for the three-dimensional nucleate boiling in water. The first snapshot is at t ≈ 0.262 s, and each subsequent snapshot at intervals ∆t ≈ 0.046 s.
Even though an in-depth investigation of this specific nucleate boiling case is outside the scope of the present study, it is evident that the proposed method can provide reliable results for such a challenging physical phenomenon. Furthermore, the current implementation is efficient enough to solve problems on high-resolution numerical grids over long integration times. More specifically, the simulation presented in this section was carried out on 1024 processors on a system based on Intel Xeon Gold 6130 CPU's for ten days, consuming approximately 0.25 million core hours.

Future improvements
In this section, a list of future modifications is compiled and motivated, with the potential to improve different aspects of the proposed methodology. The following improvements can help to enrich the physical description of the adopted model: • Wall treatment: During nucleate boiling next to a heated surface, a number of physical mechanisms are responsible for transferring energy to the forming bubble. Depending on the conditions in which boiling takes place, energy can be transferred to a growing bubble through the micro-layer (thin liquid layer trapped between the bubble and the wall) and the three-phase contact line, amongst other mechanisms [85]. Both these mechanisms act on scales that are orders of magnitude smaller than the resolution requirements for the other physical mechanisms that affect the flow. Therefore, the appropriate modeling of these mechanisms (e.g.using the models presented in [86,87,88,89]) is very important for the accurate representation of flows involving bubble nucleation close to a wall.
• Realistic EOS: Simple EOS are very convenient because they can be easily manipulated and included in the numerical method in an analytical form. On the other hand, realistic EOS for industrial applications such as the IAPWS Industrial Formulation for Water and Steam [90] are very complex and their use presents significant challenges, both in terms of thermodynamical consistency of the numerical method and computational efficiency. For such complex equations of state the presented numerical model could be coupled to a table-look up method as proposed for instance in [91,92] to achieve fast and accurate thermodynamic calculations.
Furthermore, additional numerical techniques can be incorporated to improve the accuracy and efficiency of the computational methodology: • Managing interface thickness: Even though the thermodynamic consistency is retained on the vapor-liquid interface, the interface typically becomes thicker with time in DI simulations. This inherent drawback of diffuse interface methods hinders the detachment of small bubbles from larger structures, as discussed in Section 4.3.
To treat this issue, two main approaches exist: (i) use of interface compression techniques, where carefully constructed source terms (also called regularization terms) are introduced into the equations (e.g. [24,25]), and (ii) construct a sharper colour function from the more diffused mass or volume fractions and use this sharper function to calculate interface terms such as surface tension (e.g. [30]). Employing either approach can help retain the thickness of the interface at an appropriate size.
• Improved RK3 method: The solution of the Helmholtz equation carries the biggest computational cost compared to the other algorithmic tasks. In its present form, the proposed methodology invokes the Helmholtz solver every RK3 sub-step, i.e. three times per time step. By attempting to extend the treatment of [93] and [94] to the pressure-based formulation of the present study, the Helmholtz problem will be solved only once per time step, reducing the overall computational cost of the numerical solution significantly.

Conclusions
In this study, a novel pressure-based methodology has been presented for the solution of a four-equation two-phase diffuse interface model, capable of solving low-Mach flows with mass transfer processes. The four-equation model results from the kinetic, mechanical and thermal relaxation of the general seven-equation Baer-Nunziato model, with the addition of extra terms to account for the effects of viscosity, surface tension, thermal conductivity and gravity. Mass transfer is modeled as a Gibbs free energy relaxation term.
The key characteristic that makes the proposed methodology stand out from the current state of the art is its pressure-based nature, which results in the solution of a Helmholtz equation for the pressure. This feature allows the utilisation of scalable and efficient solvers, able to exploit the full potential of high performance computing systems. With such capabilities, complex cases can be simulated with unprecedented resolution, giving a new insight into the underlying physical mechanisms.
The methodology was verified in a number of different cases, involving single-and twophase configurations with large density ratios, under both compressible and incompressible conditions, with and without mass transfer. In addition to a very good agreement reached with relevant reference data, a second-order accurate solution was demonstrated in a range of Mach numbers. Moreover, the ability of the method to conserve mass and energy was demonstrated numerically in different scenarios.
Finally, the potential of the proposed methodology to simulate challenging compressible two-phase flows with mass transfer was demonstrated with the three-dimensional nucleate boiling simulation in water. Initially, a vapor film was developed at the bottom heated wall, which in time broke up into a torus-like structure at the center and other smaller structures along the edges of the bottom wall. Eventually the torus structure generated four large bubbles, before releasing their vapor content to the top of the cavity.
Appendix A. Derivation of the relaxed pressure and temperature equilibrium model In the limit of instantaneous velocity, pressure, and temperature equilibrium, the 7equation of Baer-Nunziato model [15] is reduced to a 4-equation model. This relaxation procedure is presented here, following the asymptotic technique used by Murrone and Guillard [22] (see also [95]) to derive the 5-equation Kapila et al. model [20] from the 7equation model. The notation followed is the same as the main body of the manuscript, except from symbols µ, ϕ and τ which are redefined in the context of this Appendix. First, we write the Baer-Nunziato system in the variant of Saurel-Abgrall [13] in terms of the vector of primitive variables w ∈ R 7 specified below as, Here, p I and u I are the interface pressure and velocity. Vector Ψ (w) contains all the relaxation source terms describing transfer processes that we consider in the limit of instantaneous equilibrium: the velocity relaxation term λ(u 2 −u 1 ), the pressure relaxation term µ(p 1 − p 2 ), and the thermal relaxation term θ(T 2 − T 1 ). The relaxation parameters are redefined as λ =λ τ , µ =μ τ and θ =θ τ . Vector Φ(w) contains all the other source terms (e.g. gravity). The focus is on the behavior of the solutions of Eq. (A.1) in the limit τ → 0 + . It is expected that these solutions are close to the set U = {w ∈ R 7 ; Ψ (w) = 0}. Furthermore, it is assumed that the set of equations Ψ (w) = 0 defines a smooth manifold of dimension 4 and that for any w ∈ U, a parameterization Ξ (the Maxwellian) is known from an open subset Ω of R 4 on a neighborhood of w in U. For any v ∈ Ω ⊂ R 4 the Jacobian matrix dΞ v is a full rank matrix, moreover, the column vectors of dΞ v form a basis of ker(Ψ (Ξ(v))) [22]. Based on the above, matrix C ∈ R 7×7 can now be defined as, where dΞ 1 v , . . . , dΞ 4 v are the column vectors of dΞ v and {V 1 , V 2 , V 3 } is a basis of the range of Ψ (Ξ(v)). Based on the observations above, the matrix C is invertible. Another 4 × 7 matrix P can be defined, comprising of the first 4 rows of the inverse C −1 . With the use of matrix P the following results can be obtained (see [22]), P dΞ v = I 4 and P Ψ (Ξ(v)) = 0, (A. 3) where I 4 denotes the 4×4 identity matrix. In order to obtain a reduced velocity, pressure and temperature equilibrium model, solutions in the form w = Ξ(v) + τ z are pursued, where z is a small perturbation around the equilibrium state Ξ(v). Using this form into the system (A.1) one obtains, By multiplying the above equation by P , using (A.3), and neglecting terms of order τ , the following reduced model system is obtained, where v ∈ R 4 and A r (v) ≡ P A(Ξ(v))dΞ v ∈ R 4×4 . In the limit of instantaneous velocity, pressure and temperature relaxation, u 1 = u 2 = u I = u, p 1 = p 2 = p I = p, T 1 = T 2 = T , the vector of the variables of the reduced pressure-relaxed model can be expressed as, The Jacobian dΞ v ∈ R 7×7 of the Maxwellian is expressed as, A basis {V 1 , V 2 , V 3 }, V k ∈ R 7 , k = 1, 2, 3, for the range of Ψ (Ξ(v)) ∈ R 7×7 is found as, (A.9) The matrix C ∈ R 7×7 (A.2) can then be constructed, inverted, and matrix P ∈ R 4×7 can be obtained by taking the first 4 rows of C −1 . Finally, the reduced limit 4-equation model is obtained from (A.5) by calculating the new matrix A r (v) ≡ P A(Ξ(v))dΞ v ∈ R 4×4 and the new source term Φ r (v) = P Φ(Ξ(v)) ∈ R 4 . Note that A(w) and Φ(w) are evaluated in the equilibrium state Ξ(v) in (A.7). Following this procedure, the reduced limit 4equation model takes the form (here assuming Φ = 0), ∂ t a 1 + ∇ · (a 1 u) + S (3) a − a 1 ∇ · u = 0, (A.10a) T − T ∇ · u = 0, (A.10b) ∂ t (ρ u) + ∇ · (ρ u ⊗ u) + ∇p = 0, (A.10c) ∂ t p + u · ∇p + ρc 2 ∇ · u = 0, (A.10d) where, The effects of viscosity, thermal conductivity and mass transfer are described by source terms in the governing equations. The formulation of the four-equation model with this enriched physical description in terms of the conserved variables is given by Eq. (1). Here we show how to derive from (1) the expressions the source terms in the formulation (10). Based on (1), we can write: with a transformation matrix, where φ k and ζ k are defined in Eqs. (13f) and (13g) respectively and, Therefore, the source terms of w (RHS(w)) are determined as, a M + S   Taking into account the expressions for the derivatives of E k , along with E 1 − E 2 = ρ 1 h 1 − ρ 2 h 2 and h k = (c 2 k − χ k )/Γ k , the coefficients of the source terms are obtained as, where, (φρ) T = a 1 φ 1 ρ 2 + a 2 φ 2 ρ 1 , φ v = a 1 φ 1 + a 2 φ 2 , (B.11a) (ζρ) T = a 1 ζ 1 ρ 2 + a 2 ζ 2 ρ 1 , ζ v = a 1 ζ 1 + a 2 ζ 2 , (B.11b) Note that the same results could be obtained by writing the source appropriate term Φ in the Baer-Nunziato equations before the asymptotic procedure described in the previous Appendix.