Unified formulation of the momentum-weighted interpolation for collocated variable arrangements

Article history: Received 6 September 2017 Received in revised form 20 March 2018 Accepted 20 August 2018 Available online 24 August 2018

Momentum-weighted interpolation (MWI) is a widely used discretisation method to prevent pressure-velocity decoupling in simulations of incompressible and low Mach number flows on meshes with a collocated variable arrangement. Despite its popularity, a unified and consistent formulation of the MWI is not available at present. In this work, a discretisation procedure is devised following an in-depth analysis of the individual terms of the MWI, derived from physically consistent arguments, based on which a unified formulation of the MWI for flows on structured and unstructured meshes is proposed, including extensions for discontinuous source terms in the momentum equations as well as discontinuous changes of density. As shown by the presented analysis and numerical results, the MWI enforces a low-pass filter on the pressure field, which suppresses oscillatory solutions. Furthermore, the numerical dissipation of kinetic energy introduced by the MWI is shown to converge with third order in space and is independent of the timestep, if the MWI is derived consistently from the momentum equations. In the presence of source terms, the low-pass filter on the pressure field can be shaped by a careful choice of the interpolation coefficients to ensure the filter only acts on the driving pressure gradient that is associated with the fluid motion, which is shown to be vitally important for the accuracy of the numerical solution. To this end, a force-balanced discretisation of the source terms is proposed, that precisely matches the discretisation of the pressure gradients and preserves the force applied to the flow. Using representative test cases of incompressible and low Mach number flows, including flows with discontinuous source terms and two-phase flows with large density ratios, the newly proposed formulation of the MWI is favourably compared against existing formulations and is shown to significantly reduce, or even eliminate, solution errors, with an increased stability for flows with large density ratios.

Introduction
The coupling of pressure and velocity is a key difficulty of simulating incompressible flows and has been a central topic of computational fluid dynamics (CFD) for the past decades [1][2][3]. The difficulties associate with the pressure-velocity coupling can be illustrated by assuming an isothermal, incompressible flow, which is governed by the momentum equations (1) and the continuity equation where ρ is the density, u the velocity, p is the pressure, τ is the shear stress tensor, S are the source terms, t is time and x is the coordinate axis. Aside from the question of how to solve the strongly coupled pressure and velocity fields, the governing equations of a three-dimensional incompressible flow only provide three independent equations for four unknowns (three velocity components plus pressure), which makes the formulation of an equation for pressure based on the governing flow equations non-trivial and has lead to a variety of segregated [1,3,4] and coupled [5][6][7] algorithms. Furthermore, discretising the pressure gradient on the one-dimensional equidistant mesh shown in Fig. 1 using central differencing yields where x is the mesh spacing. The pressure gradient at node P is, crucially, not dependent on the pressure value at node P , irrespective of the algorithm applied to solve the governing equations. Consequently, the governing equations permit two independent pressure fields in a chequerboard pattern [3,4] as a valid solution to the discrete equations, a result that naturally extends to higher dimensions. Pressure-velocity decoupling is a discretisation issue typically associated with incompressible flows. When compressible flows are considered, most numerical frameworks use density as a primary variable, while pressure is determined indirectly via an appropriate equation of state. Although such density-based algorithms are the method of choice when the compressibility of the flow is appreciable, they are ill-suited for flows with low Mach numbers [2,8], in particular in the incompressible limit. Motivated by the desire to compute flows at all speeds with the same numerical framework, a number of pressure-based algorithms for flows at all speeds have been developed, e.g. [9][10][11][12]. However, the insignificant compressibility of flows with low Mach number admits pressure-velocity decoupling in the compressible flow solution on meshes with collocated variable arrangement.
Historically, pressure and velocity were coupled by staggering the points at which pressure and velocity are evaluated, the staggered variable arrangement, as proposed by Harlow and Welch [13], with velocity typically evaluated at the centres of the cell-faces, while all other variables are evaluated and stored at the cell centres. A staggered variable arrangement enforces a natural coupling between pressure and velocity, and yields a very compact stencil for the pressure gradient that drives the velocity at the adjacent cell centres through the momentum equations. There is no doubt that for Cartesian meshes, a staggered variable arrangement is efficient and effective. However, as CFD has matured as a tool, it has found ever more frequent application to analyse flows in complex geometries, represented by unstructured meshes, for which the application of a staggered variable arrangement is difficult and may include complex corrections to account for meshes of relatively poor quality [14,15]. This difficulty, in conjunction with the bookkeeping overhead associated with staggered variable arrangements [16], has motivated the development of discretisation methods for collocated variable arrangements, in which all variables are stored at cell centres, that prevent the pressure-velocity decoupling ensuing as a result of the scenario presented in Eq. (3). Notable methods that allow robust computations on meshes with collocated variable arrangement are the momentum-weighted interpolation (MWI), based on the work of several researchers in the early 1980s [17] and widely attributed to Rhie and Chow [18], the artificial compressibility method [19] and one-sided differencing [20], of which MWI is by far the most widely used at present [21].
The principle of the MWI, also frequently referred to as pressure-weighted interpolation or Rhie-Chow interpolation, is to evaluate the velocity at the faces based on weighting coefficients that are derived from the discretised momentum equations, including pressure gradients. By construction, the MWI emulates a staggered variable arrangement, introducing a cell-to-cell pressure coupling and implementing a low-pass filter acting on the third and higher derivatives of pressure [2,3,22,23] to suppress pressure-velocity decoupling, while preserving the second-order accuracy of traditional finite-volume methods [2,3,24,25]. The MWI as originally proposed by Rhie and Chow [18], however, only considers the coupling between pressure gradients with the advection and shear stress terms of the momentum equations, neglecting contributions of the transient term, source terms and originating from underrelaxation. A range of modifications to the original formulation of the MWI have been proposed that account for underrelaxation [26,27] and transient terms [16,23,[28][29][30][31][32][33]. Recently, Xiao et al. [12] showed that neglecting transient terms in the MWI for the simulation of unsteady problems results in a dispersion error for the propagation of pressure waves in compressible flows.
Including source terms in the MWI was discussed by Rahman et al. [34] with the motivation to maintain a strong pressure-velocity coupling, which may otherwise be masked when the source terms are large, due to the direct coupling between source terms and the pressure gradient. Subsequently, van Wachem and Gopala [35] and Mencinger and Zun [36] demonstrated that the inclusion of source terms follows directly from the governing equations, by presenting coherent derivations of the MWI from the momentum equations, and demonstrated the efficacy of the proposed formulation using multiphase flows, especially two-phase flows with surface tension that yield sharp discontinuous source terms. Building upon this work, Denner and van Wachem [37] presented an MWI formulation for flows with source terms, and including transient contributions, on arbitrary meshes. When computing flows with source terms, the MWI necessitates a force-balanced discretisation [24] to avoid the production of spurious velocities as a result of a mismatch of the discretisations applied to the pressure gradients and the source terms. In a force-balanced discretisation the pressure gradients and the source terms are discretised with equivalent methods, so that the forces applied to the flow by the source terms can be precisely balanced by the corresponding pressure gradients. More recent work [25,38,39] has focused on source terms arising in porous media, which follows a similar procedure as the inclusion of source terms in multiphase flows. Curiously, much of the published work on MWI remains focused on Cartesian meshes, although modifications required for arbitrary meshes have already been proposed [24,28,36,37,[40][41][42]. This spread of separate developments, typically focusing on a single aspect of MWI, has led to a large number of subtly differing approaches; to this date, a unified and consistent formulation of the MWI is not available.
The objective of this work is the derivation of a unified and consistent formulation of the MWI, which is applicable to arbitrary meshes (structured and unstructured meshes), with discontinuous source terms and varying fluid densities, suitable for the simulation of single-phase and multiphase flows. The MWI is derived from the discretised momentum equations, which leads to a consistent formulation of the MWI and provides a firm theoretical basis. The presented analysis and numerical results show that the key property of the MWI is a low-pass filter enforced on the third and higher derivatives of pressure, including a cell-to-cell pressure coupling, which suppresses oscillatory solutions, while maintaining second-order accuracy with respect to the mesh spacing. In order for this filter to be retained in flows with source terms, the discretisation has to ensure that the low-pass filter is only applied to the driving pressure gradient that is associated with the fluid motion, by carefully accounting for the source terms in the MWI. To this end, a force-balanced discretisation of the source terms is proposed, that precisely matches the discretisation of the pressure gradients for smooth as well as discontinuous source terms, and preserves the force applied to the flow. A range of representative test cases are used to scrutinise the efficacy of the proposed formulation of the MWI, and to compare it to previously published formulations of the MWI. The proposed formulation of the MWI is shown to provide a robust pressure-velocity coupling, even for flows on meshes of poor quality, and for flows with discontinuous source terms, as well as discontinuous density changes of up to six orders of magnitude, with similar or reduced errors compared to existing MWI formulations.
The applied numerical frameworks for incompressible and compressible flows are briefly introduced in Section 2. In Section 3, the MWI is derived from the discretised momentum equations for arbitrary meshes and validated on structured and unstructured meshes. An extension of the MWI for the inclusion of source terms is proposed in Section 4 and the density weighting of the MWI for flows with discontinuous changes in density is discussed in Section 5. Based on the presented step-by-step analysis of the MWI, a unified formulation of the MWI is proposed in Section 6. The article is concluded in Section 7.

Numerical framework
The governing equations for incompressible flow, Eqs. (1) and (2), are discretised using the finite-volume method, with the semi-discretised equations for cell P , shown schematically in Fig. 2, given as where subscript f denotes the faces of cell P , A f is the area of face f , n f is the unit normal vector of face f (outward pointing with respect to cell P ), ϑ f = u fn f is the advecting velocity, S P are the discretised source terms, and V P is the volume of cell P . In this study, the transient term in the momentum equations, Eq. (4), is discretised using the first-order or second-order backward Euler schemes, while the face velocity u j, f is evaluated using the central differencing scheme with an implicit correction for mesh skewness [43]. Deriving a consistent discretisation for the advecting velocity based on the MWI, where u i, f is interpolated from the adjacent cell centres, and constructing a force-balanced discretisation of the source term S P are the main objectives of this study. The MWI presented in Section 3 and its extensions to flows with source terms and density discontinuities in Sections 4 and 5, respectively, are tested using the fully-coupled finite-volume framework for single-phase and multiphase flows on arbitrary meshes of Denner and van Wachem [37]. The low Mach number flows presented in Section 3.6 are simulated using the pressure-based finite-volume framework for single-phase flows at all speeds of Xiao et al. [12]. For compressible flows the momentum equations are and the continuity equation is In addition to the momentum and continuity equations, compressible flow is governed by the energy equation where h = c p T + u 2 /2 is the total enthalpy, c p is the specific heat capacity at constant pressure, and T is the temperature.
Without loss of generality, heat conduction is neglected in this work. The considered fluid is assumed to be an ideal gas with the density given by the equation of state with γ = c p /c v the heat capacity ratio and c v the specific heat capacity at constant volume. In particular, the continuity equation, Eq. (8), is discretised as [12] ∂ρ ∂t following a Newton-linearisation of the advection term, where superscript (i) denotes values that are implicitly solved for and superscript (i − 1) denotes deferred values from the previous iteration. This Newton-linearisation of the advection term of the discretised continuity equation, Eq. (11), allows simulations without MWI, contrary to the employed incompressible framework. The interested reader is referred to the work of Xiao et al. [12] for a detailed description of this numerical framework and discretisation.

Momentum-weighted interpolation
The MWI is derived from the momentum equations with the aim of providing a consistent formulation of the advecting velocity ϑ f , Eq. (6), for arbitrary meshes and to analyse the general properties of the MWI. Using an appropriate approximation for the values at face centre f , such as a linear interpolation of values at adjacent cell centres, a first-order Euler scheme to discretise the transient term and neglecting source terms, Eq. (4) is given at cell centre P as where t is the applied time-step, superscript O denotes values from the previous time-level, subscript F represents the neighbour cell of cell P that is adjacent to face f , as schematically illustrated in Fig. 2, and a is the sum of the coefficients of the advection term and the shear stress term arising from the discretisation applied to the momentum equations. By defining Eq. (12) can be rewritten as Note that in the analysis presented here, a backwards Euler scheme with first-order accuracy is used for time integration; the extension to higher-order accurate transient schemes is straightforward [32,44]. Assuming Eq. (17) can be similarly formulated for any control volume, an equivalent equation is written for cell F and, in analogy to a staggered control volume, at face f With the aim of obtaining an expression for the velocity u j, f , ũ j, f is approximated by linear interpolation as where l f is the linear interpolation coefficient. Note that this definition for ũ j, f is chosen so that the expression for u j, f is time-step independent, i.e. the steady-state solution should not contain any terms that are dependent on the time-step and the MWI should recover the steady-state solution if a transient flow reaches a steady state [45]. Substituting Eqs. (17) and (18) into Eq. (20), the face velocity u j, f , given by Eq. (19), becomes where 2 f denotes a value at face f obtained by linear interpolation from the adjacent cell centres. With t → ∞ for a steady-state solution, c → 0 and Eq. (21) reduces to which is independent of the time-step t.

Velocity interpolation
The face velocity u i, f is interpolated with a linear interpolation is the interpolation coefficient obtained by inverse distance weighting, indicated by the superscript "(idw)", and q P f and q F f are the distances between cell centre P and face f and between cell centre F of the adjacent cell and face f , respectively. The accuracy of the interpolated face velocity u i, f reduces when the interpolation point f along the vector connecting the cell centres P and F does not coincide with the centre of face f , see Fig. 2, commonly referred to as mesh skewness. In order to retain the accuracy of the linear interpolation, a gradient-based correction is added to the interpolation of the face velocity [22,46], with the face velocity following as where r f is the vector from interpolation point f to face centre f , see Fig. 2. The velocity gradient ∂u i /∂x j | f is typically obtained by linear interpolation (e.g. using inverse distance weighting) from the velocity gradients at the adjacent cell centres. Note that the precise type of interpolation is not critical and the correction of mesh skewness is optional, as it does not have a direct influence on the efficacy of the MWI in providing a robust pressure-velocity coupling, although the interpolation given in Eq. (25) is desirable with respect to the accuracy of the interpolation on non-Cartesian meshes [21,46].

The pressure gradients
The discretisation and interpolation of the pressure gradients in the MWI has a direct influence on the low-pass filter properties with regards to the pressure field, which is widely considered to be the key characteristic of the MWI [3,24,30,34,37], and the efficacy of the associated pressure-velocity coupling. First, the interpolation on an equidistant mesh is discussed, and the low-pass filter on the pressure field is derived, to illustrate the key concepts of the MWI. This is followed by a generalisation of the interpolation and the low-pass filter to non-equidistant meshes.

Equidistant mesh
The applied numerical framework is based on a finite-volume method, so that the straightforward discretisation of the pressure gradients on the one-dimensional equidistant mesh shown in Fig. 1 follows as As the mesh is equidistant, a linear interpolation of the pressure gradients at cells P and E with a 1/2-weighting is given as Inserting the discretised pressure gradients given in Eqs. (26)- (28), including the interpolation given in Eq. (29) and neglecting the weighting term d f for simplicity, the pressure term of Eq. (22) at face f is given as For comparison, the third-order derivative of pressure at face f is given as which shows that the pressure term in Eq. (30) is proportional to the third-order derivative of pressure, It is this relationship that dampens out non-physical pressure oscillations on meshes with collocated variable arrangement [3,4,24,30,34]. Approximating the pressure gradient at face f with standard finite differences as in Eq. (26), provides a spatial cell-to-cell coupling of the pressure field and matches the discretised pressure gradient that would be employed if f would correspond to a control volume in a staggered variable arrangement. Moreover, the pressure term is proportional to x 2 , see Eq. (32), and, hence, the second-order accuracy of the finite-volume framework is preserved [2]. The extension to multiple dimensions is straightforward by computing the cell-centred gradients using the divergence theorem, given for cell P as where f are all faces adjacent to cell P , and analogously for cell E.

Non-equidistant meshes
The choice of interpolation coefficient for the pressure gradients on non-equidistant meshes is a controversial issue in the literature and has not been conclusively settled. The use of a linear interpolation with weighting coefficients based on the mesh geometry is frequently advocated [10,11,16,28,32,36,41,42,47], e.g. inverse distance weighting [28,32,36,41] or volume weighting [47], for the interpolation of the pressure gradients. However, previous studies [3,24,30,37] suggested the use of the 1/2-weighting given in Eq. (29) also for non-equidistant meshes, in order to retain the filter properties of the MWI. Pascau [16] also suggested the use of a harmonic average for all interpolations in the MWI, but did not further elaborate on the suitability of harmonic averaging.
Consider the example illustrated in Fig. 3, where the mesh spacing suddenly changes by a factor of x E / x P = 5 in the cells adjacent to face f . Applying inverse distance weighting to interpolate the cell-centred pressure values to faces e, f and g, the derivatives required for the pressure term of the MWI follow as and the pressure term of the MWI, which is also interpolated with inverse distance weighting, becomes The third derivative of pressure in this case is so that the pressure term of the MWI is Similar relationships can be found in the same manner for any ratio x E / x P and it can, therefore, be concluded, that the pressure term of the MWI is formally equivalent to the corresponding third derivative of pressure, if inverse distance weighting is applied consistently for all interpolations of pressure and its gradients. If the cell-centred pressure gradients are, however, evaluated with face values obtained with inverse distance weighting, see Eqs. (35) and (36), but interpolated in the MWI with 1/2-weighting, the resulting pressure term, for instance given for is not proportional to the corresponding third derivative of pressure, Eq. (38), and does not formally provide the low-pass filter on the pressure field. Yet, if the 1/2-weighting is applied throughout, i.e. for the interpolation of p e , p f and p g as well as in the interpolation of the cell-centred pressure gradients in the MWI, the filter on the pressure field, then given by Eq. (32), would be retained, albeit at the cost of reducing the accuracy of the pressure gradient evaluation. This reduced accuracy of the pressure gradient introduces an error in the momentum equation that, based on the Taylor expansion where superscript "(1/2)" denotes interpolation with 1/2-weighting, is proportional to the distance x and, thus, increases linearly with increasing cell-size ratio of adjacent cells. The influence of these different interpolation techniques is further analysed from a practical viewpoint using numerical results in Sections 3.6.1 and 3.6.4.

Weighting coefficients c and d
In Eq. (21) the weighting coefficients c and d appear in various forms, at face f as well as interpolated to face f from the values at adjacent cells P and F . This can be further simplified by observing that the pressure term of the MWI has to vanish and u f = u f , if the pressure gradient is constant or varies linearly (assuming a steady state),

an initially transient flow assumes a steady-state solution.
To ensure u f = u f if the gradient of pressure is constant or varies linearly, the coefficient of the interpolated face velocity u f has to be unity. Hence, (43) and similarly for the coefficients of the previous time-step the pressure term of the MWI is given as Therefore, for the pressure term of the MWI to vanish for a constant or linearly varying pressure gradient, d f becomes which is satisfied by the approximation While the approximation given in Eq. (50) has been applied in a somewhat ad-hoc manner by other researchers [11,23,25,30,32,37,42], generally without justification beyond the difference d P − d F is assumed to be small [30,34], the above analysis shows that it is necessary to obtain a physical solution.

MWI formulations
Inserting the approximations defined in the previous section, cf. Eqs. (43), (44), (48) and (50) whered The corrections to the interpolated velocity provided by the pressure term and the transient term in Eq. (51) vanish if the pressure gradient is constant or varies linearly, and if the flow assumes a steady state. Similarly, applying the approximations for d f given in Eqs. (48) and (50) into Eq. (22), the face velocity derived from the steady momentum equations follows as which corresponds to the MWI formulation proposed by Rhie and Chow [18].
It is important to note that the MWI formulation of Rhie and Chow [18] is derived from the steady momentum equations, i.e. ∂u/∂t = 0. This has sometimes lead to misunderstandings in the literature [21,28,30], where the steady MWI formulation was applied based on the coefficients of the transient momentum equations. Instead of applying coefficient d as defined in Eq. (14), i.e. based on the coefficients of the advection term and the shear stress term arising from the discretisation of the momentum equations, these formulations apply thus, including the transient coefficients as well, by dividing Eq. (12) through the coefficient of u j,P on the left-hand side of Eq. (12). Applying the approximations given in Eqs. (48) and (50) to obtain d * f , the face velocity is then defined without the transient term, i.e. using the formulation derived from the steady momentum equations by Rhie and Chow [18], as The coefficient d * f is then responsible for a time-step dependency of the MWI, with pressure-velocity decoupling reported for small time-steps t [21]; the pressure term vanishes for small t because lim t→0 d * f = 0. Choi [28] remedied this timestep dependency by applying d * f using an MWI formulation consistently derived from the transient momentum equations, starting with Eq. (12), to yield which was shown to be time-step independent [28]. The MWI formulation in Eq. (56) is almost identical to the formulation given in Eq. (51); both are derived from the transient momentum equations, and the difference between coefficients d f , Eq. (52), and d * f , Eq. (54), is minute and for practical applications irrelevant. However, unlike the transient formulation of the MWI given in Eq. (51), the formulation of Choi [28] given in Eq. (56) does not take into account temporal changes in density ρ for the coefficient of the transient term, which may lead to errors in a flow with strong changes of density (e.g. transonic flows) [12].

Advecting velocity
Noting that the face velocity u f defined using the MWI appears in the discretised governing equations, Eqs. (4) and (5) for incompressible flows and Eqs. (7)-(9) for compressible flows, as the dot product with the face normal vector n f , an advecting velocity ϑ f = u fn f can be defined (cf. [3,37,42]). This advecting velocity is given for the MWI formulation presented in Eq. (51) as and similarly for the other MWI formulations discussed in Section 3.4.
When the vectors ŝ f and n f are not parallel, as for instance in the example given in Fig. 2, the pressure gradients at cell centres are misaligned to the pressure gradient at the face, because where s f is the distance between cell centres P and F ( s f = x on an orthogonal mesh where ŝ f =n f ). Consequently, the pressure term of the MWI is no longer guaranteed to constitute a low-pass filter with respect to the third and higher derivatives of pressure. The pressure filter can be restored by applying a deferred-correction approach [22,48], as previously applied to the MWI by Zwart [40], Ham and Iaccarino [49] and Denner and van Wachem [37], and similarly proposed by Ferziger and Perić [3], decomposing the product ∇p fn f into an orthogonal and a non-orthogonal part. The pressure gradient at face f is then defined as where α f is the scaling factor of the decomposition. Inserting Eq. (59) into the pressure term of Eq. (57), the pressure term for non-orthogonal meshes follows as Hence, the correction ensures that the entire pressure term of the MWI is projected onto the vector ŝ f connecting the adjacent cell centres, a prerequisite to retain the filter properties of the MWI on arbitrary meshes, as shown in Section 3.6. Three basic decompositions for the non-orthogonal correction are readily available, determined by the choice of the scal- ing factor α f , as illustrated in Fig. 4. The straightforward choice is to apply α f = 1, with which the filter properties of the MWI are independent of the angle between n f and ŝ f . Zwart [40] and Ham and Iaccarino [49] suggested to use α f =n fŝ f , which reduces the weight of the pressure filter with increasing non-orthogonality of the mesh. In contrast, Denner and van Wachem [37] applied α f = (n fŝ f ) −1 , with which the weight of the pressure filter increases with mesh non-orthogonality. Previous studies [37,40] chose the definition of α f with reference to the literature on the deferred correction of diffusion terms [22,50,51]. Although the underpinning motivation is the same, the deferred-correction approach for diffusion terms aims at mitigating a directional bias imposed on diffusive transport, while the motivation of the deferred-correction approach for the MWI is to retain the filter properties of the pressure term. To this end, the scaling factor α f becomes merely a weighting factor on the pressure term of the MWI, as seen in Eq. (60). The influence of the choice of α f on the efficacy of the MWI is further discussed in Section 3.6.3. In summary, applying the approximations presented above, the advecting velocity for steady-state and transient flows on arbitrary meshes is given as where inverse distance weighting should be applied for the interpolations. Applying this advecting velocity in the discretisation of the governing equations, Eqs. (4) and (5), assuming source terms (separately discussed in Section 4) are negligible, enforces the spatial coupling of the pressure field that is otherwise missing. A careful choice of the approximations results in a low-pass filter that targets third-order and higher oscillations of pressure, that appear in a decoupled pressure field.

Numerical experiments
Four representative test cases are considered to assess the characteristics and properties of the MWI: the propagation of acoustic waves in Section 3.6.1, the propagation of a non-monochromatic pressure pulse in Section 3.6.2, the flow in a lid-driven cavity in Section 3.6.3, and Taylor vortices in an inviscid fluid in Section 3.6.4. The propagation of acoustic waves allows a detailed analysis of the effect of the pressure interpolation and the filter properties of the MWI, while the non-monochromatic pressure pulse highlights the importance of the transient term of the MWI. The lid-driven cavity demonstrates the versatility and robustness of the proposed formulation and the Taylor vortices enable an in-depth analysis of the influence of the MWI on the kinetic energy conservation. Unless stated otherwise, the advecting velocity ϑ f given in Eq. (61), based on the transient MWI formulation, Eq. (51), is applied.

Propagation of acoustic waves
The propagation of acoustic waves is simulated on one-and two-dimensional non-equidistant meshes. Three different formulations of the pressure gradients in the MWI and in the momentum equations, based on the analysis presented in Section 3.2, are considered: 1. The pressure at faces and the discretised pressure gradients in the MWI are interpolated with inverse distance weighting, abbreviated "p-idw, MWI-idw"; 2. The pressure at faces is interpolated with inverse distance weighting, while the discretised pressure gradients in the MWI are interpolated with 1/2-weighting, abbreviated "p-idw, MWI-1/2"; 3. The pressure at faces and the discretised pressure gradients in the MWI are interpolated with 1/2-weighting, abbreviated "p-1/2, MWI-1/2".
It is important to remember at this point, as previously discussed in Section 3.2, that the same discretised cell-centred pressure gradients are applied in the momentum equations and the MWI, and that the interpolated pressure at faces is used to determine these discretised cell-centred pressure gradients via the divergence theorem, see Eq. assumed to be inviscid, meaning the amplitude of the acoustic waves is not attenuated by viscous stresses. Since the density and velocity amplitudes are small, with ρ 0 ρ 0 and u 0 u s,0 , the theoretical pressure amplitude of the acoustic waves follows from linear acoustic theory [52] as p 0 = u s,0 ρ 0 u 0 = 4.32 Pa.
The one-dimensional domain is represented by a mesh with a sharp change in mesh spacing at x = 0, similar to the mesh shown in Fig. 3, changing from a small mesh spacing x S to a large mesh spacing x L = λ 0 /20. Fig. 5 shows the pressure profiles on the meshes with x L / x S ∈ {5, 20} for the three considered formulations of the pressure gradients.
In all cases, the pressure profile is not visually affected by the choice of the pressure gradient formulation or the cell-size ratio of the mesh, with the predicted amplitude of the pressure wave in the range 4.30 Pa ≤ p ≤ 4.35 Pa, which is in excellent agreement with the theoretical value of p 0 = 4.32 Pa. Interestingly, the pressure amplitude of the acoustic waves does not diminish as they propagate through the domain, indicating that the proposed formulation of the MWI does not introduce spurious pressure contributions or damping that alters the pressure field, provided that the pressure waves are appropriately resolved in space and time. However, the profiles of the velocity gradient, shown in Fig. 6, reveal distinct differences between the pressure gradient formulation that uses exclusively inverse distance weighting ('p-idw, MWI-idw'), which does not exhibit any dependency on the cell-size ratio x L / x S , and the other two considered formulations of the pressure gradient ('p-idw, MWI-1/2' and 'p-1/2, MWI-1/2'), which exhibit a considerable error at the position of the change in mesh spacing (x = 0). In fact, Fig. 7 shows that this error grows linearly with x L / x S , as expected from Eq. (42).
Simulating the propagation of these acoustic waves on the hybrid quadrilateral/triangular two-dimensional mesh shown in Fig. 8  Although, contrary to the corresponding one-dimensional case, this error is not entire eliminated with the 'p-idw, MWI-idw' formulation, which is attributed to mesh skewness, the error is significantly reduced.  In the formulation using 'p-idw, MWI-1/2', the error originates in the MWI, because the pressure term of the MWI is not equivalent to the corresponding third derivative of pressure and, hence, violates the filter of the MWI, cf. Eqs. (38) and (41). For the formulation with 'p-1/2, MWI-1/2', the pressure term of the MWI is formally equivalent to the corresponding third derivative of pressure and satisfies the low-pass filter on the pressure field. However, due to the interpolation of the face values of pressure with 1/2-weighting, the cell-centred pressure gradients adjacent to the change in mesh spacing are inaccurate and, therefore, introduce an error in the momentum equations. Interestingly, the error introduced in the MWI by the 'p-idw, MWI-1/2' formulation and the error introduced in the momentum equations by the 'p-1/2, MWI-1/2' formulation have identical magnitudes, see Fig. 7. This implies that the formulation of the MWI presented in this section is consistently derived from the momentum equations.   Despite the seemingly significant theoretical differences between the considered pressure term formulations, see Section 3.2, the impact on the accuracy of the results is very modest. Although the propagation of acoustic waves is a very sensitive test case with respect to the applied numerical methods [12,45], as even small inconsistencies or a lack of convergence lead to a visible change in the amplitude and speed of the waves, the propagation of the acoustic waves is predicted with high accuracy with all considered formulations. This suggests that the interpolation coefficients of the linear interpolation of the cell-centred pressure gradients in the MWI is not a primary factor for a robust pressure-velocity coupling on meshes with reasonably smooth changes of mesh resolution.

Propagation of a pressure pulse
The      Fig. 10a, the pressure pulse is simulated accurate at all considered Co. Using the MWI formulation without the transient term, however, leads to an appreciable error in the prediction of the pressure profile, see Fig. 10b, which was similarly observed by Xiao et al. [12]. In fact, without the transient term in the MWI, the amplitude of the pressure error increases steadily with Co, as seen in Fig. 11. Although a small error at the back of the pulse can also be observed when the transient term is included in the MWI, in particular for Co = 0.5, the magnitude of the error is substantially smaller. Consequently, neglecting the transient term in the MWI ensues a dispersion error in pressure signals.
The transient term of the MWI is, thus, essential for an accurate, and largely time-step independent, prediction (assuming an appropriate time resolution, typically Co < 1) of the propagation of pressure waves and disturbances.

Lid-driven cavity
A lid-driven cavity, schematically shown in Fig. 12, is simulated to demonstrate the impact of the different MWI formulations. The two-dimensional domain has the dimension 1 m × 1 m, with the top wall moving at a constant velocity of u w = 1 m s −1 . As indicated in Fig. 12, a no-slip condition is applied at the top wall, while all other walls are assumed to  On all three meshes a stable and oscillation-free result is obtained. Furthermore, the profiles of both velocity components u  and v, shown in Fig. 15, are in very good agreement on all meshes, demonstrating the efficacy of the proposed formulation of the MWI on arbitrary meshes. Note that the scaling factor of the non-orthogonal correction for the results obtained on the triangular mesh and the non-orthogonal hexahedral mesh, see Figs. 14b and 14c, is α f = 1; simulations with α f =n fŝ f and α f = (n fŝ f ) −1 yield results with no appreciable difference for this case and are, thus, omitted for clarity. The pressure along the y-centreline of the domain is shown in Fig. 16a on Cartesian meshes with different resolution, indicating that the pressure profile converges with increasing mesh resolution. In fact, as observed in Fig. 16b, the discretisation of the advecting velocity with MWI does not affect the second-order accuracy of the applied finite-volume method, confirming the theoretical analysis in Section 3.2. When the simulation on the non-orthogonal mesh is restarted applying the advecting velocity ϑ f as given by Eq. (57), i.e. without the projection of the pressure term in the MWI presented in Eq. (60) to correct for mesh non-orthogonality, yields significant pressure oscillations after only one time-step t = 5 × 10 −3 s, as observed in Fig. 14d, and the solution algorithms diverges a few time-steps later. This shows clearly the substantial improvement of accuracy and stability provided by the non-orthogonal correction presented in Eq. (60), and demonstrates the consequences when the low-pass filter of the MWI is severely compromised.
The utility of the MWI to eliminate pressure-velocity decoupling becomes strikingly apparent when the lid-driven cavity is considered with a compressible fluid; with and without MWI. The previous simulations are modified such that the fluid has a heat capacity ratio of γ = 1.4 and a specific heat capacity at constant volume of c v = 720 J kg −1 K −1 , which approximately corresponds to the properties of air at room temperature. The flow has an initial pressure of p 0 = 10 5 Pa and an initial temperature of T 0 = 347.22 K, so that the density is ρ 0 = 1 kg m −3 . Hence, the speed of sound is u s,0 = √ γ p 0 /ρ 0 = 374.17 m s −1 , which corresponds to a Mach number of M = u w /u s,0 = 2.67 × 10 −3 . Fig. 17a shows the pressure contours at steady state using the advecting velocity ϑ f given in Eq. (61), on the equidistant Cartesian mesh with 50 × 50 cells. As expected, unphysical pressure oscillations as a result of pressure-velocity decoupling are absent and the pressure distribution is in excellent agreement with the incompressible result shown in Fig. 14a. Restarting the compressible simulation with this result but omitting the MWI in the formulation of the advecting velocity, so that ϑ f = u i, fni, f , clearly discernible pressure oscillations develop, as seen in Fig. 17b. Note that the result shown in Fig. 17b is an instantaneous so that Eqs. (57) and (62) are identical for β = 1, allows to highlight the influence of the pressure-velocity coupling provided by the MWI. Fig. 19 shows the pressure contours at steady-state for β = 0.1 and β = 0.01 with Mach number M = 2.67 × 10 −3 . The onset of pressure-velocity decoupling is clearly visible with β = 0.1 in the upper right quadrant of Fig. 19a, where pressure oscillations can be observed that are absent with the unmodified MWI (β = 1), shown in Fig. 17a. The pressure-velocity decoupling amplifies for even smaller values of β, as observed for β = 0.01 in Fig. 19b. Contrary to the result shown in Fig. 17b, the results shown in Fig. 19 are at steady state, i.e. the pressure oscillations are stable and do not further grow. With respect to the pressure term of the MWI, the coefficient β in Eq. (62) has a similar effect as the scaling coefficient α f of the non-orthogonal correction. Considering that β has to be reduced significantly for pressure-velocity decoupling to set in, it is not surprising that the different definitions of the scaling coefficient α f discussed in Section 3.5 have little influence on the results, since even for a severely deteriorated mesh with an angle between n f and ŝ f of 60 • , the product of n fŝ f is only 0.5.

Taylor vortices
Two-dimensional Taylor vortices in an inviscid fluid are simulated to analyse the dissipation of kinetic energy by the MWI. The conservation of kinetic energy is a fundamental property arising from the conservation of mass and momentum, i.e. the governing flow equations, and is associated with a robust solution [53]. Following the work of Ham and Iaccarino [49], the domain has the dimensions 2 m × 2 m and the initial conditions are u = − cos(π x) sin(π y), v = sin(π x) cos(π y), p = − where is the volume of the computational domain, should be constant. However, using this test-case, Ham and Iaccarino [49] identified an unphysical numerical dissipation of kinetic energy associated with the pressure term of the MWI. Comparing the error in kinetic energy  [49]; the MWI clearly dissipates kinetic energy. After an adjustment immediately following the start of the simulation, the kinetic energy remains constant using the advecting velocity formulated without MWI, indicating that the numerical dissipation of kinetic energy is negligible. The oscillations   Fig. 21b arise from a periodical compression and decompression of the flow due to its significant compressibility, which is undamped because of the lack of physical dissipation (the fluid is inviscid and heat conduction is neglected). Applying the advecting velocity with MWI as given in Eq. (57), ε kin is several magnitudes larger and monotonically increasing, but similar for both considered Mach numbers. Note that, for the shown case, pressure-velocity coupling sets in for M 0.01.  57), as well as based on the steady MWI formulation given in Eq. (53), which corresponds to the formulation proposed by Rhie and Chow [18]. The error in kinetic energy ε kin converges with third order under mesh refinement for both MWI formulations. This third-order convergence may be at first sight surprising, given the pressure term of the MWI is proportional to x 2 , as shown in Section 3.2, and considering that Ham and Iaccarino [49] reported second-order convergence of the error in kinetic energy for the same test-case. However, Ham and Iaccarino [49] did not consider the coefficient d as defined in Eq. (14), which is d ∝ x. Hence, the product of d f and the pressure term is and similarly if    the errors produced by the transient and steady MWI formulations. For a flow field that evolves slowly in time, such as the quasi-steady Taylor vortices (the flow is steady apart from the dissipation introduced by the MWI) considered here, this is to be expected. In Section 3.6.1, the applied coefficients of the interpolation of the cell-centred pressure gradients in the MWI are shown to have a small but appreciable influence on the flow field, with the inverse distance weighting providing superior results compared to 1/2-weighting. A similar effect on the flow field can be observed for the Taylor vortices, shown in Fig. 24 on a mesh with a sharp change of mesh resolution ( x L / x S = 5), with inverse distance weighting and 1/2-weighting of the cellcentred pressure gradients in the MWI; the velocity gradient has a noticeable discontinuity with the 1/2-weighting where the mesh spacing changes. Fig. 25 shows the spatial error in kinetic energy ε kin at t = 1 s on the same non-equidistant mesh with different mesh resolutions, using the two considered interpolation methods. The error in kinetic energy ε kin is almost identical with inverse distance weighting and 1/2-weighting, and both converge with third order, similar the results obtained on the equidistant Cartesian meshes in Fig. 22a.

Source terms
The advecting velocity defined by Eq. (61) is applicable to steady-state and transient flows on arbitrary meshes and, as long as sources terms vary smoothly, this is sufficient for many applications. However, in the presence of source terms that are discontinuous or, more generally, have large gradients, previous studies [35][36][37] have shown that the effect of these source terms on the pressure field have to be accounted for in the MWI.
The reason for including these source terms in the MWI can be illustrated by assuming a quiescent flow, with an external force applied by means of a source term S . The momentum equations, Eq. (1), in semi-discretised form reduce for this case to The corresponding advecting velocity, Eq. (61), which is ϑ f = 0 since the flow is quiescent, becomes a relationship previously used by Rahman et al. [34] for the inclusion of source terms in the MWI. However, Eq. (68) can only be satisfied if the source term results in a uniform or linearly varying pressure field. Hence, in order for the discretisation to be truly force-balanced, meaning that the discretisation of the pressure gradients and the source terms are equivalent, a source term S has to be constructed that can match the discretised pressure gradient in all circumstances. Equation (66) suggests that source terms cause a pressure gradient, an observation also exploited in previous studies [25,[35][36][37], in addition to the pressure gradient associated with the underlying flow ∇ p, henceforth called the driving pressure gradient, with the pressure gradient being It is the pressure gradient associated with the velocity field, i.e. the driving pressure gradient, that is relevant for the pressure-velocity coupling, while all other contributions to the pressure gradient, i.e. source terms, should be excluded. The advecting velocity including source terms should, hence, be Therefore, the discretisation of the source terms has to precisely match the discretisation of the pressure gradients to avoid spurious corrections of the MWI that manifest as unphysical fluid accelerations.

Discretisation of source terms
The analysis of the pressure term in Section 3.2 shows that, to preserve the low-pass filter on the pressure field, the cell-centred pressure gradient has to be evaluated using the divergence theorem, given by Eq. (33). This can be reformulated in terms of p f as with w f = (1 − l f )n f A f . Using a finite-volume method, the last term on the right-hand side of Eq. (71) is by definition zero and, thus, Eq. (71) is equivalent to Eq. (33). If the interpolation of the pressure at face f is amended by additional correction terms, such as the gradient-based skewness correction in Eq. (25), the general formulation for the pressure gradient reads where k f is the correction coefficient of the interpolation, e.g. k i, f = r f ∇p fni, f A f for mesh skewness. The discretisation of the source terms should follow the same template as the discretisation of the pressure gradients, with the cell-centred source terms discretised as similar to the definition of the pressure gradient at face centres given in Eq. (58). The discretisation of the source term given in Eq. (73) is consistent with the discretisation of the pressure gradients and, therefore, ensures a force-balanced discretisation. However, this discretisation may modify the actual force applied to the fluid by the discretised source term. The interpolation of the source term at face f from the cell-centred values at P and F has to ensure that the discretised source terms apply the correct force to the flow. Returning to the simplified onedimensional, quiescent flow and assuming, for now, that Eq. (66) is satisfied by the discretised source terms and pressure gradients, the pressure difference p f = p F − p P is given by with q P f and q F f the distance between cell centre P and face f and between cell centre F and face f , respectively. Since the force applied to the flow by the source term is integrated over the distance (or area/volume in two/three dimensions), the contribution of the source term to the pressure gradient increases with the distance between cell centre and face centre. Consequently, the force applied by the discretised source terms is preserved with distance weighting, given as contrary to the inverse distance weighting in Eq. (24).
Adding the discretised source terms as defined above to Eq. (61), the advecting velocity is defined as with S f obtained by Eq. (76), while the interpolation of ∇p f and S f has to be conducted in the same way, following the explanation given in Section 3.2 for pressure. In order to preserve the low-pass pressure filter of the MWI, the source terms are projected along vector ŝ f in the same way as the pressure term described in Section 3.5. For Eq. (78) to be equivalent to Eq. (70), and for the discretisation to be force-balanced, the discretised source term S P also has to be applied in the discretised momentum equations, see Eq. (4), so that the momentum equations for the quiescent flow discussed above, Eq. (66), become

Numerical experiments
The proposed discretisation of the source terms and the robust pressure-velocity coupling provided by the MWI in flows with source terms is tested for discontinuous source terms in one-dimensional flows on equidistant and non-equidistant meshes and in a two-dimensional flow on a hybrid quadrilateral/triangular mesh, as well as for a spherical drop with surface tension in a three-dimensional domain.

One-dimensional flow
In both cases, the exact solution satisfies and the error of the computed solution is where u in is the velocity at the domain-inlet. The velocity errors ε u (x) for both cases after one time-step are shown in Fig. 26 using different source term treatments; MWI without source terms (abbreviated "MWI"), MWI with source terms S included as-is (similar to [34], abbreviated "MWI-S"), and MWI with source terms S that are discretised as proposed in Section 4.1 (abbreviated "MWI-S "). Both cases demonstrate that a discrepancy between the discretisation of the pressure gradient and the discretisation of the source term causes an artificial acceleration that leads to an incorrect velocity. Using the source term discretisation proposed in Eq. (73), the fluid does not accelerate, matching the exact solution. Away from the discontinuities, all three approaches agree with the exact solution. The ramped case demonstrates that a linear change in pressure gradient, here generated by a linear variation in source term, does not affect the pressure-velocity coupling and low-pass filter of the MWI, as discussed in Section 3.3, while velocity errors ensue at the points where the source term varies non-linearly (x = 0.25 m and x = 0.75 m), if the source term is not discretised consistently.
To test the effect of a non-equidistant mesh and demonstrate the efficacy of the proposed interpolation of the source terms using distance weighting, proposed in Eq. (76), a source term is applied at the central cell P of a one-dimensional domain. The applied non-equidistant mesh has an increasing mesh spacing towards the centre of the domain, schematically shown in Fig. 27. As in the previous case, continuity dictates a constant velocity in the domain, while the source term leads to a change in pressure, with the pressure difference being  As the source term is only applied in cell P , the discretised pressure of the exact solution varies linearly in cell P , but is constant in its neighbour cells, with ∇ p W = ∇ p E = 0. Fig. 28a shows the pressure field computed with distance weighted interpolation and inverse distance weighted interpolation applied to the source term, alongside the exact solution. In both cases the source terms are discretised consistently with the pressure gradients, as proposed in Eq. (73). Consequently, errors in the velocity field are negligible in both cases. However, the applied interpolation has a significant effect on the computed pressure field, as observed in Fig. 28a. The pressure difference is predicted accurately when the source terms are interpolated using distance weighting, as proposed in Eq. (76), while the pressure difference is underpredicted when inverse distance weighting is applied. Varying the cell-size ratio between cell P and its neighbours W and E, the pressure difference remains in excellent agreement with the exact result when the proposed distance weighted interpolation is applied, as observed in Fig. 28b. Applying inverse distance weighting for the interpolation, however, the error in pressure difference compared to the exact result, shown in Fig. 28b, increases linearly with the ratio x P / x F .

Two-dimensional flow
The most attractive aspect of the collocated variable arrangement, and hence MWI, is the ease with which arbitrary meshes can be handled. This means the discretisation of source terms proposed in Section 4.1 must also be applicable on arbitrary meshes. A two-dimensional mesh with dimensions 1 m × 0.2 m consisting of (Cartesian) quadrilateral cells and an embedded region of triangular cells, extending in the range 0.25 m ≤ x ≤ 0.75 m and 0.05 m ≤ y ≤ 0.15 m, shown in Fig. 29, is used to test the source term discretisation on arbitrary meshes. The triangular cells combine the effects of changes in cell size, skewness and non-orthogonality. Using the stepped source term described by Eq. (80), the source term covers the same x-range as the triangular section of the mesh and extends over the complete height of the domain. The flow is introduced with a uniform velocity at the domain-inlet, with free-slip boundary conditions applied to the bottom and top walls. As the source term is constant over the height of the domain, the velocity is expected to remain one-dimensional and uniform. The velocity error ε u along the centreline is shown in Fig. 30. As for the previously discussed test cases, the proposed discretisation and distance weighted interpolation does not yield any noticeable errors. Not including the source terms or including the source terms without taking special care of the discretisation, leads to significant errors in the velocity field. Fig. 30. Velocity errors obtained after one time-step for the stepped source term in a two-dimensional flow, plotted along the centreline of the domain with different source term treatments: MWI without source terms (abbreviated "MWI"), MWI with source terms included as-is (abbreviated "MWI-S"), and MWI with source terms that are discretised as proposed in Section 4.1 (abbreviated "MWI-S ").

Drop with surface tension
The importance of including source terms in the MWI becomes particularly apparent when considering interfacial flows with surface tension. Using the popular Continuum Surface Force (CSF) model [54], the source term representing surface tension is where σ is the surface tension coefficient and κ is the curvature of the fluid interface. The interface indicator function γ is, for instance, given by the Volume-of-Fluid (VOF) method [55], where the local volume fraction γ = 0 and γ = 1 represent the two interacting fluids, respectively, and the interface is present in cells with 0 < γ < 1. Following Eq. (66), surface tension yields a discontinuous change, and hence large gradients, in pressure at a curved fluid interface, for which any imbalance results in substantial errors that may invalidate the simulation results. Although the previous studies of Mencinger and Zun [36] and Denner and van Wachem [37] have addressed this particular issue in detail, the influence of different MWI formulations on the balance between surface tension and the pressure field is briefly revisited here for completeness.
A spherical drop in mechanical equilibrium with initial diameter d 0 = 10 −3 m is considered. The velocity and pressure fields are initialised with u 0 = 0 and p 0 = 0, respectively, and gravity is neglected. To single out the effect of surface tension, both fluids are inviscid and have a density of ρ = 1 kg m −3 . The surface tension coefficient is σ = 0.25 N m −1 and the exact interface curvature κ = 4000 m −1 is applied at the interface. The drop is situated at the centre of a cubical domain of size 2d 0 , which is represented by an equidistant Cartesian mesh with mesh spacing x = d 0 /25. The applied time-step t satisfies the static capillary time-step constraint [56]. Because the drop is in mechanical equilibrium, the pressure difference as a result of surface tension is given by the Young-Laplace law [57], given as p = p in − p out = σ κ = 1000 Pa, where p in and p out = p 0 refer to the pressure inside and outside the drop, respectively. Furthermore, the flow should remain quiescent, since the only force acting on the fluid (the force due to surface tension) is balanced by the pressure field. A compressive VOF method [58] is used to advect the interface with the underlying flow. Fig. 31 shows the normalised pressure error and the velocity magnitude |u|. The pressure error is negligible, with ε p < 10 −13 , if the source term representing the surface tension is accounted for in the MWI. However, not including the source term in the MWI leads to large pressure errors of up to ε p ≈ 50%. This pressure error translates into spurious currents, see Fig. 31b, which are absent if the source term is included correctly in the MWI.

Density discontinuities
In addition to varying or discontinuous source terms, a feature of many multiphase flows are discontinuous density fields, which results in discontinuous pressure gradients and leads, in turn, to oscillatory solutions or failure to reach a solution [37]. Rearranging the momentum equations, Eq. (1), as an expression for the pressure gradient, shows that the pressure gradient is proportional to the density. Due to the use of linear interpolation in the discretisation of the pressure field, this discontinuity cannot be represented by the discrete representation of the pressure gradient, leading to discrepancies between the discretised momentum equations and the equation for the advecting velocity at the face.
To illustrate the problem that arises, a one-dimensional, inviscid, incompressible two-phase flow, with the two bulk phases separated by a sharp interface, subject to a constant acceleration in the absence of source terms is considered. Because both fluids are incompressible, the velocity field should be spatially uniform. Hence, the discretised momentum equations reduce in this case to where a is the spatially uniform acceleration of the flow. Similar to Eq. (87), the discretised momentum equations, Eq. (88), shows that the pressure gradient is proportional to the cell-averaged density and is, therefore, discontinuous where the density is discontinuous. In the limit of large density ratios, the pressure gradient in the heavier fluid is significantly larger than in the lighter fluid, i.e.
and the pressure term in the MWI of the advecting velocity, Eq. (61), is where l f is the interpolation coefficient. As a result, the discrete pressure gradient is underpredicted in the heavier phase and overpredicted in the lighter phase, which leads to an artificial acceleration of the flow in the vicinity of the interface. In the case of extremely large density ratios, the large and unphysical force applied to the lighter phase may lead to divergence of the solution algorithm [37].

Density weighting in the MWI
Denner and van Wachem [37] proposed to weight the pressure gradients in the MWI by the corresponding density, with the pressure term in the MWI of the advecting velocity, Eq. (61), becoming where the cell-centred pressure gradients are interpolated with 1/2-weighting and with the face density evaluated by harmonic averaging ρ f = 2/(ρ −1 P + ρ −1 F ). The generalisation of this density weighting is straightforward, with the pressure term of the MWI given as where l f is the interpolation coefficient, and the face density defined by Note that the interpolation coefficients of the cell-centred pressure gradients add up to unity, which is crucial for a consistent and bounded interpolation. In the limit of large density ratios, for instance with ρ F /ρ P → ∞, the face density becomes and the pressure term in Eq. (92) follows as Thus, the cell-to-cell change in pressure tends to the value corresponding to the minimum pressure gradient on either side of the density discontinuity. This has a stabilising effect in the discretisation, because instead of applying a force that is too large in the lighter fluid, leading to large accelerations, a force that is too small is applied to the heavier fluid. Consequently, the errors associated with a discontinuous change in density are substantially reduced and the numerical solution remains stable.

Application to two-phase flows
An accelerating incompressible two-phase flow in a one-dimensional domain is simulated, in which the bulk phases with different densities are separated by a sharp interface. The flow is initially quiescent and the lighter fluid occupies the entire domain. The heavier fluid is introduced at the inlet with velocity where a is a constant acceleration and t is time, and the discretised momentum equation is given by Eq. (88). In theory, the solution is a spatially uniform, time-varying velocity given by Eq. (97), with a discontinuity in the pressure gradient field at the interface due to the discontinuity in density. As considered in Section 4.2.3, the two-phase flow is represented using the VOF method [55] and a compressive VOF method [58] is used to advect the interface with the underlying flow.
where the value of ρ(x) is determined by the interface location x (t) = at 2 /2 as Fig. 32a shows the pressure gradients of the numerical solution obtained with and without density weighting for a two-phase flow with a density ratio of ρ H /ρ L = 50. Both cases exhibit errors in the computed pressure gradient. The exact solution is two constant gradients with a discontinuity at the interface, whereas over-and undershoots are evident in the numerical results. However, with the density weighting in the MWI the maximum magnitude of the error is <10% compared to the exact value, whereas the maximum magnitude of the error increases to more than 50% without the density weighting in the MWI. The ensuing errors in the velocity field, shown in Fig. 32b, confirm that the errors are diminished when the density weighting is applied in the MWI. Fig. 33 shows the same case, but for two fluids with density ratios ρ H /ρ L = 10 3 and ρ H /ρ L = 10 6 . Only results obtained with the density-weighted MWI are shown in Fig. 33, because a converged numerical solution without density weighting is not available for these cases. Despite the large density ratios, the pressure gradients are predicted accurately when the density weighting is applied in the MWI, with only small errors. The presented results, notably a comparison of Figs. 33a and 33b, suggest that the magnitude of the error associated with the pressure gradient increases with increasing density ratio. In fact, as seen in Fig. 34, the error magnitude of the computed pressure gradient increases linearly with density ratio. Hence, the relative error in pressure introduced by the density-weighted MWI is independent of the density ratio. Furthermore, the presented results suggest that the error of the computed pressure gradient is largely contained in the heavier phase, which is desirable as it reduces spurious accelerations of the flow, and, hence, in conjunction with the linear relationship between density ratio and associated errors, indicates that the density-weighted MWI ensures a stable result for a wide range of density ratios. To this end, Denner and van Wachem [37] demonstrated a significant improvement of convergence for two-phase flows with density ratios of up to ρ H /ρ L = 10 9 using this density weighting. Furthermore, Denner and van Wachem [56] reported stable results for a two-phase flow with a density ratio of ρ H /ρ L = 10 24 , without any noticeable errors, far exceeding density ratios of typical multiphase flows.   This formulation provides a robust pressure-velocity coupling and a low-pass filter on third and higher derivatives of the driving pressure on arbitrary meshes, it is time-step independent and satisfies the steady-state solutions of steady-state ( t → ∞) as well as initially transient ( t is finite) problems. Furthermore, this formulation provides stable results for any density ratio and it reduces to the normal velocity at cell faces, ϑ f = u fn f , for steady-state solutions with constant or linearly changing driving pressure.

Conclusions
When simulating flows in and around complex geometries, the discretisation of the governing equations is greatly simplified by using a collocated variable arrangement. In simulations of incompressible and low Mach number flows this gives rise to pressure-velocity decoupling, with the characteristic chequerboard pressure field, if a straightforward discretisation is employed [1]. The momentum-weighted interpolation (MWI), typically attributed to have been introduced by Rhie and Chow [18], is a widely used method to couple pressure and velocity in collocated variable arrangements and as a remedy for pressure-velocity decoupling. However, in the current literature there are a number of varieties of MWI, and it is so far unclear what the optimal formulation is.
In this paper, a unified formulation of the MWI for arbitrary meshes has been derived based on physically-consistent arguments, including extensions for discontinuous source terms and discontinuous changes in density. The presented stepby-step derivation and analysis of the MWI has been used to develop theoretical justifications for the discretisation of velocity, pressure and source terms, including the applied interpolation and weighting coefficients, under the main assumption that MWI enforces a low-pass filter acting on the discrete pressure field, thereby imposing a direct relationship between neighbouring pressure values that suppresses oscillatory solutions. This theoretical analysis has been further supported with numerical results of representative test cases on arbitrary (structured and unstructured) meshes, demonstrating the impact of the MWI in general as well as the impact of the low-pass filter enforced by the MWI. The proposed MWI formulation has been demonstrated to yield a third-order convergence under mesh refinement with respect to the conservation error of kinetic energy. Although the conservation of kinetic energy has been shown to be time-step independent for MWI formulations derived from the steady and the transient momentum equations, the transient term of the MWI has been found to be essential for an accurate prediction of transient pressure perturbations.
With regards to discontinuous or strongly varying source terms, only the driving pressure gradient, i.e. the pressure gradient associated with the flow, should be coupled to the velocity field. Failing to account for source terms in the MWI can lead to decoupled solutions and artificial accelerations of the fluid in cases which include sharp gradients and discontinuities of source terms, as demonstrated in the presented results. The proposed reconstruction of the discrete source term provides an exact balance with the discretised pressure gradient, a so-called force-balanced discretisation, and conserves the force applied to the flow, on arbitrary meshes. Furthermore, the application of a density weighting in the MWI has been analysed from a theoretical perspective and shown to have a stabilising effect on flows with large density ratios. Without such a treatment, the fluid would be accelerated without limit, resulting in divergence of the solution algorithm. MWI can also play a vital role in simulating low Mach number flows using pressure-based algorithms to overcome the weak pressure-density coupling when the compressibility of the flow is negligible, as demonstrated by the presented results.
In summary, MWI is very effective in maintaining pressure-velocity coupling in simulations of incompressible and low Mach number flows on meshes with collocated variable arrangement, but the effect of external forces, and of the discrete approximations itself, on the discretised pressure gradient have to be accounted for carefully to obtain physically realistic results and robust solutions. In all considered cases, the proposed MWI has been shown to offer superior accuracy and stability compared to the considered alternatives, in particular with regards to meshes with large non-orthogonality and in flows that are subject to discontinuous source terms or large density differences.