A weakly Compressible, Diffuse-Interface Model for Two-Phase Flows

We present a novel mathematical model of two-phase interfacial flows. It is based on the Entropically Damped Artificial Compressibility (EDAC) model, coupled with a diffuse-interface (DI) variant of the so-called one-fluid formulation for interface capturing. The proposed EDAC-DI model conserves mass and momentum. We find appropriate values of the model parameters, in particular the numerical interface width, the interface mobility and the speed of sound. The EDAC-DI governing equations are of the mixed parabolic–hyperbolic type. For such models, the local spatial schemes along with an explicit time integration provide a convenient numerical handling together with straightforward and efficient parallelisation of the solution algorithm. The weakly-compressible approach to flow modelling, although computationally advantageous, introduces some difficulties that are not present in the truly incompressible approaches to interfacial flows. These issues are covered in detail. We propose a robust numerical solution methodology which significantly limits spurious deformations of the interface and provides oscillation-free behaviour of the flow fields. The EDAC-DI solver is verified quantitatively in the case of a single, steady water droplet immersed in gas. The pressure jump across the interface is in good agreement with the theoretical prediction. Then, a study of binary droplets coalescence and break-up in two chosen collision regimes is performed. The topological changes are solved correctly without numerical side effects. The computational cost incurred by the stiffness of the governing equations (due to the finite speed of sound and the interface diffusion term) can be overcome by a massively parallel execution of the solver. We achieved an attractively short computation time when our EDAC-DI code is executed on a single, desktop-type Graphics Processing Unit.


Introduction
The usual way of modelling the low speed and incompressible flows, called later on as the truly incompressible approach, is based on the law of momentum conservation and the assumption that the speed of sound c s , in comparison to the convective velocity scale, is high enough to be considered as infinite. The assumption c s = ∞ implies the elliptic character of the system. The pressure field is a solution of the Poisson-type equation and is no longer treated as a thermodynamic quantity but rather as a source of momentum which enforces a kinematic constraint on the velocity field, i.e. ∇ ⋅ = 0 . The advantages of this traditional approach are: (1) only one velocity scale is present and (2) the density is indeed constant. The second feature is especially advantageous in the case of two-phase flows where the densities of particular fluids differ-the density varies only at the fluidfluid interface, while it remains constant in the bulk of each phase. There are, however, disadvantages in the numerical context, both in single and multiphase flows: the algorithmic complexity stemming from the operator splitting and the necessity of performing some iterative subprocesses (due to elliptic nature of the governing equations). It is well known that the truly incompressible models do not take a significant advantage of the parallelisation of the computational process, especially when graphics processing units (GPU) are used as the computing devices. Only fewfold speed-ups (in comparison to a single desktop CPU) are reported in the literature and the main reason is the elliptic type of the equations.
The situation changes significantly if one assumes a finite speed of sound. The governing equations are then of the mixed parabolic-hyperbolic type. They may even become purely parabolic as is the case of the Entropically Damped Artificial Compressibility (EDAC) model for single-phase flows, see Sect. 2. This is advantageous for the algorithmic and numerical handling. Unfortunately, due to disparity of the velocity scales of the convective transport and the pressure waves propagation, which is required to approach the incompressible flow limit, the governing equations are stiff. The resulting time step restrictions can be counterbalanced by efficient parallelisation of the computations. On the other hand, the acoustic (compressible) effects cause some problems when applied to two-phase flows; these issues are comprehensively addressed in this work.
In the present paper we propose a novel mathematical model for the simulation of twophase interfacial flows. The flow is described with the assumption of finite speed of sound by means of the EDAC model introduced by Clausen (2013). In our opinion (as argued in Kajzer and Pozorski 2018a) a better-suited name of the approach would be EDWC where WC stands for "weakly compressible". Concepts similar to EDAC were also proposed in Borok et al. (2007), Ohwada and Asinari (2010) and Toutant (2017). An interesting discussion and comparison of these approaches, called also the general pressure equation methods, has recently been presented in Shi and Lin (2020) and Dupuy et al. (2020).
It has to be noted that there exist alternative, well established WC approaches to the modelling of incompressible flows. The Lattice Boltzmann Method (LBM) (Succi 2001) which solves the mesoscopic level equations was proven to be an accurate and efficient tool for the simulation of single-and multiphase flows, see e.g. Schoenherr et al. (2011), Moqaddam et al. (2016). In the LBM, the Boltzmann equation in a discretised phase space is solved and the discretisation is done using (most often) uniform grids, or lattices, and a finite set of the propagation velocities. It can be shown that the velocity and density fields, computed from the resolved probability density functions, on the macroscopic level satisfy the mass and momentum conservation equations in the weakly compressible regime. One can also solve the macroscopic equations directly using the finite difference or finite volume methods as recently presented in Bigay et al. (2017) and Vittoz et al. (2019). The hyperbolic character of the mass conservation equation requires the upwind-type schemes for the spatial discretisation which can suffer from the accuracy deficiencies and numerical diffusion in the case of low-Mach number flows; however, some improvements in this topic have been proposed (Thornber et al. 2008). Also, when very low Mach numbers need to be imposed, the accuracy could be impaired due to the appearance of short pressure waves. Since the present work is devoted to the twophase flows we also mention the Smoothed Particle Hydrodynamics (SPH) approach (Monaghan 2012). This method is Lagrangian and meshfree; therefore, it is well-suited for the description of the fluid-fluid interfaces and free surfaces, see for example Olejnik and Pozorski (2020). SPH is most often (due to simple implementation) used with the assumption of weak compressibility. On the other hand, SPH is a time-consuming method in comparison to mesh-based approaches and, unfortunately, is only first-order accurate.
For the description of the two-phase systems we use here the so-called one-fluid formulation. A comprehensive classification of the two-phase flow models together with a short description of the one-fluid approaches can be found in Tryggvason et al. (2011) and Mirjalili et al. (2017). The one-fluid formulation needs either explicit tracking of the fluid-fluid interface represented by Lagrangian markers or implicit interface capturing using an additional field which is governed by model-specific equations. The class of interface capturing approaches contains the Level-Set (LS) methods, the Volume of Fluid (VoF) method and the Phase Field Models (PFM). The relationship between LS, VoF and PFM was discussed by Wacławczyk (2017). In our opinion, the main feature that distinguishes PFM from LS and VoF is that the PFM approach takes into account, at least approximately, the physics underlying the interface dynamics (including the topological changes) while LS and VoF are based on purely geometrical considerations.
In this work we consider the diffuse-interface (DI) model that originated in the work of Sun and Beckermann (2007) and was rendered mass-conservative by Chiu and Lin (2011). The latter model variant (after a re-arrangement of the coefficients) is also referred to in the recent literature as the conservative Allen-Cahn (CAC) equation, see e.g. Fakhari et al. (2017), and indeed can be considered as a member of the PFM family.
For the first time the EDAC model with DI technique was proposed in a basic variant, along with rudimentary numerics, in our earlier work (Kajzer and Pozorski 2018b). Here, a considerably generalised model is presented together with a comprehensive account on the discretisation schemes that are developed specifically for the EDAC-DI equations. Recently, a WC approach coupled with the original model of Chiu and Lin was reported by Matsushita and Aoki (2019). In their approach, the flow solver was based on the method of characteristics and the DI model was complemented with the standard LS equation (its solution was used to improve the surface tension computation). This makes the numerical methodology quite sophisticated as compared to our EDAC-DI model which is relatively simpler to solve. This paper is organised as follows. We first put forward the EDAC-DI model and find the expressions for model coefficients. Next, we present the proper numerical methods for solving the governing equations. Then, the results of three-dimensional simulations of steady droplet and binary droplets collisions, including the coalescence and break-up regimes, are presented. Finally, we briefly discuss the computational efficiency of the flow solver implemented for the execution on GPU. As we only show a few examples of twophase flows, the present work should be perceived as a "proof of concept" of the EDAC-DI approach.

The EDAC Flow Model
Due to its favourable features (Kajzer and Pozorski 2018a), we use the EDAC model as the flow solver. For its detailed derivation for single phase flows the reader is referred to Clausen (2013). Here, for the sake of brevity, we only recall the EDAC equations as needed for further considerations: where p is the pressure, is the velocity, 0 and 0 are the (constant) fluid density and kinematic viscosity, respectively, and is the vector of body forces. Equation (1) is the usual momentum conservation equation, i.e. the Navier-Stokes equation of the single phase flow, and Eq. (2) establishes the EDAC model. It is derived from the entropy balance equation (hence the phrase "entropically damped" in the model name). One assumes that the speed of sound c s ≫ | | max , i.e. the Mach numbers Ma = | | max ∕c s are low. Nevertheless, the incompressibility constraint ∇ ⋅ = 0 no longer holds. Therefore, despite the very acronym of EDAC, its affiliation to the family of artificial compressibility (AC) methods is in our opinion misleading. In the classical AC approach of Chorin (1997) the concept has been used for steady flows and when the computations converge, the resulting velocity field is divergence-free. For unsteady flows one introduces an artificial time and iterates the discretised momentum and pressure equations at each physical time step till ∇ ⋅ = 0 is satisfied. Clearly, this is not the case of EDAC approach. On the other hand, the density is (artificially) assumed constant which is physically consistent (asymptotically) at Ma → 0.
Clausen's main idea consists in the introduction of a physically justified (derived from the temperature diffusion) term 0 ∇ 2 p . Therefore, the EDAC model is of purely parabolic type. This is important from the numerical point of view since centred spatial discretisation schemes without addition of numerical or artificial diffusion can be used if the grid is sufficiently fine. Indeed, this was proven by Kajzer and Pozorski (2018a). The pressure diffusive term significantly reduces the noise in the velocity divergence field, bringing the flow closer to the incompressible limit. As for the implementations of EDAC to date, apart from some laminar and turbulent flow cases reported in the original paper of Clausen (2013) and the direct numerical simulation of wall-bounded turbulence (Kajzer and Pozorski 2018a), the model has recently been studied in the context of Large Eddy Simulation (Delorme et al. 2017).
To be able to handle two-phase flows in the weakly-compressible approximation, we propose a more general form of the EDAC model, permitting variable density and viscosity (due to material properties of the phases): These equations will be used in the following to obtain the final one-fluid formulation of the two-phase interfacial flow model. Again, the only difference in comparison to the truly incompressible approach is that the divergence-free condition, ∇ ⋅ = 0 , is replaced by the EDAC pressure equation, Eq. (5). Importantly, the eigenstructure of the system remains unchanged, i.e. information propagates with the maximum speed of c s + | | max .
As it was mentioned in the Introduction, one can argue whether the EDAC model belongs to AC or WC methods in the case of single-phase flow since the density variations are then neglected. In Eqs.
(3)-(5), the mass conservation is also explicitly taken into account, admitting local density variations that stem both from the two-phase nature of the system and the compressibility effects. Therefore, it is a true member of the WC family. Importantly, the pressure is explicitly evolved in the considered model so there is no need of any equations of state.

The Diffuse-Interface Model
The PFM approach is often associated with the Cahn-Hilliard (CH) equation, see e.g. Magaletti et al. (2013). The CH model has a strong physical foundation: it expresses a conservation law and does not involve any geometrical information, like direction normal to the interface. On the other hand, the CH equation is a nonlinear 4th-order partial differential equation (PDE), so it is difficult to solve numerically. Explicit schemes suffer from very short timesteps dictated by the high order spatial derivative of the so-called order parameter and therefore implicit schemes have to used. Moreover, issues according to the spurious shift of the solution arise and spurious shrinkage of the droplets or bubbles occurs. An alternative, easier in numerical handling, is the Allen-Cahn equation, which is a 2nd-order PDE. The serious drawback of that model is the lack of mass conservation. Therefore we have decided to put forward another approach. The very starting point is the PFM-type equation proposed by Sun and Beckermann (2007). This equation involves only 1st and 2nd spatial derivatives but requires the vectors normal to the interface. However, the model does not conserve mass. This deficiency was remedied by Chiu and Lin (2011) who introduced the conservative PFM given by the following 2nd-order PDE: In the above equation ∶ d ↦ [0, 1] is the order parameter (which we will also call the phase indicator function); it takes values 0 and 1 in the regions occupied by the separate fluids. The interface is understood here as a finite width transitional band where 0 < < 1 ; its actual position is determined by the isoline (or isosurface in 3D) = 1∕2 . The direction normal to the interface is defined by a unit vector = ∇ ∕|∇ | . If we neglect the advection term the dynamics is similar to the Burgers equation: the compressive term (1 − ) tends to create a "shock" and the diffusive term smoothes it. The steady state solution of Eq. (6) in one spatial dimension is given by Chiu and Lin (2011): where is the local coordinate in the direction perpendicular to the interface located at ̃ . The constant has the dimension of length and controls the width of the interfacial region. Obviously, the real world values of the order of 10 −9 m can not be used in any feasible numerical solution of macroscopic flow problems. Therefore, one sets the value of proportionally to the spatial resolution x . The width of the interface has a significant influence on the stability and accuracy of the solution since it controls the smoothness of the density and pressure gradients across the interface. On the other hand, setting a large value of reduces the effective resolution of the simulation. The second model parameter in Eq. (6) is . It obviously has the dimension of velocity and controls (in the absence of advection) the rate of approaching a steady state. From the literature it is, however, not clear what are the proper values of from the physical point of view. Let us rewrite Eq. (6) in the form which is also referred to as the conservative Allen-Cahn equation (Fakhari et al. 2017): where = 4 and ̂= ∕4 . It is now clear that ̂ is a diffusion coefficient (it has dimension of m 2 ∕s ), sometimes referred to as the mobility. This is somehow improper since the dimension of mobility is, in fact, m 3 s∕kg (Magaletti et al. 2013). To retrieve the dimensional consistency let us recall the Cahn-Hilliard equation (Magaletti et al. 2013): where the physical quantities are the fluid mobility and the surface tension coefficient . Clearly, the coefficient 3 Let us now focus on the choice of the mobility value which is not straightforward. Physically-sound values of the mobility are of the order of 10 −17 m 3 s∕kg . Although they reflect the microscopic nature of the interfacial phenomena, such low values can not be handled numerically in an otherwise macroscopic model and, alike for , higher values have to be set. Importantly, there exist a value of mobility which is optimal when the CH equation is used for the two-phase flow simulation and we will follow the work of Magaletti et al. (2013). Although their analysis was performed for the CH equation, Eq. (9), one can expect that the result gives at least reasonable approximation of a proper mobility value in the CAC equation since it reveals the same dynamics. The dimensionless mobility is * = ∕(UL 2 ) , where U and L are characteristic velocity and length scales. Magaletti et al. (2013) recommended * = 3Cn 2 with Cn = ∕L being the Cahn number. Therefore, we get: Using the above relation in Eq. (10) and recalling that = 4 , we get: This way, we obtain a physically-sound and no longer adjustable estimation of the parameter in Eq. (6) as = 12U . It should be emphasised that this value is much higher than in the original work of Chiu and Lin and in the recent paper of Mirjalili et al. (2019) who set = U . The probable reason is that the DI model considered there was coupled to truly incompressible flow solvers and high values of led to prohibitively short time steps when explicit time integration was used. Moreover, high values of quickly lead to severe unphysical deformation of the interface if improper spatial discretisation is applied, as we will show later. To obtain our final model we will use Eq. (12) and for convenience we will drop the overbar from the symbol .
Let us make an additional comment on the original conservative level set (CLS) method of Olsson and Kreiss (2005) since it is quite popular and, in a sense, related to the DI model we consider here. The CLS equations are: where is the artificial time of the so-called re-initialisation stage, Eq. (14). This model can be considered as a splitting of the advection and DI operators in Eq. (6) by setting = t and taking one re-initialisation step. Although the normal vectors used in the CLS are fixed during the re-initialisation process ( 0 denotes the normal vectors obtained from distribution at the most recent physical time step), the number of steps in the artificial time with given acts in a similar way like the mobility: the long-time integration of the re-initialisation equation corresponds to high mobility value. Therefore, insufficient re-initialisation time after each physical time instant can lead to non-physical results and requirement of special numerical treatment, as exemplified by artificial rupture of the gas film in the simulation of droplets collisions, cf. Amani et al. (2019). Olsson et al. (2007) solved this issue by proposing another CLS model in which the diffusive term is projected on the direction normal to the interface, resulting in a formulation based on purely geometrical considerations.
Another remark should be made here: using an iterative subprocess in artifical time between subsequent physical time steps (the level-set methods, also some algebraic VoF approaches, see e.g. So et al. 2011) rises questions about the accuracy since the distribution of the density is modified while the distribution of momentum remains unchanged.

The EDAC-DI Model
The main idea of the one-fluid formulations is that the local fluid properties are computed (directly or after some mollification as in sharp interface methods, e.g. VoF) from the phase indicator function: where 0 , 1 and 0 , 1 are the densities and viscosities of the phases, respectively. Without the loss of generality we can assume 1 > 0 . For brevity, we will also use the notation [ ] = 1 − 0 for the jump of the density across the interface. Provided that the densities differ, we can write the trivial inverse relation of and to be used next:

3
Using the above relation, from Eq. (6) we can obtain an "interfacial" continuity equation: where = ∇ ∕|∇ | . Obviously, considering Eq. (7), the stationary solution is: The reasons to solve the equation for rather than the one for are threefold: (1) one avoids the multiplication of possible dispersive numerical errors by an arbitrarily large factor [ ] ; (2) there are less arithmetic operations to be performed in the solution algorithm; (3) the interpretation of the governing equations is straightforward. It should be noted that the formulation (17) is less general than solving Eqs. (6) and (15) since it does not permit the fluids density ratio 0 ∕ 1 = 1 which, on the other hand, is less common in practice. We emphasise that the actual interface is now defined in a natural way as the isoline (isosurface) corresponding to = ⟨ ⟩ with ⟨ ⟩ = ( 1 + 0 )∕2 being the mean of the material densities. Apart from the interfacial region, Eq. (17) should reduce to the usual mass conservation law. However, due to the weak-compressibility effects it is possible that ∇ ≠ also apart from the interface. The DI terms and the surface tension acting in the "spurious" interfacial region (i.e. locations where ∇ ≠ is caused by compressibility) could further non-physically steepen the density gradient (or, in case of low surface tension, excessively smear the compressible density variations). Also, for the sake of physical consistency, the EDAC pressure diffusion should not act in the interfacial region: there is no reason to artificially smooth the pressure jump across the interface. Moreover, if the pressure diffusion is not switched off on the interface, a constant shift of mean pressure appears, which is fed by the surface tension. A question arises whether applying the EDAC pressure diffusion is justified at all. The answer is positive since in the well-resolved flow regions it will reduce the action of the limiters used in the numerical approximation.
Summarising, one needs to distinguish between the interfacial region and the density fluctuations in the bulk of the phases stemming from compressibility. Therefore, we propose the following EDAC-DI model, cf. Eqs. (17), (4) and (5), as the final formulation: where st is the surface tension force (we do not take the body forces into account for brevity) and is a switch that takes the value of 1 in the interfacial region and 0 apart from it. The density variations due to compressibility are proportional to Ma 2 . The same holds for the phase indicator function, i.e.
∼ Ma 2 and variations of such magnitude appear where ≈ 0 and ≈ 1 . Consequently, unlike in Eq. (6), the constraint 0 ≤ ≤ 1 no longer holds in the WC implementation. Before we explicitly define the switch , let us introduce auxiliary, clipped fields and ̃: The switch is defined as follows: where is a threshold to determine whether a given location is attributed to the interfacial region. Keeping in mind the discussion following Eqs. (19)- (21), we set = Ma 2 . The coefficient has to be chosen carefully. If it is too low the surface tension and the DI terms can possibly act also in the bulk of the particular phases. On the other hand, too high value of will artificially narrow the region where the surface tension and DI terms should act, leading to oscillations of the pressure field and excessive diffusion of the density. We have found that = 1 was a proper and rather safe choice in all the cases considered here; however, one has to remember that it is dependent on the numerical method used (in particular the numerical diffusion). The switching procedure can be summarised as follows: (1) get rid of the compressible overshoots (i.e. regions of mass accumulation) where ∼ 1 and the undershoots (mass rarefactions) where ∼ 0 ; (2) remove (by setting ̃ ) the undershoots where ∼ 1 and the overshoots where ∼ 0 ; (3) mark the interface as the locations where <̃< 1 − . Notice that this procedure is performed explicitly and cell-wise so neither the auxiliary fields nor the switch have to be stored.
It is clear that in the limit Ma → 0 , when compressible density variations vanish, no switch is needed and the DI terms can be applied in the whole domain keeping the physical consistency. This would result in the parabolic type of the mass conservation equation in the whole domain allowing us to use centred spatial schemes (without any upwinding or artificial diffusion). The dispersive errors would be damped by the DI diffusive term since the cell Peclet numbers Pe = U x∕(12U ) ≪ 1 , as ∼ x . On the other hand, in conjunction with the truly incompressible solvers whose efficiency is based on the absence of sound waves, an implicit time integration has to be applied to the DI equation, Eq. (19), due to the stiff diffusive term, see Sect. 4 for details. The divergence-free velocity field is obviously not our case and a more sophisticated numerical technique than centred schemes for the spatial discretisation has to be used (Sect. 3.1).
The surface tension force is defined by Tryggvason et al. (2011): where is the identity tensor, S is the surface Dirac delta function related to the interface and is the surface tension coefficient (assumed constant in this work). The local dynamic viscosity, in order to avoid negative values, is computed according to Eq. (15) using ̃ as the order parameter; the same is done when computing the EDAC pressure diffusion coefficient. We emphasize that the clipped values ̃ are used only for the evaluation of local viscosity coefficients; the density value taken in the computation of other terms is unchanged, so the mass and momentum are conserved. The normal vectors used for the computation of the DI compressive term and the surface tension could also be determined using ̃ , but since the overall behaviour of the solution is sensitive to the numerical errors in the computation of gradients, the use of clipped field ̃ should be avoided ( or are to be used).
In the flow cases that we consider in this work the maximal velocity of the flow varies significantly during the simulation (e.g., due to the capillary waves created by the topological changes of the interface). Choosing the reference velocity scale U is therefore problematic. Setting U = | | max based on theoretical considerations before the simulation can lead to, e.g., overestimation of the speed of sound and the mobility. Therefore, we adaptively set the reference velocity as the maximum speed in the domain with 10% safety factor, i.e. U(t) = 1.1 × | (t)| max ; the speed of sound and the mobility are also set adaptively using U(t). This allows us to keep Ma almost fixed during the simulation and assure a proper value of without significant overestimation. This is important since very low Mach numbers should be avoided due to excessive numerical diffusion and creation of short pressure waves of high magnitude; the role of was discussed earlier. Additionally, this allows one to achieve higher computational efficiency since the time step is also set adaptively.

Summary of the Model Parameters
The model parameters that have to be set are: the Mach number Ma, the numerical interface width , and the coefficient . The choice of was already discussed and we will always use = 12U . During the investigation we found that = 0.75 x is a reasonable choice with respect to the accuracy and resolving power; this results in the interface being resolved using ∼ 8 grid cells. The choice of the Mach number is mainly dictated by the computation of the surface tension: to avoid significant pressure oscillations in the vicinity of the interface one has to take into account the locations where |∇ | > 10 −3 ∕ x . Setting Ma=0.03 we obtain the switch = 1 at the positions where the phase indicator variations satisfy > 9 × 10 −4 .

Numerical Method
To numerically solve Eqs. (19)- (21) we use the method of lines: the equations are rewritten in the semi-discrete form with a prior discretisation of spatial operators, resulting in a system of ordinary differential equations. Let us present the spatial discretisation first.

Spatial Discretisation
For the single-phase flow the EDAC equations are of purely parabolic type and, therefore, one can use centred spatial discretisation without introducing artificial or numerical viscosity if the grid is sufficiently fine, as it was proven by Kajzer and Pozorski (2018a). The case of the EDAC-DI model is different: for physical consistency we switch off the DI terms (including the diffusion) of the mass conservation equation in the bulk of the phases and the EDAC pressure diffusion is switched off in the interfacial region. This results in a locally hyperbolic character of the equations and, in the absence of diffusion, the centred schemes would lead to unphysical oscillations in the solution. Moreover, later in this work we consider the simulation of colliding liquid droplets immersed in gas. Even though the droplet diameters are very small (less than 1 mm ) the resulting Reynolds numbers are of the order of 10 3 . The situation gets worse at time instants just before the collision, when velocity rises locally in the region of gas film significantly above the initial velocity of droplets and the film has to be resolved using only few grid cells. We were not able to perform a stable simulation in this case using the central schemes for the discretisation of the hyperbolic part of the governing equations. Therefore, we have chosen to use the scheme proposed by Kurganov and Tadmor (KT) (2000) which is a simple, Riemann solver-free, upwind-type finite volume scheme for convection dominated problems. Below we present the details. The discretisation method is shown here for the case of two spatial dimensions (2D). The extension to three dimensions (3D) is straightforward, except for some terms which are also specified in 3D case. We will use the uniform, cartesian grids with cell edge length denoted by x . All variables are arranged in collocated manner and are stored at the cell centers. Let = (u, v) and = (n x , n y ) . For clarity we rewrite the governing equations as follows: Above, we already replaced the surface tension, Eq. (25), by the Continuum Surface Stress (CSS) model (Tryggvason et al. 2011) and we use s ≈ |∇ |∕[ ] . The CSS model advantage is that it conserves the total momentum which is important in the flow cases dominated by inertia. The main drawback of CSS (with respect to the popular CSF approach) is the appearance of stronger numerical side effect, the so-called parasitic currents, see Tryggvason et al. (2011). We analyse this phenomenon in Sect. 4.
We emphasise that the density gradients and normal vectors used for the computation of the surface tension are marked by a hat symbol. They need to be distinguished from the ones used in the DI equation (without hat symbol) since they are discretised in different ways. In the following we will use the approximation of local kinematic viscosity at the cell faces given by: The local dynamic viscosity and the switch are computed in the same way.
The viscous flux in the momentum equation is discretised using centred differences and for the off-diagonal terms the arithmetic mean of cell-centred values is used: The density and pressure diffusive terms are discretised by analogy with Eq. (30).
The CSS term in Eq. (27) is discretised using the following approximations: 1 3 and The normal vectors are computed as follows: Namely, the density gradients on the cell faces are computed as the average of the cell centred values while the gradient norms are computed as the average value of norms and not the norm of the average gradient; please note that � |∇ | ≠ |∇ | . We have found that this discretisation of the CSS model allows to avoid spurious topological changes of the interface during coalescence process and results in the lowest magnitude of the parasitic currents among the other combinations of schemes. This observation corresponds, to some extent, to the results of Jamet et al. (2002). The KT scheme (Kurganov and Tadmor 2000) for the mass and momentum convective fluxes reads, respectively (we present only the x-direction): where a i+ 1 2 ,j = c s + max(|u L i+ 1 2 ,j |, |u R i+ 1 2 ,j |) and the reconstructions of the left and right states of a generic variable f (here, , u or p) are done by the MUSCL procedure: where x f i,j is the x-component of the gradient of f computed using a chosen limiting procedure. The above discretisation of the momentum convective term is, however, improper in our application: the density gradient in the interfacial region, that stems from the material properties and not from the flow, should not add diffusion to the momentum flux. Therefore, using the elementary identity x (fg) = f x g + g x f and the mean value of the left and right states, we can fix this problem by setting: while keeping consistency at the locations in the bulk of the phases. The pressure flux is also computed using the KT scheme with the MUSCL reconstruction. Let us now present the discretisation of the EDAC pressure equation. The required left and right states of the velocity components and the pressure are already computed for the discretisation of the momentum fluxes. The acoustic and convective terms read: It should be emphasised here that we do not apply the switch for the numerical diffusion, although we turn off the EDAC pressure diffusion term in the interfacial region. This is required to control the pressure oscillations during the topological changes of the interface. We now discuss the discretisation of the DI terms of Eq. (26). The discretisation of the momentum and mass fluxes consists in approximation of the flow fields (ρ, u and p) at the cell faces and then the numerical fluxes are computed using these face values. Such a procedure can be also applied using centred interpolation of variables, leading to the so-called skew-symmetric or split forms of the nonlinear fluxes, see, e.g., Pirozzoli (2010), Kennedy and Gruber (2008). This approach, however, should not be applied to the DI compressive term. Instead, we should use the so-called divergence form: compute the compressive fluxes at the cell centres and then average them on the cell faces. Namely, setting and we compute while the following split form should be avoided: When using the above form, the normal vector at the cell face could be computed using a multidimensional stencil, like it is done for the computation of the off-diagonal terms of (38) the viscous stress tensor, or by averaging the neighbouring cell centred values. Unfortunately, independently of the way the normal vectors are computed, the split forms lead to a spurious deformation with tendency to "squaring"-a circular interface becomes a square. It was found that computation of the normal vectors at the cell centres also needs special attention as far as shape preservation of the advected structures is considered. The first option is to discretise the density gradients in the dimension-by-dimension (DBD) manner using the discrete operators introduced above: However, this leads to a spurious deformation of the interface, see Fig. 1. These deformations are more significant with increasing and decreasing . We found that in 2D cases it is much better to use a genuinely multidimensional (MDIM) stencil and compute the The non-isotropic behaviour of the numerical solutions of PDEs with the DBD spatial discretisation strategy was described by Shukla and Giri (2014) who recommended the use of genuinely MDIM schemes with isotropic truncation errors. Here, we also recall the proper formula for the discretisation of the normal vectors in 3D since it is not a straightforward analogy to the 2D case: where, e.g., and the interpolation on the cell face is done using 2nd-order averaging formula. Notice that, in practice, since we are interested in the cell centred values, the above formulation is simplified to a form similar to 2D case. For example, to compute the derivative in the x-direction one will not use cells marked by i at all-only i ± 1 will be involved. Morover, the 3D formula is, fortunately, simpler than its straightforward extension of 2D case and does not require the full 27-cell-stencil nor the use of the 2D Simpson rule on the cell faces. In the present paper we apply these formulae to the computation of the DI normal vectors and the gradients required by the Multidimensional Limiting Process (MLP) (Hubbard 1999) that we use for the flux approximation. The MLP procedure in 2D is described in "Appendix". The necessity of using the multidimensional reconstruction will be proven in Sect. 3.3.

Temporal Discretisation
The time advancement of the solution is done by means of the 2nd-order, strong stability preserving Runge-Kutta (SSPRK2) method, see e.g. Gottlieb and Shu (1998). Our choice is dictated by the rich variety of physical phenomena governed by the EDAC-DI equations and we look for overall robustness of the numerical methodology. An additional advantage is that this method requires only two registers per unknown variable which is important in GPU computing (due to relatively low global memory available). The stability of an explicit time integration schemes for EDAC-DI equations in d-spatial dimensions is ascertained when: where generalised diffusion coefficient = max( 0 , 1 ) (stemming from the momentum and pressure equations) or = (stemming from the DI equation). The maximum allowable time steps also depend on the applied type of the spatial discretisation. When the surface tension is taken into account, another time step restriction established by Brackbill et al. (1997) is: It stems from the CFL condition based on the estimated speed of capillary waves, c 2 = k∕( 0 + 1 ) , where k is the highest resolved wave-number; for the 2nd-order, threepoint central scheme we have k = ∕ x . However, as pointed earlier, we set the speed of sound and adaptively by tracking the maximal velocity of the flow which also takes the capillary waves into account. Therefore, the restriction on the time step given by Eq. (50) can be neglected. An interesting discussion on the time step constraints stemming from the surface tension is presented by Denner and van Wachem (2015).
Let us compare the relevant time step restrictions. In this work we consider low viscosities, therefore we apply the CFL condition and the DI diffusion related time step limits: and respectively, where c = ∕ x is the proportionality coefficient. The ratio of these time steps is: For the parameters used, we have r ≈ 1.05Cr∕C diff . This means that one should apply the time integration methods for which the maximum allowable values of Cr and C diff are in relation Cr < C diff . This is the case of most of the Runge-Kutta methods while, e.g. the Adams-Bashforth methods should be avoided. On the other hand, when using large c , the diffusive time step constraint will be always more severe than CFL. In such case one could consider the use of, e.g., the Runge-Kutta-Chebyshev methods which have a large stability region along the real axis.
When the advection-diffusion equation is considered and the maximal t CFL and t DI are of similar magnitude, as in our case, the maximal allowable Cr is a function of the Peclet number, Cr = Cr(Pe ) . In the presence of advection, the maximal allowable C diff is also lower than in the pure diffusion case. The time integration methods with Cr less sensitive to the diffusion coefficient are proposed by Torrilhon and Jeltsch (2007), however the (51) improvement is done at the cost of additional storage. We have found that, in the considered range of the model parameters, stable computations with SSPRK2 are possible when Cr ≤ 0.5 and C diff ≤ 1.2 , resulting in r < 1∕2 so it is still the speed of sound that limits the time step.

Assesment of the Numerical Method for DI Equation
First, we analyse the influence of the chosen discrete forms of the compressive part of the DI flux and normal vectors approximations on the shape preservation. The test consists in solving Eq. (12) without advection in a periodic square domain [0, 1] × [0, 1] . The initial condition is the circular interface of diameter D = 0.5 , set by adapting Eq. (7) to two dimensions. The resolution is set that D∕ x = 50 . The time step was set according to the remarks in the previous subsection assuming U = √ 2 and the computations are performed till t = 10 . The results are shown in Fig. 1. It is clear that the split form, Eq. (43), gives very bad results independently of the way we compute the normal vectors, while the divergence form, Eq. (42), performs much better. Importantly, if an improper discretisation is used, it is observed that the spurious deformations become stronger with decreasing and increasing (not shown). On the other hand, using large leads to a loss of the effective spatial resolution and should be avoided. Fortunately, the divergence form with multidimensional scheme for the discretisation of the normal vectors offers enormous improvement over the other forms and allows us to use arbitrarily high value of without significant non-physical deformation of the interface.
To investigate the proper scheme for the advection term we tested, along with the already mentioned MLP reconstruction, three popular one-dimensional limited slope reconstructions: superbee (SB), monotonised-central (MC) with limiting constant 1.3, and the optimised MUSCL (OM) reconstruction of Leng et al. (2012). We consider two cases: advection in the direction aligned to grid lines and in the direction 45 • oblique to the grid lines. In the first case U = 1 and in the second case U = √ 2 . To mimic the situation of coupling the DI equation with EDAC equations we take the speed of sound stemming from Ma=0.03 in the diffusive (upwind) part of the KT fluxes. The results for the one-dimensional limiters are shown in Figs. 2 and 3. Clearly, the SB and MC limiters perform very poorly: even after a short time a severe deformation is observed. The OM scheme performs better but still introduces a significant distortion. The MLP scheme outperforms all the other schemes, however a slight deformation is visible in the case of the grid-aligned velocity field. We used the least compressive variant of MLP with = 1 , see "Appendix". One can eliminate the spurious deformation by making the limiter a bit compressive by setting = 1.1 but, on the other hand, we noticed stability issues when this was applied to the EDAC equations. Therefore, in the following we will always use = 1 . We note that similar results are obtained in the truly incompressible formulation; obviously, some tuning of the limiting parameter is then required.
To obtain meaningful quantitative information we simulated the so-called Rider-Kothe benchmark case (Rider and Kothe 1998). The case is well known and we do not include the very details here, see e.g. Olsson and Kreiss (2005). The domain is again a unit square. The off-centre circular interface of diameter D = 0.3 is rotated (in counterclockwise direction) in a shearing velocity field. The velocity field is reversed after t = T = 1 and the solution should revert back to the initial state at t = 2T . The test was performed at four spatial resolutions using the grid built of N = 80 , 160, 320 and 640 cells in each direction, resulting in D∕ x = 24, 48, 96 and 192 and Cn = 1/32, 1/64, 1/128 and 1/256, respectively. The results are shown 1 3 in Fig. 4. It is seen that to prevent the spurious topological change the Cahn numbers should be ∼ 0.01 or less.
In Table 1 we report the outcome of quantitative analysis. We consider the standard L 1 error of the solution and the L A error which measures the area and shape deviation of the contour corresponding to = 0.5 , see Olsson and Kreiss (2005): where H is the step function and exact is the analytical solution which is also imposed as the initial condition (the errors are computed at t = 2T).
For the low resolutions the order of convergence is higher than 2 since the solutions differ qualitatively. When the resolution is sufficient to prevent the spurious topological change the order of convergence is lower than 2 which is expected since we keep the interface width constant (with respect to the grid spacing x). Fig. 2 The results of translation test using different limiters for the advection term (the abbreviations are expanded in the text) in the case of the grid-aligned velocity field. The line styles and colors are coded in the same way as in Fig. 1 Summarising, the considered DI model with the numerical discretisation detailed in Sect. 3 is accurate in terms of shape and volume preservation but requires relatively fine grids. Avoiding the spurious interface deformations is critical when the surface tension plays a significant role in the flow.

Assesment of the Numerical Method for Variable Density EDAC Equations
The numerical methodology proposed in this work differs significantly from the one we used previously for the DNS of the channel flow using EDAC model (Kajzer and Pozorski 2018a). One can expect significant numerical diffusion from the discretisation presented in Sect. 3.1. For the comparison purposes we took the well known case of Taylor-Green vortex (TGV) at Re = 1600 . As the reference we take the DNS data (Jammy et al. 2016) obtained by solving the compressible Navier-Stokes equations with Ma = 0.1 , on a grid built of N = 512 cells in each direction. We solved Eqs. (3)-(5) using the numerical method Fig. 3 The results of translation test using different limiters for the advection term (the abbreviations are expanded in the text) in the case of the grid-oblique velocity field. The line styles and colors are coded in the same way as in Fig. 1  Fig. 4 The Rider-Kothe test: results obtained by solving Eq. (12). The contours corresponding to = 0.05, 0.5 and 0.95 are shown. The black dashed lines denote the DI solution at t = T and t = 2T , the gray solid line is the reference solution where for t = 2T we used the initial condition, and at t = T (only the contour = 1∕2 is shown) we used the solution of the pure advection equation obtained on 1000 × 1000 grid with the 5th-order monotonicity preserving (MP5) scheme (Suresh and Huynh 1997)  presented above. In Fig. 5 we compare the time evolution of the kinetic energy obtained at two spatial resolutions, N = 192 and 384, at Ma = 0.1 and Ma = 0.03 . The speed of sound was set adaptively during the simulation according to the remarks in Sect. 2.3 and the update was done every 100 timesteps. Clearly, the numerical diffusion is significant and depends on the Mach number. Using N = 384 we did not achieve satisfactory agreement with the DNS reference data computed with N = 256 (not shown). It should be emphasized that the DNS of interfacial flows should assure spatial resolution high enough to resolve the interface on a length scale that is equal or lower than the Kolmogorov scale l K . In our case this would result in l K > 8 x leading to tremendous computational effort but, on the other hand, the turbulence would be then sufficiently resolved, even by a strongly diffusive scheme. However, it seems that one should consider the use of adaptive mesh refinement (AMR) and hybrid spatial discretisation (by blending low and high order schemes) to enhance the resolving power since high numerical diffusion is not needed except in the interfacial region.

Results
We have chosen two benchmark cases (Sects. 4.1 and 4.3) to validate the new model and numerics in their basic apects. We study the coalescence and break-up of liquid droplets immersed in gas in two collision regimes. The proper behaviour of the system in those cases is not straightforward to obtain in many simulation strategies, see e.g. Moqaddam et al. (2016) and Amani et al. (2019). It should be emphasized that both coalescence and break-up are always under-resolved, independently of the model applied. The physically-sound thickness Fig. 5 The evolution of the kinetic energy in TGV case of the gas film between the droplets just before collision and thin filaments created before the break-up are significantly smaller than available grid spacing. Moreover, the maximal resolved speed of capillary waves is limited by the grid spacing, see Sect. 3.
Additionally (Sect. 4.2), we conducted the simulation of a single steady droplet to validate the model against the Young-Laplace law and to estimate the numerical errors emanating as the so-called parasitic currents. In all the simulations the maximal velocity was updated every 100 time steps, which corresponds to physical time interval of the order 10 −5 s.

Head-on Collision of Water Droplets in Air
We first simulate the head-on collision of two equal diameter ( D = 0.5 mm ) water droplets immersed in air. The computational domain is a periodic box with the edge length equal 3D. Initially, the droplets' centres are separated by 1.5D. We have chosen the Weber number (based on the initial droplets velocity U 0 ) We = 1 D(2U 0 ) 2 ∕ = 10 . In this case no secondary break-up should occur upon coalescence (Ashgriz and Poo 1990). Additionally, at this relatively low Weber number, the thin liquid lamella will not appear and therefore the simulation is sufficiently resolved (in the context of interfacial structures) for the whole analysed time period.
The densities of water and air are taken as 1 = 10 3 kg∕m 3 and 0 = 1.226 kg∕m 3 , respectively. The viscosities are 1 = 1.137 × 10 −3 kg∕(m s) and 0 = 1.78 × 10 −5 kg∕(m s) . The surface tension coefficient is = 0.0728 N∕m . The droplets initial velocity (aligned with the x axis) is U 0 = 0.60 m∕s . The resulting Reynolds number is Re = 1 D(2U 0 )∕ 1 = 530 . The simulation is conducted till t = 1.5 ms. We set the spatial resolution to 192 3 , so the cell size is x = 7.8125 × 10 −6 m . Initially the droplets are resolved using 64 grid cells per diameter. The grid spacing results in the maximal capillary waves speed c = 5.4 m∕s.
In Fig. 6 we present an overview of the simulation. After the collision and merging, the resulting droplet expands until the surface tension balances the inertia. We did not observe any spurious break-up at the maximum expansion stage (see the right panel of Fig. 10). Then, the "cylindrical" droplet collapses and expands along the initial velocity direction. Importantly, the evolution of the interface corresponds well to the experimental results of Ashgriz and Poo (1990) (although our simulation is done at lower Weber number).
In Fig. 7 we present the interface shape at t = 0.608 ms. Clearly, the axial symmetry is very well preserved, even in the regions of high curvature. This confirms that the spatial discretisation schemes work correctly, cf. Figs. 1 and 2.
The very moment when the droplets coalesce is shown in Figs. 8 and 9. Most importantly, the density distribution is physically correct, i.e. no gas is trapped inside during the collision. To see the effect of the choice of , as compared to the recommended value of 12U (see Sect. 2.2), we also tested = 6U which indeed results in improper behaviour: tiny air bubbles become entrapped in the resulting droplet and a region of spurious rarefaction of liquid phase appears (not shown). Although the spatial discretisation we use is highly diffusive, some pressure oscillations are still present in the narrow region between the droplets just before collision. On the other hand, they quickly disappear after the coalescence is completed. The switch acts properly: the surface tension and the DI terms are active only in the interfacial zone. The absolute value of the divergence of the velocity field is high in the gas phase at the locations neighbouring the interface. This is an unwanted effect of applying which sharply switches off the surface tension. On the positive side, the variations of ∇ ⋅ are only local. A larger weakly-compressible structure is visible in the under-resolved gas film escaping from between the droplets (see Fig. 9: third row, middle column). The values of the velocity divergence are quite high but one has to remember that the variations of the density due to the compressibility effects scale as ∼ t ∇ ⋅ and in our simulations t ∼ 10 −8 s which results in ∕ ∼ 1% as shown later.
In Fig. 10 the pressure field cross-sections are shown at later time instants when the interface is strongly curved. The pressure is smooth and, as expected, it achieves extremal values near the interface where the curvature is the highest.
Finally, in Fig. 11 we report some quantitative results. In the left panel the history of the maximal velocity in the flow is presented. Three peaks are visible: the first one corresponds to gap draining (the gas film ejection from the region between the droplets); the second and third ones are due to the capillary action caused by the highly curved interface, compare with Figs. 6 and 10. In the right panel of Fig. 11 we show the history of normalised extreme density values and the volume bounded by the interface. The results are presented for t > 2 × 10 −6 s for clarity. From that moment on, the discontinuous initial condition for the velocity which imposes a significant compressible effect (the gas phase density locally decreases by ∼ 3% ) is relaxed. At later times, the minimum density of gas decreases locally (in the vicinity of the interface) by not more than 1.6% . The density of liquid phase is less prone to the cumulation and its maximal deviation from 1 does not exceed 0.1% . The volume of phases bounded by the surface = ⟨ ⟩ is well preserved and varies by ± 0.5% during the simulation.

Steady Water Droplet in Air
Usually, the test of steady droplet is performed in a domain bounded by the solid walls with no-slip boundary conditions but here we use the periodic domain. This choice makes, in fact, the domain infinite so (due to the parasitic currents) there is no reason to perform On the other hand, shorter simulations can also bring some meaningful outcome. We kept the settings from the previous case and considered three spatial resolutions: 128 3 , 192 3 , 256 3 which corresponds to ∼ 43 , 64, and ∼ 85 grid cells per droplet diameter. The simulations were performed until t = 2 × 10 −4 s . The reference velocity scale U was fixed during the simulations and its value was set to be (approximately) the initial droplets speed from the test of head-on collision, i.e. U = 0.6 m∕s.
In Fig. 12 we show the time averaged pressure profiles normalised by the pressure jump from the Young-Laplace law, p YL = 4 ∕D , which in our case is equal to 582.4 Pa. For the time averaging, the data were taken at 100 equally spaced time instants, however some acoustic effects are still visible at the highest resolution. The pressure levels inside the droplet agree well with the theoretical prediction, the maximal deviations are less than 1% and decrease with increasing spatial resolution. On the other hand, sharp local pressure peaks appear in the liquid phase in the vicinity of the interface. This effect is caused by the switch and is more significant at higher resolution when the numerical diffusion is weaker.
In Fig. 13 we present the history of maximal velocity of the spurious currents expressed by the capillary number Ca = 1 | | max ∕ and the total kinetic energy k. The maximal velocity of the spurious currents does not reveal a monotonic behaviour with increasing spatial resolution (at least in the considered range of x ). This is explained as the interplay of the numerical error of the CSS discretisation and the numerical diffusion of momentum. The maximal velocity attains a constant level quite quickly but, unfortunately, it is higher for higher spatial resolution. Interestingly, the case D∕ x ≈ 43 is quite different from the two others-it results in the highest velocity which, additionally, does not "saturate" during the simulation. It seems that in this case the numerical errors stemming from the discretisation of the surface tension are high when related to the amount of numerical diffusion. On the other hand, the higher the spatial resolution, the lower the kinetic energy. This means that, although the velocity magnitudes can grow with increasing resolution, the region affected by the spurious currents becomes smaller. We have also found that the magnitude of the parasitic currents depends significantly on the interface width (the value of ).
The deficiencies reported in this section stem from the sharp switching between the interfacial and bulk regions which was confirmed by setting = 0.1Ma 2 (in the case of steady droplet this also assures the modelling consistency). As an alternative to CSS one could consider the use of the conservative and well-balanced surface tension model proposed recently by Abu-Al-Saud et al. (2018) which is, however, much more complicated than CSS and is not straightforward to be applied in our model. The investigation on the influence of the value and the application of other surface tension models is warranted.

Off-Centre Collision of the Carbohydrate Droplet in Nitrogen
In this section we present a more dynamic case of binary collision of carbohydrate droplets immersed in nitrogen. The main purpose of this simulation is to verify whether the droplet break-up process is simulated correctly. We exactly set the conditions of the experiment of Qian and Law (1997) (case "o" there). The domain size is now 8D × 3D × 3D . The directions of initial velocities of the droplets are parallel and the off-centre parameter is 0.68D. The densities of liquid and gas phases are 1 = 758 kg∕m 3 and 0 = 1.138 kg∕m 3 , respectively. The viscosities are 1 = 2.128 × 10 −3 kg∕(m s) and 0 = 1.787 × 10 −5 kg∕(m s) . The surface tension coefficient is = 0.026 N∕m . The droplets diameters are D = 0.38 mm . The Weber number is equal to 61 and the Reynolds number is 314. The initial velocity of droplets is U 0 ≈ 1.17 m∕s . The simulation is conducted till t = 2.3 ms. We set the spatial resolution to 640 × 240 × 240 , so the cell size is x = 4.75 × 10 −6 m . Initially the droplets are resolved using 80 grid cells per diameter. The grid spacing results in the maximal capillary waves speed c ≈ 4.76 m∕s.
In Fig. 14 we present an overview of the simulation. Till t = 1.84 ms our results agree very well with the experiment, however, due to high numerical diffusion the process is slowed down. The interface evolution is proper since it is not sensitive to the Reynolds number, as pointed by Ashgriz and Poo (1990). At time instants 1.84 ms < t < 2.1 ms , the agreement with the experimental data, where the filament breaks into four droplets, is not good. This is the effect of insufficient spatial resolution. At the final analysed time the interface topology matches well the experimental result when three small droplets are present.
In Fig. 15 the detailed view on the appearance of the first break-up event (at t ≈ 1.4 ms ) is provided. Although the interface curvature is high at this location, the density field remains sharp after the the break-up is completed and no mass diffusion is visible; this artifact would occur if too low value of the coefficient were set.
Finally, in Fig. 16 we show the history of the extremal density values and the volume bounded by the interface. The maximal and minimal density behave similarly to the previous case of the head-on collision. Since the deformation is much stronger, some loss of volume bounded by the interface occurs during the stretching phase after collision. However, the volume is well preserved overall and less than 1% volume is lost.
To summarise the results presented in this section, we point out that: (1) the topological changes of the interface are modelled correctly, no spurious gas entrapment inside the liquid phase nor significant density diffusion are observed; (2) the density variations stemming from the weak compressibility do not exceed 1% which is a commonly accepted criterion to treat the flow as incompressible; (3) the volume bounded by the interface is well preserved; (4) the computed pressure jump across the interface  The off-centre collision of droplets. The history of the normalised extremal density values together with the volume V bounded by the surface = ⟨ ⟩ . V 0 denotes the initial volume of the droplets is in good agreement with the theoretical prediction; (5) modelling of the surface tension with CSS in the presence of the switch causes a locally oscillatory behaviour of the solution; (6) the model parameters' values (in particular ) we used did not allow to obtain the convergence to vanishing velocity field in the case of the steady droplet.

Computational Efficiency
The weakly-compressible flow modelling is inefficient unless the computations can be massively parallelised. This is the only possibility to overcome the severe time step restrictions stemming from the acoustic effects. In the sequential execution one should rather use truly incompressible models.
For this reason, all the 3D computations presented in this work were performed using our in-house EDAC-DI code, dedicated for the execution on GPU. The code was written using NVIDIA-CUDA C API. As a computing device we used a typical desktop computer equipped with a single NVIDIA GTX 1080 TI GPU.
Let us summarise the total memory requirements for the presented algorithm. To store the flow variables, 5 arrays at two time levels are required. To store the normal vectors and gradients of the variables, 18 arrays have to be used. This gives in total 28 arrays. It is still quite economical when compared to the LBM, which for the widely used D3Q19 variant requires 2 × 19 + 4 = 42 arrays for the single phase flow (!). It has to be emphasized that if one uses the traditional, one-dimensional MUSCL reconstruction (using, e.g., the OM scheme which performed also quite well for short-time simulations), the storage requirements of the EDAC-DI solver significantly decrease to 13 arrays, so more than twice. Obviously, one does not need to store the MDIM-limited gradients of primitive variables resulting from the MLP reconstruction, but then they have to be computed twice for each grid cell which significantly increases the execution time since it is the most demanding part of the overall algorithm.
In Table 2 we report the wall-clock time of the two simulations of colliding droplets together with the memory use. In our opinion the computational effort is very attractive keeping in mind that the simulations were performed using a desktop computer.

Summary and Future Work
In this paper we have presented a computational approach for two-phase interfacial flows using a novel mathematical formulation based on the diffuse-interface technique in a weakly compressible framework. Proper values of the model coefficients are found; in particular, the mobility is determined by applying an analogy to the Cahn-Hilliard equation. A robust and efficient numerical method to solve the model equations is proposed. The solution of the DI equation is almost free from the numerical interface distortions. The model has been validated against the steady and colliding liquid droplets immersed in gas, including the topological changes. The outcome of the simulations agrees well with the theoretical and experimental results. The additional advantages of the presented EDAC-DI model are: (1) the use of straightforward and simple numerical solution procedure, (2) very efficient parallelisation of the computational process, and (3) low memory requirements. The need of high spatial resolution, required by the diffuse interface modelling, may be considered as a disadvantage with respect to sharp interface methods such as VoF. Yet, the DI model can deal with the interfacial phenomena without the use of artificial, modelspecific treatments, usually required by the sharp interface methods. In our opinion, the proposed approach to interfacial flows is an attractive alternative not only to other, well established weakly compressible models, such as those used in LBM or SPH, but also to the traditional truly incompressible approach. As stated in the Introduction, the presented paper provides a "proof of concept" of the EDAC-DI model. More work still has to be done and we list the next steps that, in our opinion, make a natural follow-up. From the modelling point of view, the most important issue is the formal analysis and numerical quantitative assessment of the choice of parameter , alike done by Magaletti et al. (2013). The applied CSS model of the surface tension and the chosen implementation of the switch introduce numerical issues and possible alternatives should be investigated. To broaden the range of EDAC-DI applications, an effort on imposing other boundary conditions has to be done, with focus on the wetting of solid surfaces. From the numerical point of view, the following aspects seem to be important: (1) the investigation on more efficient numerical methods, in particular the time integration, that would allow a wider range of the DI model parameters to be used; (2) the use of hybrid schemes to reduce the numerical dissipation apart from the interface region; (3) AMR may be considered to apply our model to more complex flow cases, keeping reasonable computational effort; (4) since AMR is not well suited for the GPUs, the brute-force approach by multi-GPU implementation could possibly appear as a reasonable alternative.