A Consistent and Implicit Rhie–Chow Interpolation for Drag Forces in Coupled Multiphase Solvers

: The use of coupled algorithms for single ﬂuid ﬂow simulation has proven its superiority as opposed to segregated algorithms, especially in terms of robustness and performance. In this paper, the coupled approach is extended for the simulation of multi-ﬂuid ﬂows, using a collocated and pressure-based ﬁnite volume discretization technique with a Eulerian–Eulerian model. In this context a key ingredient in this method is extending the Rhie–Chow interpolation technique to account for the unique ﬂow coupling that arises from inter-phase drag. The treatment of this inter-ﬂuid coupling and the fashion in which it interacts with the velocity-pressure solution algorithm is presented in detail and its effect on robustness and accuracy is demonstrated using 2D dilute gas–solid ﬂow test case. The results achieved with this technique show substantial improvement in accuracy and performance when compared to a leading commercial code for a transonic nozzle conﬁguration.


Introduction
The complexity of phenomena in multiphase flows is reflected by the variety of numerical methods used for their simulation: for free surface flows interface tracking methods, such as VOF [1][2][3], level-set [4,5] or front-tracking methods [6,7] provide an accurate representation of the interface that separates the phases. The Euler-Lagrange approach is used for more complex flows with multiple, complex interfaces and when a detailed analysis of local processes is needed, such as for the analysis of nucleation in wet steam [8,9]. However, the Euler-Lagrangian formulation suffers from many limitations when scaling to large industrial applications with complex geometries [10][11][12]. For such flows the Euler-Euler approach, where a set of averaged conservation equations are solved for all continuous and discrete phases, along with a series of closure models [11,13], is the prefered approach. The main algorithm for these types of simulation has been some variant of the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm [14] extended to solve multiphase flows through the addition of a supporting algorithm that deals with the inter-fluid terms known as the Inter Phase Slip Algorithm (IPSA) [15,16]. This approach has allowed for the simulation of complex flows; however, it suffers from low robustness and convergence problems [17]. This degradation of performance is easy to understand since in multiphase flows, the inter-fluid coupling and the phasic velocitypressure coupling have to be resolved for convergence, with an algorithm that uses a semiimplicit approach to resolve the velocity-pressure coupling of one phase at a time. While some remedies have been proposed to address the inter-fluid coupling as in Spalding's Partial Elimination Algorithm (PEA) algorithm [15,16] or Lo's SImultaneous solution of Non-linearly Coupled Equations (SINCE) algorithm [18], the solution of multiphase flows is still much more expensive than is warranted by the increase in the number of equations.
In a number of articles [19,20] the authors have presented a fully coupled approach to resolve the velocity-pressure coupling that arises in single fluid flow. In this approach, the momentum and continuity equations are assembled and solved as one block resulting in a substantial increase in performance and robustness. The coupling of velocity and pressure variables in one system of equations presents many opportunities when solving multiphase flows. For example the inter-fluid drag can be discretized in a fully implicit fashion as in [13]. However, by neglecting the phasic velocity pressure coupling, the approach of Kunz falls short in terms of performance [17].
In this paper the authors propose to treat all the main couplings present in multiphase flows, that is the inter-fluid coupling and the phasic velocity-pressure coupling, in an implicit fashion by treating the pressure and the phasic velocity fields as one coupled system of equations. At the core of this method is the treatment of the pressure equation that arises from the phasic mass conservation equations; for this end, a novel, implicitly consistent momentum interpolation technique is used based on a fully coupled formulation of the momentum equations. This formulation not only removes the partially explicit treatment of the drag forces, but it is derived in such a manner that simplifies its implementation in existing numerical frameworks for the simulation of multiphase flows.

Governing Equations
The governing conservation equations are summarized here for later reference. Index k refers to the phase index. Since this article focuses on two phase flows, the two phases will later be named α and β.

Phasic Energy Equation
Each component of the single fluid Navier-Stokes equations is multiplied with a prefactor r, known as volume fraction, which accounts for the volume occupied by each phase. Interphase source terms are summarized in the expressions S k c , S k m and S k e along with any other source terms such as, e.g. buoyancy.

Traditional Segregated Multiphase Algorithms
As for single phase flows, the majority of publications in the field of numerical analysis is based on segregated solution techniques. These algorithms solve the governing equations sequentially. This solution strategy has several advantages, e.g., the ease of implementation. However, it is mainly driven by low memory requirement, which was the bottleneck of computational fluid dynamics in its early days. The segregated treatment, however, decouples the velocity from the pressure during inner iterations. A numerical decoupling of this physically strong connection results in poor convergence behavior, an effect which increases for large scale problems. The authors have therefore developed a framework that simultaneously solves the pressure and momentum equation, i.e. combines the three momentum and the pressure equation in a coupled system of linearized equations in a block-matrix structure. Details on the coupled framework for single phase flows can be found in [19,20].
For multiphase flows, the segregated solution strategy is not only problematic due to the pressure velocity decoupling. Similar issues arise through a variety of interfacial source terms which couple the various phases. To introduce the coupling terms, the linearized momentum equation is given in Equation (5) for phase α.
Term A accounts for momentum transfer due to mass exchange between phase α and β and term B originates from interfacial drag. For segregated algorithms, discretization of these terms is limited to a partially implicit formulation, keeping coupling to the velocities of other phases fully explicit on the right hand side. Rearranging any implicit dependency on the phasic velocity on the left hand side leads to Equation (6) [21]. The superscript * refers to previous iteration values.
From term B in Equation (5), it can be seen that the drag coefficient multiplies the slip velocity. Numerical problems arise, when the slip velocity is small. The drag coefficient is increasing exponentially, towards low slip velocities.
For two-phase flows, this issue was addressed by the Partial Elimination Algorithm (PEA) [16], which is briefly explained here. The momentum equation for the two phases α and β are written as given in Equation (5). In a next step, the PEA algorithms forms a new equation through a summation of Equation (5). This result is then used to find expressions for u α p,i and u β p,i . Inserting these expressions in the two original momentum equations as given in Equation (5) eliminates the interfacial drag term and increases the diagonal dominance of the equation system. Improved stability of the numerical solution procedure is thus achieved. The PEA algorithm is, however, only applicable to two phase flow situations. An alternative formulation is given through the SImultaneous solution of Non-linearly Coupled Equations (SINCE) [18]. This algorithm requires to solve a N p xN p equation system for each cell with N p being the number of phases. For both algorithms, PEA and SINCE, it needs to be mentioned, however, that the interfacial mass transfer term is still treated semi implicitly. Detailed overview concerning the efficiency of different interface coupling algorithms is given in [22]. An extension of the coupled solution strategy to multiphase flow problems would abandon the need for complicated stabilization problems originating from semi implicit treatment of the interfacial source terms. The development of such a solution strategy is therefore explained in what follows.

Novel Coupled Multiphase Framework
As stated above, segregated solution techniques solve the governing equations in a sequential form. This prevents a thorough coupling of the variables in the solution vector of the linearized system of equations. A pressure-based two phase Euler-Euler system needs to solve three momentum equations for each of the two phases and a pressure equation, hence seven equations are solved. Considering the standard segregated solution method, after linearization, each of these equations can be written in scalar form as given in Equation (7).
This scalar variables of x are the components of the velocity of phase α (u α x , u α y , u α z ), velocity of phase β (u β x , u β y , u β z ) and the pressure (p). Hence, the only implicit contributions are terms depending on the current scalar solution variable, any coupling to other solution variables is explicitly accounted through the source b c . The proposed formulation solves this drawback by extending the scalar assembly to a matrix form as shown in Equation (8).
The coefficient matrix A, the solution vector x and the right hand side b are written as: This work can therefore be understood as an extension to what is described in [19,20] to multiphase flows. Compared to the other solution methods for segregated algorithms, this 7 × 7 block matrix structure allows for a fully implicit implementation of interfacial source terms such as the drag. This completely abandons the need for complicated stabilization methods as described above.

Coupled Rans-Equation Assembly for Two Phase Flows
This section summarizes the discretization of the Navier-Stokes equations as used in the presented framework. The single phase version is already presented in [19,20,23]. Extension to real gas state equations is given in [24]. The following lines therefore highlight the differences when moving to multiphase flow analysis.

Momentum Equation, Fully Coupled Drag
The formulation of multiphase momentum equations was given in Equation (3). Extending the inter phase source term S k m to account for drag allows to rewrite the equation for the two phase α and β as given in Equations (10) and (11), here assuming steady state.
With the drag coefficient usually taking the form given in Equation (12). The coefficient c d is then modeled using well known drag models such as the one of Schiller-Naumann [25] with the slip velocity u s,i defined as u Using the introduced coupled 7 × 7 block matrix structure, the drag can be treated fully implicit as shown in Equation (13). For traditional sequential methods, the coefficients highlighted in blue contribute explicitly to the right hand side. The proposed coupled approach therefore improves the convergence behavior of the iterative procedure.

Continuity Equation
The basis of the pressure-based continuity equation starts with the definition for the conservation of global mass. Summation of the phasic equations leads to Equation (14) with its semi-discrete form in Equation (15).
A detailed derivation of the pressure equation based on Equation (15) is provided in [24] for real gas flow in single phase analysis. The main steps are repeated here to maintain readability. The density uses a first order linearization, introducing the implicit dependency on pressure (Equation (16)). For the non-linear convective part, a Newton linearization is used (Equation (17)). Finally, the actual time step velocity is expressed using Rhie-Chow interpolation (Equation (18)) with D k i,j being the inverted coefficient matrix of the momentum equation. A detailed discussion on this special interpolation is given in a later section, including the derivation of a novel fully implicit formulation for multiphase flows.
The superscript n stands for the actual time step index, superscript * indexes values from previous inner iterations. The notation Φ c f describes an interpolation of the variable Φ from point c to point f . Using the above given expressions, the multiphase mass conservation equation in a pressure-based form is found as given in Equation (19).
Considering the coefficient matrix given in Equation (9), the second part of term P1 and term P4, contributes to coefficient a pp . These coefficients are also present for segregated algorithms and therefore again colored red. The improved robustness of the suggested algorithm comes from term P2. This term links the pressure equation to the velocity field, setting this algorithm apart from the segregated approach. The contribution is therefore highlighted in blue in Equation (20). For segregated algorithms these terms would have been assigned explicitly to the right hand side, explaining the superiority of the chosen coupled approach.

Momentum Interpolation Techniques
The Rhie-Chow interpolation [26] is critical to the derivation of the pressure equation when using a collocated variable arrangement. It ensures that checkerboarding, unphysical oscillations are suppressed by defining an interpolation of the momentum equation that retains the strong physical dependency between the face velocity and the pressure field at the cell centers. In effect the Rhie-Chow interpolation, also denoted by Momentum Interpolation, mimics the staggered variable arrangement. After a brief introduction of the main features of the Rhie-Chow interpolation, its acurrent application to multiphase flow simulation is shown to be inconsistent, and a novel, fully consistent momentum interpolation technique is presented that completely suppress checkerboarding, unphysical oscillations from multiphase flow simulations.

Standard Momentum Interpolation
The Rhie-Chow interpolation is used to interpolate the velocity field to the cell face in a way that imitates a staggered variable arrangement. Starting with the momentum Equation (21) we can define a pseudo momentum equation defined at the face centers (22). Any remaining terms are summarized in H as given in Equation (23).
where the value A f , V f and H f are interpolated from the cell values expressed as Φ c f Equation (24). This compares with the simply interpolated momentum equation Equation (25).
Subtracting Equation (24) from Equation (25) leads to an expression providing a momentum-consistent representation of the velocity on the face, see Equation (26).
Considering the interpolation and multiplication to be commutative, the following assumptions are made.
In the following, the inverted coefficient matrix A c −1 f will be labeled D f . The final expression for the face velocity therefore reads as: The Rhie-Chow interpolation was extended over the years to additional terms in the momentum equation such as the relaxation factor [3,27], the time step size [28] and body forces [29]. An approach was presented by [30,31] to overcome the problems of relaxation factors and time step dependence.
For multiphase flows, it is the interfacial drag that has to be accounted for, especially since it can become the dominant term [32]. The treatment is best understood by reviewing the definition of the discretized momentum equation in Equations (21) and (22). Any term that is only represented in H is removed in the original derivation since it is assumed that H f = H c f . Any dominant force, has to be removed from this term and handled separately. Two methods to achieve this are presented in what follows.

Standard Decoupled Multiphase Momentum Interpolation
For multiphase flows, the main body force is governed by the drag between the two phases, neglecting for the moment possible mass exchange between the phases. The drag term in the momentum equations usually takes the form given in the last term of Equations (30) and (31).
The standard methods known from literature start by moving the contribution to the actual phase velocity from the drag to the left hand side [32,33], leading to the semi-implicit formulation given in Equations (32) and (33) Again, starting with the cell-based and a pseudo faces-based momentum equation, the drag term can be expressed as Equations (34) and (35).
This formulation presents a consistent extension of the momentum interpolation technique to multiphase flows. However, the explicit contribution on the right hand side can still present problems when drag is dominant. This issue is addressed in the next fully implicit formulation.

Proposed Coupled Multiphase Momentum Interpolation
In this section, a novel, fully-coupled formulation of the Rhie-Chow interpolation for multiphase flows is presented. The derivation starts with the coupled formulation of the momentum equations. For segregated algorithms, the coupling to the second phase is prohibited as shown in Equation (36). It must thus be accounted for explicitly and its contribution is kept constant and depending on previous iteration values. The coupled formulation, shown in Equation (37), removes this drawback based on its capability to store cross coupling terms in the coefficient matrix.

Segregated momentum:
A α A procedure similar to the original Rhie-Chow interpolation is applied. The staggered approach is mimicked by expressing the set of momentum equations on the faces, as in Equation (38), while the interpolation from the cell-centered equation is written as Equation (39).
As for the standard Rhie-Chow interpolation, Equation (38) is subtracted from Equation (39). With some algebraic manipulations, the face velocities is be written as Equation (40).
The face velocities can now be re-written individually for each phase as: The mesh consists of 10'000 cells and has a total length of 20 meters indicated as L in Fig. 1. The   153 considered setup is what is known as "Dilute Gas-Solid Flow" from [36]. At the inlet, the continuous   Figure 2a shows the dispersed phase velocity from inlet to outlet. The squares show the analytical solution, with the fine dotted line being the results obtained with the developed framework. With the numerical results being on top of the analytical solution, it can be concluded that the proposed approach is consistent in that respect. A quick assessment of the convergence behavior for this simple case is given in Figure 2b. The graph includes results published in [36], demonstrating the capabilities of the proposed approach. A smooth convergence is achieved with both phases showing similar convergence behavior.

Analytical 2D Case for the Assessment of the Developed Rhie-Chow Formulation
The key parameter of any momentum interpolation technique is its capability to suppress unphysical checkerboarding fields. Therefore, the following testcase provides a very simple geometry that provides an example of the quality of the achieved solution when compared to the standard Rhie-Chow formulation. The geometry is given in Figure 3. The key parameter of any momentum interpolation technique is its capability to sup nphysical checkerboarding fields. Therefore, the following testcase provides a very simple geo hat provides an example of the quality of the achieved solution when compared to the sta hie-Chow formulation. The geometry is given in Fig. 3.  Two identical fluids are entering the computational domain with the only difference being the inlet velocity. The continuous fluid used an inlet velocity of 10 m/s and the dispersed phase is assigned a velocity of 9.99 m/s. With both fluids being computed inviscid and laminar to reduce the uncertainty in the source for any spurious oscillations, the only force that could lead to an inconsistency in the formulation of the momentum interpolation is the drag itself. This force is modeled using Schiller-Naumann [25] using a hard coded viscosity value of µ = 1.83e −7 kg ms . The solutions of the volume fraction fields are presented in Figure 4. In the image to the left, the solution using standard Rhie-Chow interpolation is presented. As can be seen from the image to the right, the developed algorithm is ably to suppress any unphysical oscillations.

Two-Phase Flow in a Transonic Nozzle Configuration
Validation against commercial code results is carried out based on a geometry published by the Jet Propulsion Laboratory (JPL) located at the California Institute of Technology [37]. The flow is a dilute two-phase flow in a converging-diverging nozzle. Two-phase flow investigations have been published by [38][39][40]. The overall geometry is given in Figure 5 and is discretized using 4655 structured hexahedra elements. At the inlet, the total conditions are set to p c 0 = 10.34 bar and T c 0 = 555 K for the continuous phase. The dispersed particle phase uses the same static pressure at the inlet as the continuous gas phase with a mass fraction of Ψ d = 0.3. At the outlet, supersonic conditions are applied, accelerating the dilute two-phase flow to transonic conditions through the nozzle. The energy equation is solved for the continuous gas phase using ideal gas state equation, while the particles are assumed isothermal and incompressible with a density of ρ d = 4004.62 kg/m 3 . The two fluids are assumed inviscid with the viscosity in the drag term being defined using Sutherland's law as given in Equation (51).
As shown in [39], two different particle sizes have been chosen. The first particles have a radius of 1 µm and the second simulation uses particles of r p = 10 µm. Figure 6 shows the results for the smaller particles and Figure 7 the results with the bigger particles. Both results compare well with commercial code results. Unphysical oscillations as visible at the symmetry plane of the commercial code results are not present for the proposed method. As expected, the lighter small particles follow the strong curvature of the convergingdiverging nozzle section, the heavier particles tend to generate a particle free zone in the diverging section of the nozzle.   . Particle radius r p = 10 µm, r p ∈ 7.2 · 10 −5 · · · 0.001 .
A comparison of the root mean square (RMS) error of the residual for the case of r p = 1 µm is shown in Figure 8. Both simulations are started with uniform initialization, upwind advection scheme, a time step of ∆t = 1e −3 s and have been run using double precision. Double precision is requested from the commercial code whenever multiphase calculations are performed. This originates from the fact, that the volume fraction can attain very small values but still be of importance to the numerical results. While both codes initially close the outlet due to flow reversal occurring from the non-initialized conditions, the in-house framework recovers quickly from this situation. The problematic convergence behavior of commercial code up to iteration 600 is related to this closed faces at the outlet. However, after having recovered, it takes another 600 iterations to arrive at the supersonic conditions. Only thereafter, the simulation starts to convergence quickly. The coupled framework was introduced in Section 4. A method is proposed that solves the two sets of momentum equations and the continuity equation in a 7 × 7 block matrix structure. It was later explained in Section 5.1 how this framework allows to treat inter phase drag implicitly in the coefficient matrix. The coupled drag formulation is therefore compared to traditional semi-implicit, segregated treatment of the drag term. In order to to so, the coefficients marked in blue in Equation (13) are entirely moved to the right hand side of the linearized equation system. The comparison therefore only differs in the treatment of the drag term, any other coupling is kept and the setup of the test case is identical. The results shown in Figure 9 clearly demonstrate the advantages of the proposed solutions. The RMS-residuals of the implicit formulation shown in the figure to the right outperform the explicit formulation.

Conclusions
A new formulation of the classical Rhie-Chow interpolation method was presented. The proposed solution takes advantage of a coupled formulation of multiphase momentum equations. The new formulation accounts implicitly and consistently for the drag body forces. In addition to yielding improved convergence and more accurate results, this coupled approach can be readily implemented in existing codes for multiphase flow simulation. The main modifications is in the computation of the inverted coefficient matrix, while computational memory requirements remain unchanged. The algorithm is first compared against analytical solution, showing the consistency of the implementation. The second testcase demonstrates the benefits of the consistent, coupled formulation compared to the original Rhie-Chow formulation. Finally, a transonic nozzle configuration was chosen and results were compared against commercial code to demonstrate the quality of the framework. The coupled formulation was then compared against the traditional segregated formulation of the drag term. Superior convergence behavior of the proposed approach was demonstrated. With the current derivations being limited here to two-phase flows, extension to any number of fluids is straightforward.
Author Contributions: The extension of a pressure-based coupled numerical algorithm to multiphase flow physics was mainly driven by L.H. The final target of the thesis is a framework, able at predicting two-phase flow physics for steam turbines. It is, however, based on the excessive work done by L.M. and M.D. in the development of such a coupled framework and thus also supported by these two authors. E.C. and D.M.V. provided a variety of test cases and material for the validation of the presented algorithms and supervised the progress. All authors contributed to writing and editing this paper. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors gratefully acknowledge the financial contribution provided by the Swiss National Science Foundation (SNF: Grant-Number 175900).