Basic verification of a numerical framework applied to a morphology adaptive multifield two‐fluid model considering bubble motions

A morphology adaptive modeling framework is derived that is able to handle computationally efficiently dispersed as well as resolved interfacial structures coexisting in the computational domain with the same set of equations. The Eulerian multifield two‐fluid model is combined with the compact momentum interpolation method for multiple phases, which has been proposed in the literature as an extension to the Rhie‐Chow pressure‐velocity coupling. Additionally to the interfacial drag force, the virtual mass force is consistently accounted for in the model. Utilizing a specialized interfacial drag formulation, large interfacial structures can be described with the presented method in a volume‐of‐fluid‐like manner, additionally to the disperse description. The strong phase coupling due to the drag closure model in interfacial regions is resolved with a partial elimination algorithm, which is adapted to work in an approximate manner for more than two phases via a sum formulation. The presented model is implemented in the C++ library OpenFOAM and solver performance is compared with results obtained with the homogeneous model approach in two cases of a single rising gas bubble for two‐ and three‐dimensional space, respectively. Additionally, for both three‐dimensional cases, the results are compared with experimental data. Finally, the presented method's capability of representing dispersed and resolved interfacial structures at the same time is demonstrated with two test cases: a two‐dimensional gas bubble, rising in a liquid, which is laden with micro gas bubbles, and a two‐dimensional stagnant stratification of water and oil, sharing a large‐scale interface, which is penetrated by micro gas bubbles.

kinds of reactors. In such fields, reliability and safety are of great interest. Flows occurring in the aforementioned applications are characterized by a high level of complexity. Multiple different morphology regimes appear, even in facilities of simplest geometries, such as pipe flows, where segregated, bubbly, annular, or slug flows play a role. The level of complexity is further increasing, if countercurrent flows or wave breaking phenomena are of interest. From that arises a high demand for simulation techniques for predicting gas-liquid flow behaviors with different flow morphologies sufficiently exact, with feasible computational expense and less knowledge about the occurring flow regimes prior to the investigation. One of the biggest challenges in describing multiphase flows is, in analogy to turbulent flows, the wide range of spatial scales, which span a range between the molecular length scale for the thickness of an interface to the dimensions of an industrial facility. There are many situations, in which a large-scale interface (a boundary between two immiscible phases), for example, a free surface, meets disperse interfacial structures, such as bubbles or droplets, or situations, in which disperse interfacial structures of different phases, size or shape interact with each other. For each of those individual flow morphologies, there are multiple modeling approaches to choose from, which are characterized by different levels of accuracy and computational efficiency.
Interfaces, which exist at a scale above grid scale, are resolved on the computational grid. Wörner 1 gave a comprehensive overview over numerous different interface-resolving models that have been proposed in order to capture the physics of large-scale interfaces, for example, Lagrangian methods such as marker and cell 2 or front-tracking 3 as well as Eulerian approaches, for example, level set, 4 conservative level set, 5 volume of fluid (VOF) with 6 or without 7 interface reconstruction, or phase field methods. 8 All of these methods are typically formulated as a homogeneous model, which means that immiscible phases share a common velocity field and physical properties are considered as mixing properties depending on phase fractions.
Interfacial structures, which are smaller than grid scale-the disperse phase-cannot be resolved and have to be described by appropriate models instead. One way to achieve this is the application of an Euler-Lagrange method, where interfacial structures, such as bubbles or droplets, are individually described by point elements and get transported with an equation of motion, which considers all body forces as modeled source terms. 9 Another way is to represent continuous and disperse phases as interpenetrating continua-the Euler-Euler point of view-in an ensemble averaged manner, 10 where interfacial closure models are formulated as continuous forces and include assumptions about element size, number density, and shape. Interfacial forces are accounted for by including volume force terms into the momentum equation.
It is desirable to describe both interfacial flow structures, which are smaller than computational grid size, as well as such, which are larger, in a single-computational domain. Therefore, interface resolving and nonresolving models have to be combined in a hybrid model approach and the appropriate individual models are applied, depending on the local flow morphology. Lagrangian point particle methods have been combined with a level-set method, 11,12 as well as with a VOF method. 13 Additionally, there is a large number of attempts to combine Euler-Euler models for dispersed flows with VOF or VOF-like models with or without interface reconstruction (see Bestion et al 14 for a comprehensive classification). VOF-like methods with interface reconstruction and disperse Euler-Euler model were blended in a two-fluid formulation by Boualouache et al 15 and Yan and Che, 16 and with switching between one-and two-fluid formulation byČerne et al. 17 Lopes et al 18 added an additional disperse Eulerian phase to a VOF model without interface reconstruction with all continuous phases sharing a single-velocity field. Special attention was paid to stratified flows and large interfaces with modeling of interface features below grid scale. This has been realized with VOF-like methods without reconstruction in two-fluid formulation. [19][20][21][22][23][24] Contrary to that, a VOF-like approach without interface reconstruction was combined with Euler-Euler modeling for dispersed flows strictly limited to two fields, representing both physical phases, for example, a gas and a liquid phase. [25][26][27][28][29][30] Within the frame of this modeling approach a blending method ensures that the two-fluid model behaves either as a VOF-like model or as disperse Euler-Euler model, depending on the local flow morphology. To achieve this, interfacial momentum exchange terms have to be formulated for all possible regimes, namely continuous liquid with dispersed gas, continuous gas with dispersed liquid, and continuous gas with continuous liquid. The latter combination describes the interaction between both phases in the region of a pronounced, shared interface. For instance, Wardle and Weller 30 combined an artificial interface compression term for phase fraction transport with a residual drag formulation in order to prevent resolved interfaces from smearing by numerical diffusion. With this method, a slip velocity between both phases sharing an interface is still allowed in that region. Finally, model contributions of all three regimes are summed up, weighted by a regime indicator function, which is based on the ratio of gas and liquid volume fraction. This kind of formulation offers a compact solution for a number of different flow regimes but also implies limitations. For instance, interactions between multiple disperse phases or between a disperse phase and two different continuous phases cannot be covered by such a hybrid model, because only two phases with individual velocity fields are considered. To account for more complex flow situations, disperse phases need to be simulated as individual numerical phase, additionally to the numerical phases for regions, where continuous flows occur. One possible combination is continuous water, continuous air, air bubbles, modeled as disperse air, and water droplets, modeled as disperse water, in the same computational domain. Tomiyama and Shimada 31 and Tomaselli and Christensen 32 proposed models, where both gas and liquid continuous phases share a single-velocity field and every disperse phase is considered to have its own velocity field. Such a method is able to handle an arbitrary number of disperse phases, but only one pair of gas and liquid phases. In contrast to that, Mer et al 33 38 paid special attention to an accurate description of resolved interfaces with a two-fluid model, whereby both authors used conservative level-set methods. As stated by Štrubelj et al, 38 from a mathematical point of view it is convenient to use a two-fluid formulation for both VOF-like and disperse Euler-Euler models in order to solve for a unique set of equations in the whole computational domain. Furthermore, the choice of a two fluid model in principle allows for a nonzero slip velocity between two continuous phases at an interface. Gauss et al 39 showed that such an interfacial slip allows for the description of large interface structures, such that the rising velocity of a single bubble can be predicted correctly, even on numerical grids, which are to coarse to fully resolve physical processes on all scales.
The present work takes the model of Hänsch et al 35 as a conceptual starting point with the goal to describe all continuous and disperse phases as individual numerical phases in a computationally efficient way. A numerical framework is presented, which incorporates a multifield two-fluid model, which is able to account for resolved interfaces. In principle, the two-fluid model allows slip between two phases sharing an interface, for example, in the context of underresolved interfaces. However, within the scope of the present work, a no-slip condition between continuous phases is preserved. The solution procedure is implemented in the OpenFOAM C++ library and the source code is made public by Meller et al (2020). 40 To achieve numerical robustness and sufficient accuracy, the compact momentum interpolation method (CMI) for multiphase flows, which was proposed by Cubero et al, 41 is extended to consistently take into account the interfacial momentum exchange closure of virtual mass force.
The interfacial drag force of Štrubelj and Tiselj 42 is applied to pairs of continuous phases, to describe resolved interfaces. The strong phase coupling resulting from that is approximately resolved with a partial elimination algorithm (PEA), 43 which is formulated for an arbitrary number of phases. In order to show capability of the presented solvers to describe the interplay between large resolved interfaces and finely dispersed flow regimes, interfacial closure models for disperse flows are included, namely drag and virtual mass forces.
Physical modeling and the numerical method are described in Sections 2 and 3, respectively. Subsequently, in Section 4, the capability of the presented method is demonstrated with for a selected set of test cases considering bubble flows. Finally, results are concluded and perspectives are shown in Section 5.

MODELING
The presented model framework consists of the phase specific, ensemble averaged transport equations of the Euler-Euler two-fluid model, a set of closure models for interfacial momentum transfer terms for regimes of dispersed flows, an interfacial drag formulation to represent large, resolved interfaces in the two-fluid model and a morphology adaptive weighting procedure to apply closure models consistently to individual continuous and disperse phases.

Basic equations of the two-fluid model
The ensemble averaging procedure results in a set of mass and momentum conservation equations for each individual phase 44 Equation (1) shows the phasic mass conservation equation, which will later be utilized to realize the transport of phase volume fractions r . Time and spatial partial derivatives are denoted as t = / t and i = / x i , respectively, and u ,i as velocity field of phase . In the phase momentum conservation Equation (2), p represents the pressure, which is shared between all phases . The phase dynamic viscosity is , the phase strain rate tensor is the gravitational acceleration is g. Surface tension for resolved interfaces is modeled with the approach of Brackbill et al, 45 where denotes the surface tension coefficient for the pair of both continuous phases and with according interface curvature and interface normal vector n ,i . Volume forces are denoted by f ,i , which includes all interfacial momentum exchange closure terms.

2.2
Interfacial momentum exchange for dispersed and resolved interfacial structures in the two-fluid model In the disperse two-fluid model, all information about interfacial structures and interactions between different phases is lost due to the applied averaging procedure. Hence, interfacial momentum exchange mechanisms between both phases of each continuous-disperse phase pair have to be modeled. The following interfacial forces are usually considered in modeling of dispersed flows: 46 drag, lift, virtual mass, turbulent dispersion, and wall lubrication. In bubbly flows, drag and virtual mass forces are of great importance 9 and their consideration is required to achieve a solution procedure that is numerically stable. Therefore, only interfacial drag and virtual mass forces are accounted for, together with related closure models, which are selected according to Liao et al (2019). 46 Nevertheless, the algorithm is capable of handling all other interfacial momentum exchange forces. For drag forces, the closure model of Ishii and Zuber 47 is applied without any correction terms considering swarm effects. The closure term by Crowe et al 48 is applied to model the virtual mass force between a disperse phase d and a continuous phase c The substantial derivative of a quantity is denoted as D t, , considering phase velocity u . The model constant C VM cd is set to 0.5. For a pair of two continuous phases cc, a virtual mass force is not defined and, therefore, C VM cc is set to zero. Whenever large-scale interfaces occur, such flow structures have to be resolved within the computational grid. In order to keep the thickness of an interface finite, local velocity components normal to the interface have to be identical for both continuous phases, which share the interface. This requirement is often referred to as the interface normal condition. In a one-fluid model, which is also known as homogeneous model, both continuous phases share a single-velocity field, or in other words, all velocity components are identical, which is a more restrictive constraint than only ensuring the equality of interface normal components. Such a no-slip interface condition is reasonable, if the resolution of the computational grid is fine enough to resolve interfacial structures. 9 In contrast to a homogeneous model, the two-fluid model is able to allow for a tangential slip velocity between both phases, which is an appropriate approach to model interfaces in an underresolved manner. 39 For the sake of simplicity in the present work, the model framework is setup, such that all relative velocity components between multiple phases are identical in the region of a shared interface. In the future, the interfacial drag modeling in resolved interface regions will be made capable to properly handle a free-slip condition. In this way, the capability of the present model is demonstrated to reproduce the behavior of a homogeneous model, which has been shown theoretically by Yan and Che. 16 This can be achieved by application of an interfacial drag model, which is designed to ensure that both phases have identical velocities in the interface region.
Štrubelj and Tiselj 42 proposed a drag closure, which is formulated for a pair of two phases and The mixture density is formulated as The drag coefficient K D is formulated in terms of relaxation time r , which has to be chosen much smaller than the physical time step. 42 The model is applied to the present framework with a small value for the relaxation time r = 10 −8 Δ t in order couple continuous phases in an interfacial region, even if dispersed phases are present at the same location.

Morphology adaptive weighting of closure models
The present model framework utilizes an arbitrary number of continuous and disperse phases and each of those is dedicated to a physical phase, for example, air, oil, or water, and to an individual flow morphology, for example, continuum, bubbles, or droplets. Hence, special care has to be taken of interfacial momentum transfer modeling, whenever three or more phases occur in the same location. To account for such situations, a weighting factor F is inserted into the formulation of all individual momentum exchange terms between two phases and wheref denotes the original closure formulation. It is worth noting that this type of regime adaption is purely related to the fact that the presented hybrid model distinguishes between continuous and disperse phases. The weighting factor between any pair of continuous phases F cc is set to unity. In that way, the strong velocity coupling between velocity fields of different continuous phases via the drag model formulation of Štrubelj and Tiselj 42 is maintained in order to fulfill the no-slip interface condition. Between two disperse phases (dd), no momentum exchange terms are defined and, therefore, F dd is set to zero. For all momentum exchange terms of a pair of disperse and continuous phases (cd), the following weighting factor is applied, where C denotes the set of all continuous phases. As F cd equals unity, wherever exactly one continuous phase is present, the original closure model is recovered, even if there are two or more disperse phases, resulting in the intended behavior of the original closure model. In an interfacial region, where multiple continuous phases are existent, the prefactors F cd for all pairs of one and the same disperse phase d and any of the continuous phases c sum up unity. In that way, all momentum exchange closures, for example, the drag force, are accounted for in the right amount and the disperse phase is able to pass the interface, which will be demonstrated in Section 4.

Equations
Equation (2) results in a set of partial differential equations that is solved numerically on a collocated grid. The coupled system of equations is solved in a segregated manner by linearization and subsequent application of the Jacobi method. Thereby, the resulting equational system is solved iteratively. In the following, the notation of Cubero et al 41 is adopted, so variables at cell centers and faces are denoted with [] P and [] f , respectively, and the linear interpolation operation from cell centers to cell faces is depicted with ⌊⌋ f . The linear system of equations is expressed in matrix form The implicit diagonal matrix coefficients are referred to as a, explicit off-diagonal contributions, momentum exchange closures and volume force contributions as H, and the contribution due to pressure force is denoted as G∝ − ∇p. At this point, implicit and explicit contributions of the transient term are included in coefficients a and H, respectively. In order to suppress decoupling of pressure and velocity fields (checker boarding) due to the stencil structure of the discretization scheme for pressure gradient forces, Rhie and Chow 49 proposed a face-based momentum interpolation algorithm. The numerical flux [ ] f is calculated by formulation of Equation (9) on cell faces f and division by the implicit coefficient a: Wherever r takes a value of zero, residual contributions are added to both coefficients, a and H, in order to obtain a numerical solution from Equation (10).

Compact momentum interpolation (CMI)
For the multifluid model Cubero and Fueyo 50 and Cubero et al 41 showed that the above face-based formulation does not satisfy the four requirements of consistency: A solution procedure is considered consistent, if it fulfills the following requirements: 41 1. If ever achieved, a steady-state solution has to be independent from the size of the time step, 2. Wherever only a single phase is present, the single-phase CMI formulation needs to be recovered inherently, 3. All interfacial drag forces in the momentum equations need to cancel out, 4. The pressure forces acting on individual phases have to proportional to the forces resulting from the shared pressure.
Hence, they proposed the so-called compact momentum interpolation with the pseudo-flux where n f is the normal vector of face f. The basic idea of CMI is to split the implicit and explicit coefficients a and H for the individual phases into separate parts for convective-diffusive (without superscript), temporal (T), and interfacial drag (D) contributions, which are separately interpolated to the faces. Temporal coefficients are defined as a T = r V∕Δ t with volume V of a grid cell and time step width Δ t . Time integration is formulated with first-order Euler backward scheme, where n denotes quantities from the recent time step. The extension to second order time integration scheme is described by Shen et al. 51 As the system of equations is solved iteratively in a segregated manner, the numerical flux [ ] f of another phase is taken from the recent Jacobi iteration. In contrast to Cubero et al, 41 interfacial drag is formulated for an arbitrary number of phases by using a sum formulation in the present work. The drag coefficient is defined as a D = F K D , including weighting factor F , according to Section 2.

Consistent formulation of virtual mass
The formulation of the virtual mass force needs to be included into the solution procedure in a way that the requirements for consistency from Section 1 are met. Therefore, the virtual mass coefficient is The model of Crowe et al, 48 which is shown in Equation (4), contains material derivatives and, therefore, convective and temporal contributions. The convective contribution is and is added to both implicit and explicit convective-diffusive parts of Equation (11), a and H . The same treatment is applicable to other interfacial closure models, namely turbulent dispersion, lift, and wall lubrication forces, which are excluded in the present work for sake of simplicity. Analogous to temporal terms (T) temporal virtual mass contributions (VM,T) are treated separately with coefficient a VM,T = a VM ∕Δ t . Both convective and temporal contributions of the virtual mass force are added to Equation (11): This formulation combines the formulation of all interfacial models used in the present work with CMI, which allows for a consistent pressure-velocity coupling.

Approximate resolution of phase coupling
As stated in Section 2, drag coefficient a D is chosen to be very large in order to ensure stability of resolved interfaces, even in the presence of more than two phases at the same location. In the limit of an infinitely large a D , this coefficient dominates Equation (15) and both phase velocities u and u become identical. This results in underrelaxation of phase velocities and phase fluxes and the system of equations becomes stiff. Multiple approaches have been proposed in the literature to resolve such situations of stiff phase coupling. For an arbitrary number of phases the coupling can be fully canceled by solving the phase momentum Equation (2) for all phases in a coupled manner, as proposed by Karema and Lo. 52 A solution procedure utilizing full-block coupling of all phase velocities as well as shared pressure in combination with CMI was presented by Ferreira et al, 53 who proved the phase coupling to be resolved for a number of up to four phases for dispersed flow regimes. A different approach is the partial elimination algorithm (PEA) of Spalding, 43 which is algebraically formulated for the quantity of exactly two phases. Compared with coupled approaches, PEA can be included into any segregated solution procedure, which is available in OpenFOAM-dev. 54 The basic idea is to eliminate terms in the phasic momentum Equation (2), which contain the phase velocity or the phase flux of the other phase by substitution and rearranging of those equations. A full cancelation of phase coupling can only be achieved for exactly two phases. Nevertheless, an arbitrary number of phases occurs at the same location and, therefore, phase coupling becomes resolved in an approximate manner. The steps of the derivation 43 are applied to Equation (15) in the following way: 1. sum up Equation (15)  For the whole procedure, the sum formulation needs to be maintained considering all phases , , ≠ , ≠ { , }. This results in CMI equations with approximately resolved phase coupling with implicit drag coefficient and explicit drag contribution .
It is obvious that there are still contributions and from third phases ( , ≠ { , }) in the explicit drag contribution. This means that phase coupling is not fully resolved, whenever there are more than two phases present. The original PEA algorithm is recovered in the presence of exactly two phases. It was found that number of phases that equivalently dominate the momentum equation in terms of mass has to be limited, which is clearly the case, for example, for a single liquid phase and multiple gas phases. It is expected that typical gas liquid flow applications are within the limit of this multiphase formulation of the PEA algorithm.
In the limit of infinitely high drag coefficients for resolved interfaces (see Equation (5)), the proposed system of equations is equivalent the homogeneous model, 16 such that the behavior of an algebraic VOF is recovered with the presented method. On the other end, if for dispersed flows closure models of Ishii and Zuber 47 and Crowe et al 48 are used for drag and virtual mass, respectively, the presented method acts similar to the well-known Euler-Euler model.

Solution procedure
The presented solution procedure utilizes a projection method, 55 which consists of the following steps: 1. calculation of predicted phase velocities [u * ] P and fluxes [ * ] f from Equation (16) by omitting pressure contribution G, 2. solution of pressure equation, and 3. obtain correct velocities and fluxes from Equation (16) with pressure contribution.
In the case of incompressible flow, the pressure equation is formulated as with pressure diffusivity The entire solution procedure is summarized in Figure 1.

Discretization
The presented numerical framework is implemented into OpenFOAM, 40,54 which is an open-source C++ library based on a finite-volume method formulated for unstructured grids. The solution procedure utilizes a discretization that is second-order accurate in space and provides an implicit time integration scheme of first order. The transport of phase fractions is realized with explicit time integration of phase fraction transport Equation (1). Boundedness is maintained with the multidimensional limiter for explicit solution algorithm (MULES) 54 together with the flux limiting scheme as proposed by Van Leer 56 for the discretization of convective transport terms of the phase fractions. Due to the explicit nature of the time integration scheme for phase fractions in OpenFOAM, Courant numbers below unity are required. For the accurate prediction of transient behavior, Jasak 57 reported that an even smaller maximum Courant number is necessary. The interface compression mechanism according to Weller 7 is included in the phase fraction transport algorithm in order to keep resolved interfaces from diffusing with a model constant of c = 1.

Variables
In the cases considering resolved rising bubbles, the initial bubble diameter is referred to as D b , which is used to define any dimensionless length asL = L∕D b , where L can be any scalar or vector quantity with a length dimension, for example, a position vectorx or diameter of disperse interfacial structuresd. Furthermore, gravitational velocity and time scales are defined by U g = √ gD b and t g = √ D b ∕g, respectively, with gravitational acceleration g. Dimensionless time and velocity are defined ast = t∕t g andŨ = U∕U g , respectively.

Two-dimensional rising bubble
The first test case setup, which features a two-dimensional gas bubble, rising in a stagnant liquid, has been proposed by Hysing et al, 58 Resulting dimensionless numbers for both cases under investigation are shown in Table 1.
Results are compared in terms of evolution of bubble center of gravity, rising velocity, and circularity over time as well as of interface location at a timet = 4.2. The center of gravity of the bubble, x b , is defined as: 59 with domain Ω ⊂ R 2 . The vertical component of x b directed opposite to gravitational force is referred to as y b . The rising velocity U b of the gas bubble is defined accordingly: 59 Its vertical component is denoted as U b . In two-dimensional space, the circularity is defined as 58 with circle-equivalent diameter d cir , actual bubble area A b , and actual bubble perimeter P b . In order to assess the influence of spatial resolution, a mesh study is carried out with the domain being spatially discretized with 20 × 40, 40 × 80, 80 × 160, and 160 × 320 grid points, which results in Δx = 1∕10, Δx = 1∕20, Δx = 1∕40, and Δx = 1∕80, respectively. Time step sizes are set to Δt = 1.4 × 10 −2 , Δt = 7 × 10 −3 , Δt = 3.5 × 10 −3 , and Δt = 1.75 × 10 −3 , respectively, such that the Courant number stays Co < 0.12 for each individual mesh level. The solver is set up, such that the number of inner loops is fixed to three, whereas outer loops are done until the initial pressure residual at the beginning of an outer loops is lower than 1 × 10 −4 . Subsequently, one final outer iteration is done before proceeding with the next time step. Vertical bubble positionỹ b , rising velocityŨ b , and bubble circularity c b over time as well as the interface location (r G = 0.5) att = 4.2 are presented in Figure 2.
For the vertical bubble position, no major differences can be observed, while the rising velocity noticeably differs between mesh levels Δx = {1∕10, 1∕20} and Δx = 1∕80 at timet ≈ 1. Fort = 4.2, the rising velocity on the coarsest mesh is significantly larger than on the finer ones. Almost no differences in rising velocity can be observed between meshes Δx = {1∕20, 1∕40, 1∕80} att = 4.2. On coarse meshes, circularity values of c b > 1 can be observed in the beginning of the simulation. Mathematically, the circularity of an object cannot be larger than unity, as the circle, which has c b = 1 by definition, is the shape with the lowest possible ratio of perimeter to surface area. The simulation results can be explained with the error of the procedure for estimating the perimeter of the bubble, which is quite large on coarse grids. Nevertheless, the qualitative behavior of the bubble is recovered with all meshes used. While the resulting differences between the coarsest and the finest meshes are large, minor deviations between the results of the meshes Δx = {1∕40, 1∕80} can be observed fort > 3. The interface position att = 4.2 converges nicely with refinement of the computational grid and meshes Δx = {1∕40, 1∕80} deliver nearly identical results. Hence, mesh level Δx = 1∕40 is chosen for further investigations of this case. Although the presented results suggest that total bubble volume changes between meshes of different spatial resolution, gas volume is conserved. This is evident by integrating r G over the domain, giving a relative volume error of 1.02 × 10 −8 for 0 <t < 4.2 for the coarsest mesh. The reason, why bubbles appear smaller on a coarse mesh with the interface definition of r G = 0.5 is that the interface smears over several grid cells, which is typical for algebraic VOF methods. Figure 3 shows the results of case 1, namely vertical bubble rising velocityŨ b and circularity c b over time, as well as the interface position (r G = 0.5) at timet = 4.2 and the relative pressure residuals for each inner and outer loop for the latest time step att = 4.2. The present solver's performance is compared with the results, which Hysing et al 58 Table 2. For the maximum rising velocityŨ b,max , the present solver predicts a value, which is in between both reference values, being slightly closer to the data of Klostermann et al 60 than to the ones of Hysing et al. 58 The same accounts for the period of time, which is needed until maximum rising velocity is reached. Subsequently, the bubble slows down, which is most pronounced in the results of the present solver. Finally, att = 4.2, the rising velocity as predicted by the presented model is slightly higher than the one of Klostermann et al 60 but at the same time lower than the value predicted by Hysing et al. 58 For the circularity, the present solver predicts a minimum value of c b, min , which is below the values of the reference data from both sources. Att = 4.2, the circularity value obtained with the presented model is between the reference values. The bubble shape at this point of time slightly differs from the results of Hysing et al 58 but is very close to the ones of Klostermann et al. 60 Att = 4.2, the present solver predicts a value for the vertical bubble position that again lies between the results of both references. Regarding the relative pressure residuals over inner and outer loops of the recent time step att = 4.2, a good convergence behavior is observed inside each individual outer loop as well as over the progress of all outer loops. The four outer loops express in the four peaks and the final relative residual for the whole time step is quite low with a value below 10 −8 .  For case 2, the geometrical setup is identical to the one in the first test case. Gravitational acceleration and material properties are chosen, such that they result in the dimensionless numbers of this case, as defined by Hysing et al 58 and listed in Table 1. As for case 1, a convergence study regarding spatial resolution is carried out as well. Additionally to the meshes investigated for case 1, another mesh refinement level with a grid spacing of Δx = 1∕160 is used together with a time step size oft = 8.75 × 10 −4 . The time step sizes for all other mesh refinements levels are set according to the refinement study in case 1. The results of this mesh study regarding the bubble's position, vertical rising velocity and circularity over time as well as the interface position (r G = 0.5) att = 4.2 are shown in Figure 4.
It turns out that a higher bubble rising velocity is predicted on finer meshes, which results in the highest bubble position for the finest mesh refinement level. In the beginning of the simulation fort < 2, except for the coarsest mesh, the results show only minor differences in terms of rising velocity. Subsequently, a second local velocity maximum appears, which is most pronounced on fine meshes. This results in a persistent velocity difference between the results of the different mesh levels. With the coarsest mesh, the second local velocity maximum is not captured at all. In terms of circularity value c b , a convergent behavior can be observed with increasing spatial resolution with the largest differences just between the coarsest mesh Δx = 1∕10 and the next finer mesh level Δx = 1∕20. Comparing the interface position (r G = 0.5) att = 4.2, large differences are present between the results of different spatial resolution levels. On the coarsest mesh, the skirt in the lower part of the bubble cannot be captured. Mesh level Δx = 1∕20 results in a bubble skirt, which is of nearly uniform thickness. With increasing mesh resolution, the skirt becomes more and more slim in the center, while the lower part nearly retains its thickness. Finally, with the finest mesh Δx = 1∕160, two gas portions are already detached from the skirt att = 4.2. For further investigations, mesh resolution Δx = 1∕160 is chosen for this case. As for case 1, the results obtained with the present solver are again compared with the data of Hysing et al, 58 who used the TP2D solver and a mesh resolution of Δx = 1∕320 (originally denoted as 1/h = 640), as well as with the data obtained by Klostermann et al 60 with interFoam and a mesh resolution of Δx = 1∕160 (1/h = 320). In Figure 5, vertical bubble rising velocity and circularity over time as well as interface position (r G = 0.5) and relative pressure residuals for timet = 4.2 are shown.
Numerical values of local extrema, as well as vertical bubble positionỹ b , rising velocityŨ b , and circularity c b at timẽ t = 4.2 are reported in Table 3. For the first local maximum of the rising velocityŨ I b,max , the value obtained with the present solver again lies between the results of Hysing et al 58 60 and the present solver predict that the skirt region is still more pronounced and two slightly larger gas structures just separate from there. Those results show only minor differences, which can also be observed in the evolution of the circularity value with a final deviation of about 1.5%. The initial relative pressure residual over inner and outer loops shows a convergent behavior, even though the final level of convergence is, with a value of about 10 −6 , not as deep as in case 1. At the final time step, five outer loops are necessary to reach the criterion of 10 −4 , which expresses in the five peaks in the residuals.
In summary, the present solver is able to recover the dynamics considering a single two-dimensional rising gas bubble. Moderate deviations are observed in comparison with the results obtained with a level-set method by Hysing et al, 58 while very small differences occur in comparison with results of the VOF method without interface reconstruction by Klostermann et al. 60

Three-dimensional rising bubble
The three-dimensional case setup was proposed by Balcázar et al, 61 who investigated the performance of a level-set method and compared results with experimental data. The gas bubble is initialized as sphere with diameter D b , with the center position 2D b above the bottom. The computational domain is described as a cylindrical vessel with diameter 8D b and height 12D b . As in the two-dimensional cases, a free-slip condition is applied to the side wall of the cylindrical domain. Contrary to the setup of Balcázar et al, 61 a no-slip boundary condition is imposed at the bottom and pressure is fixed at the top boundary instead of applying periodic boundary conditions. A free gas-liquid surface is established 1D b below the upper boundary. For illustration, the computational domain with dimensions, the computational grid as well as the initial interface position are shown in Figure 6. Two cases have been selected from the ones investigated by Balcázar et al, 61 whose corresponding dimensionless numbers are shown in Table 4. Additionally to the dimensionless numbers, which are introduced in the previous section, flow regimes are categorized by the bubble Reynolds number with rising velocity U Bal b , which is the numerical result for the bubble rising velocity of Balcázar et al. 61 and (21) with Ω ⊂ R 3 . Similarly to circularity c b , the bubble sphericity s b is defined as the ratio of the surface area of a sphere with equivalent volume V b to the actual bubble surface area S b : The computational domain is discretized, such that in the center of the domain, where the bubble is located, computational cells are orthogonal and of constant grid sizing. Outside the core region, cells are getting larger the further their distance to the core is in order to save computational cost.
Influence of spatial resolution is assessed via comparison of results obtained with computational grids of spacings Δx = 1∕8, Δx = 1∕16, and Δx = 1∕32. Time step sizes are set to Δt = 1.25 × 10 −1 , Δt = 6.25 × 10 −2 , and Δt = 3.13 × 10 −2 , respectively, and the Courant number is Co < 0.12 for all simulations. As in the two-dimensional cases, the number of inner loops is set to a fixed number of three. Outer loops are run until the initial pressure residual is lower than a value of 1 × 10 −3 , before the final outer iteration is carried out. Figure 7 shows evolutions of vertical position, rising velocity, and sphericity of the bubble over time as well as bubble interface position (r G = 0.5) att = 12.53 for the different mesh levels in case 1.
On the coarsest mesh with Δx = 1∕8, the predicted bubble rising velocity decreases over time after it reaches its maximum value att ≈ 2, while the value of rising velocity is maintained on the finer meshes. Therefore, for the coarse mesh, vertical bubble position increases slower over time compared with finer spatial resolutions, which leads to a significant vertical distance att = 12.53 between the bubble on the coarse and on both finer meshes. Between Δx = 1∕16 and Δx = 1∕32, slight differences in bubble rising velocity can be observed for 1 <t < 6 with a lower value with the finer mesh, before they approach each other again. In terms of bubble sphericity s b , unphysical values larger than unity can be observed for both the coarser meshes, which has the same reason as explained for bubble circularity c b in the two-dimensional cases. Furthermore, the graph of the sphericity for Δx = 1∕8 appears to be very jagged for the same reason that the precision of the interface reconstruction for postprocessing is very poor for coarse spatial resolutions. For the finer meshes, the graphs are smooth and moderately differ from each other fort < 10, before they take very similar values. In terms of final interface position att = 12.53, the bubble shape obtained with the coarsest mesh is too wide and flat compared with the results on the other meshes. A direct comparison of bubble shapes from Δx = 1∕16 and Δx = 1∕32 shows that the shape is very similar in the upper region, while the dent between the skirt structures is less pronounced on the finest mesh. The overall differences between the results obtained with both finer meshes are moderate and, hence, mesh Δx = 1∕32 is chosen for further investigations.
The performance of the present solver is compared with the experimental results of Bhaga and Weber 62 as well as to the numerical results of Balcázar et al, 61 who used a mesh of spatial resolution of Δx = 1∕30. Additionally, the three-dimensional bubble is simulated with a solver based on the homogeneous model (interFoam, 54 ) with mesh spacing Δx = 1∕32. Bubble rising velocityŨ b and sphericity s b over time as well as interface location (r G = 0.5) at timet = 12.53 for case 1 are shown in Figure 8. It is worth noting that the data for interface position from both Bhaga and Weber 62 and Balcázar et al 61 Table 5. Fort < 4, the rising velocity predicted with the present solver is noticeably smaller than the one obtained with interFoam. In the second half of the simulation, both values approach each other. Throughout the whole simulation, the values for rising velocity reported by Balcázar et al 61 62 Concerning the interface position at final time, the numerically obtained bubble shape reported in the literature shows blunt tips for both skirt regions, which is in contrast to the results obtained with interFoam, which show a steeper rise of the interface in the inward direction of both skirt regions. With the present solver a bubble shape is predicted, which has even more pronounced ligaments of the skirt with the tip region being more sharp than the one obtained with interFoam. That is completely in line with the numerical values for bubble sphericity at final time.
Apart from the skirt regions, nearly no differences in bubble shape can be observed between the results of both solvers and of Balcázar et al. 61 The comparison with the experimentally obtained bubble shape shows a good agreement. For assessment of influence of spatial resolution in case 2, the same computational grids as for the investigation in case 1 are used together with the according time step sizes. Additionally to that, another level of mesh refinement is added with Δx = 1∕64 with a time step size of Δt = 1.56 × 10 −3 . With this setup, Courant numbers of Co < 0.18 are achieved on all mesh levels. Results in terms of bubble vertical positionỹ b , rising velocityŨ b , and sphericity s b over time as well as interface location at timet = 12.53 for case 2 are shown in Figure 9. As before, the definition of the interface location for meshes Δx = 1∕16 and finer is r G = 0.5 but for Δx = 1∕8 the interface smears out strongly and the distribution of the phase fraction of the gas is so dilute that the criterion mentioned before is met nowhere. Therefore, an interface definition of r G = 0.1 is used for that coarsest mesh level.
As it turns out, the coarsest mesh cannot capture the flow dynamics at all and the bubble breaks up into multiple small, quite dilute gas structures. Connected to that are very low values of rise velocity and vertical position compared with results obtained on finer meshes. Accordingly, the evolution of sphericity over time becomes very jagged and rapidly rising in the second half of the simulation. In terms of vertical position, rising velocity sphericity, and final interface position, all other meshes show a convergent behavior with increasing mesh resolution and only minor deviations are observed between the values obtained with meshes Δx = 1∕32 and Δx = 1∕64. Hence, mesh Δx = 1∕32 is chosen for further investigations.
For case 2, Figure 10 shows results in terms of bubble rising velocityŨ b and sphericity s b over time as well as the location of the gas-liquid interface (r G = 0.5) at timet = 12.53, which are obtained with the present solver, the homogeneous model solver, experimentally obtained by Bhaga and Table 4). Furthermore, literature reference data for interface position are scaled and shifted, analogously to case 1.
In contrast to case 1, the rising velocity shows two pronounced and one flat local maximum value,Ũ I b,max ,Ũ II b,max , andŨ III b,max , in the early phase of the simulation (t < 6). Numerical values of those velocity extrema and their time of occurrence as well as values of vertical bubble position, rise velocity, and sphericity at timet = 12.53 are presented in Table 6. Considering rise velocity and sphericity, both the present model and the homogeneous model perform very similar with relative deviations below 1%. While the homogeneous model predictsŨ b to be slightly higher in the beginning of the simulation than the presented model does, the situation is vice versa fort > 2. Both solvers predict that the bubble slowly decelerates in the second half of the simulation, making it more than 3% slower than reported by Balcázar et al. 61 The deviation of the present result from the experimental value of Bhaga and Weber 62 in terms of final rising velocity is approximately −8 %. In terms of bubble shape, the results of Balcázar et al 61 show that the bubble is nearly flat in the lower region between the ligaments of the skirt region. In contrast to that, the other solvers predict a bulge of the interface in the region of the symmetry axis. Moreover, the ligaments are predicted to be longer and more slender compared with the results of Balcázar et al. 61 The bubble shape as obtained with the presented solver has ligaments, which are slightly more blunt at the tips compared to the homogeneous model solver. In contrast to the numerical results, the experimentally obtained bubble shape shows that ligaments at the lower end of the gas structure are pinched toward the center line, while deviations regarding interface position at the upper part of the bubble are rather small.
Overall, the presented model frame work performs very similar to the solver, which is based on the homogeneous model, while deviations from the results of Balcázar et al 61 are moderate.

Two-dimensional bubble rising in liquid with micro gas bubbles
In order to demonstrate the present solver's capability to describe the interplay between large and finely dispersed interfacial structures, the two-dimensional case setup 1 from Section 4.2 is modified such that a large gas bubble rises trough liquid with a layer of micro bubbles, which consist of the same gas. In contrast to the original setup in Section 4.2, the computational domain is extended in vertical direction, resulting in dimensions of 2D b × 20D b . The spatial resolution is kept constant, such that 80 × 800 grid points are used for discretization. Micro bubbles are modeled as disperse Eulerian phase (d) with fixed bubble diameter for interfacial models of d d = 2 × 10 −3 D b . The disperse phase is initialized with a void fraction of r d = 0.1 in the region of 3 <ỹ b < 4. Disperse drag and virtual mass formulations as described in Section 2 are applied to model momentum exchange between disperse and continuous phases. All walls are modeled with free-slip conditions for the disperse phase. As in case 1 of Section 4.2, time step size is fixed to Δt = 3.5 × 10 −3 , which results in Co < 0.133. The large gas bubble is rising and deforming just as described in Section 4.2 and after it passes the region, where the micro bubbles are present, a portion of dispersed phase gets entrapped into the wake region of the large gas structure. The velocity fields of the large bubble and the distribution of entrapped micro bubbles at timet = 63 are shown in Figure 11. At this time, the bubble's vertical position isỹ b |t =63 = 17.27. It turns out that the micro bubbles distribute in close relation to the streamlines of the liquid-phase velocity outside the large gas structure and, therefore, the majority of them stays in the recirculating liquid flow for a long time. Inside the resolved gas bubble, the continuous gas phase also develops quasisteady vortex structures and, as dispersed and continuous gas phases have identical material properties, mixes with the disperse gas phase, which concentrates at the lower part of the gas bubble. From there, small disperse gas portions get sucked upward to the bubble cap by the internal coherent vortex structures. In the frame of the Euler-Euler method is not a drawback, as the phase fraction represents a probability due to the implicated averaging procedure. This demonstrates that the present solver is able to handle both subgrid and supergrid interfacial structures at the same time.
Further modeling is required to handle interactions of small bubbles with a resolved interface properly.

Two-dimensional stagnant stratification of water and oil with air bubbles
In order to further demonstrate the solvers capability to handle interactions between disperse phases and interfaces between continuous phases at the same time, dispersed air bubbles of diameter d air are introduced into a stagnant stratification of oil and water. The two-dimensional computational domain of dimensions 0.15 m × 0.5 m is discretized with 30 × 100 grid points, and at time t = 0 s, the lower half of the domain is filled with continuous water phase, whereas the upper part is filled with continuous oil phase. The top boundary is setup to provide a fixed-pressure condition. For lateral and bottom boundaries no-slip conditions are imposed for the continuous water and oil phases. For the disperse air phase, these boundaries act as free-slip walls. At the bottom region between 0.05 m < x < 0.1 m, disperse air is introduced into the domain with a volume fraction of r d = 0.26 and a wall-normal velocity of U d = 0.197 m s −1 , where both r d and U d are linearly ramped up from zero values for time 0 s < t < 0.3 s, kept constant until t = 0.6 second and again ramped down linearly to zero values for 0.6 s < t < 0.7 s. Again, drag and virtual mass closure formulations from Section 2 are applied to model interactions between dispersed air phase and both continuous water and oil phases. Time step size is fixed to a value of Δ t = 0.03 s, resulting in a Courant number Co < 0.296. It is expected that the dispersed air bubbles rise due to gravity g = 9.81 m s −2 and pass the water-oil interface smoothly. The bubble rise velocity should slightly change after passing the interface. Physical properties for the present case are reported in Table 7. The initial residuals for the solution of pressure Equation (17) over one single time step at t = 1.7 s are shown in Figure 12. A total number of i = 10 outer and j = 2 inner loops is used. It turns out that initial residuals stay nearly constant over the outer loops for i ≥ 5 with relative residuals smaller than 2 × 10 −7 for each second iteration (j = 2). Figure 13 shows the temporal evolution of spatial distribution of dispersed air volume fraction r d . After starting to rise in the water from the bottom, the swarm of air bubbles forms a mushroom-like shaped cloud, which then hits the water-oil interface at t ≈ 1.6 s, penetrates it and continues rising to the top boundary, where it exits the computational domain. The water-oil interface becomes temporary deformed by the bubble swarm and keeps sloshing for some time without being completely disrupted. After the air has left the domain through the outlet, the resolved interface almost perfectly recaptures its initial vertical position of y| t = 0 s = 0.25 m, which means that the total amount of water is conserved in the domain over the course of the simulation within an accuracy of 0.398% after t = 10 s. From the contours of the vertical component of the air velocity U a, y ( Figure 13H-N) it becomes clear that the dispersed air phase passes the sharp water-oil interface smoothly without any oscillatory behavior in the phase velocity field. Furthermore, in the air velocity field a region can be identified ahead of the bubble swarm, where air bubbles have a nearly constant rising velocity, which results from the force balance between buoyancy and drag between small bubbles and the continuous liquid around. Figure 14 shows vertical profiles of phase volume fraction r as well of vertical phase velocity components U ,y .
It is visible that the interface between water and oil phases maintains its relatively small thickness, even when dispersed air phase is present at the same location. Phase velocities U ,y of water and oil are close to be perfectly identical everywhere in the domain, which is the desired behavior and results of the choice of the interfacial drag formulation by Štrubelj and Tiselj 42 together with the small value for relaxation time r . The vertical air velocity is also nearly identical to the ones of continuous phases, wherever air volume fraction is close to zero. In regions, where air is present, its vertical velocity component is approximately 0.15 m s −1 higher than vertical velocity components of the continuous phases, which results from the buoyancy force acting on the gas bubbles and accelerating them in vertical direction. Again, all velocity profiles turn out to be very smooth, especially in the interface region, where three phases are present in this case.

SUMMARY AND CONCLUSIONS
The compact momentum interpolation method was applied to the multifield two-fluid model. Special care was taken to properly define phase coupling in regions of large, resolved interfacial structures by applying the interfacial drag formulation of Štrubelj and Tiselj. 42 The resulting stiffness of the system of multiple momentum conservation equations for the individual phases was resolved in an approximate fashion by extending the idea of partial elimination 43 to an arbitrary number of phases by maintaining a sum formulation, wherever more than two phases are present. This method worked out to resolve phase coupling in the investigated test cases with a maximum of three phases while preventing the solution procedure from heavy underrelaxation effects. Furthermore, this method was extended to consistently treat the virtual mass force, which is a momentum exchange model for description of dispersed flows in the used multifield two-fluid model. For the first time it has been shown that dispersed as well as resolved interfacial structures can be simultaneously solved with the same set of equations considering interactions between flow and interfacial structures of different length and time scales in a computationally efficient way. The next steps are to include further interfacial closure models, namely lift, wall lubrication, and turbulent dispersion forces, so that the typical set of momentum exchange closure models for dispersed flows can be applied. Furthermore, mass transfer between phases describing different morphologies of the same physical phase, such as dispersed air bubbles and continuous air, will be added, for example, in order to model the bursting of small air bubbles at a large water-air interface. Further extensions are the consideration of drag modeling regarding resolved interfaces in an underresolved manner by allowing for interfacial slip, subgrid scale modeling in the fashion of large-eddy simulations as well as modeling of further multiphase morphologies such as froth.