Momentum exchange modeling for coarsely resolved interfaces in a multifield two‐fluid model

Morphology‐adaptive multiphase models are becoming more established for the numerical description of complex gas‐liquid flows adapting dynamically to the local flow morphology. In the present study, two different numerical methods originally designed for distinct flow morphologies are combined, namely the volume‐of‐fluid and the Euler–Euler method. Both edge cases have been proven to be capable of delivering reliable predictions in the respective use cases. The long‐term goal is to improve the prediction of gas‐liquid flows, regardless of the flow regime in a specific application. To capture the system dynamics with a given grid resolution, the flow fields need to be predicted as precise as possible, while the shape of structures such as gas bubbles need to be recovered adequately in topology and shape. The goal is to obtain reliable predictions on intermediate mesh resolutions rather than relying on fine meshes requiring more computational resources. Therefore, a procedure is proposed to locally measure the degree of resolution. With this information, the hydrodynamics in the interface region can be controlled by means of a dedicated interfacial drag formulation in order to improve simulation results across several levels of spatial resolution. A modified formulation of buoyancy is proposed to prevent unphysical oscillations of vertical velocity near a horizontal interface. The functionality is demonstrated in a three‐dimensional case of a gas bubble rising in stagnant liquid and in a co‐current stratified air‐water channel flow in two‐dimensional space. The choice of these different applications demonstrates the general applicability of the proposed model framework.

may appear at the same time in different locations at given operating conditions. Typically, features of the interface, that is, the boundary between two immiscible phases, are divided into large-scale interfaces and dispersed interfacial structures, such as small bubbles, droplets or particles of varying size, and shape. Both aforementioned morphologies can occur simultaneously and may interact with each other or even evolve, such that they become transformed to each other. Industrial applications featuring such phenomena are, for instance, centrifugal pumps, 1,2 swirling flow separators, 3 valve trays, 4,5 or pressurized thermodynamic shock scenarios. 6 The numerical methods for fluid dynamic simulations are typically limited in terms of length and time scales that can be resolved. Scale-resolving approaches, such as the volume-of-fluid (VOF) method, 7 are designed to predict large interfaces in relation to the computational grid. With statistical methods, such as the Euler-Euler model, dispersed structures are described as sub-grid scale interfaces. In order to cover wider ranges of scales morphology-adaptive multiphase models (hybrid models) are being developed as a combination of both approaches. A variety of such modeling techniques for multiphase flows exists, that allow for morphology transitions, also referred to as inter-scale transfers. 8 Methods can be distinguished by the underlying averaging method, for example, Euler-Euler versus Euler-Lagrange. Moreover, the number of equations and other technical details can be used for further specification, for example, domain-wise blending of level-set and Euler-Euler models or drift-flux models, which are based on VOF. A comprehensive overview of such methods is beyond the scope at this point, but can be found in literature. 9,10 The present work focuses on approaches that are based on the two-fluid models, that is, transport equations for mass and momentum are stated for each phase. There are two major modeling strategies when modeling dispersed and resolved interfaces based on the two-fluid model: 1. the two-field approach, in which all physical phases locally adapt their representation and behavior by blending between different closure laws depending on the local conditions, and, 2. the four-field approach, where each numerical phase has a fixed representation and then adapts to the local morphology via being either present or absent as a results of explicitly defined morphology transfers.
Within the two-field approach each physical phase, for example, gas or liquid, is represented by exactly one numerical phase, that is, the number of numerical equals the number of physical phases. Hence, each numerical phase represents resolved as well as sub-grid morphologies, depending on an algebraic blending method to detect the individual morphologies. All known approaches require a blending formulation based on threshold values for absolute values as well as gradients of the void fractions of the phases in combination with a shape function, which might be of linear, hyperbolic, or polynomial type. The formulation of blending coefficients for systems of three or more phases quickly becomes very complex, such as the one in the OpenFOAM Foundation release 11 version 10 and higher. Appropriate scale-dependent interfacial closure models are applied via multiplication of the individual models with a blending factor. Thus, the inter-scale transfer of phases takes place inside one single numerical phase simply by switching the closure formulation by means of changing blending factors. Hence, such transfers also do take place due to deficiencies that are inherent to the numerical method, such as numerical diffusion of the interface modifying the effective blending behavior or the interface compression mechanism. 11 Reliable criteria for inter-scale transfers are also hard to determine. For instance quantities like interface curvature have been used for that purpose by Coste et al. 12 and by Wardle and Weller. 13 Overall, the two-field approach appears to be quite robust, but not particularly accurate and it has been used for various applications. 1,2,10,[12][13][14][15][16] At this point, the recent work of De Santis et al. 10 is discussed in more detail, which is based on the approach of Wardle and Weller. 13 The morphology is adapted based on a binary switch that is coupled to the blending framework implemented in the OpenFOAM Foundation release. That particular blending approach combines the following two closure formulations: (a) the Euler-Euler model including interfacial closures for interfacial momentum transfer for describing the sub-grid scale dynamics and (b) the large-interface formulation according to Marschall 17 for describing large-scale resolved interfaces together with the continuous surface force approach of Brackbill et al. 18 for modeling of the surface tension. The binary morphology switch takes a value of 0 for dispersed and a value of 1 for large-scale interface morphologies. A large-scale interface is indicated, if the following three conditions are met: • the interface resolution quality index IRQ 19 based on interfacial curvature exceeds a given critical value, and • the gas void fraction takes values between 1% and 99%, and • the dispersed bubble or droplet diameter exceeds the local grid spacing by a user-defined factor.
The approach is tested for various multiphase scenarios: an ethanol droplet in equilibrium, a dam break, a rising gas bubble, a dispersed flow in a vertical pipe, and a plunging liquid jet. The model is further verified for stratified and slug flow regimes by Colombo et al. 14 Overall a significant over-estimation of the interfacial shear stress is observed, which is attributed to the interfacial momentum transfer modeling. Especially, the interfacial drag modeling at resolved interfaces is the key to high accuracy in that context. More precisely, the interfacial drag formulation needs to ensure the kinematic condition at the interface, that is, interface normal components of gas and liquid velocity need to be identical, and at the same time provide an appropriate amount of friction between both phases for the transfer of momentum directed parallel to the interface. Overall the high complexity of the algebraic blending in connection with the momentum transfer at resolved interfaces makes it difficult to derive generalized closure models in the context of the two-field approach.
A more promising way to control momentum transfer at a resolved interface and, thereby, to improve the accuracy is the four-field approach. The main goal of this method is to gain explicit control over which features of the interfacial flow become resolved directly and which are modeled on a sub-grid scale. 8 This approach can be characterized by the number of numerical and physical phases. The latter quantity is simply defined by the problem statement, for example, a system of air and water contains two physical phases. Each physical phase might be represented by several numerical phases, that are each dedicated to a single morphology, for example, dispersed gas (bubbles), continuous gas, dispersed liquid (droplets), and continuous liquid. A set of balance equations for mass, momentum, energy, and so on, are then defined for each individual numerical phase. The interfacial closures between the pairs of phases are formulated based on the underlying morphology of the numerical phase. In other words, each numerical phase is dedicated to a specific morphology of a specific physical phase and, hence, no algebraic blending expression is required anymore. As the transition between morphologies is not inherent to the modeling approach, morphology transfers require explicit modeling strategies via so-called phase transfer models, some of which are proposed by Frederix et al. 8 Due to the increased number of transport equations, hybrid models based on the four-field approach are more computationally demanding compared to the blending based two-field approach. For the sake of direct comparison, the results for the edge case of a fully resolved single gas bubble rising in stagnant water 20 are presented in Figure 1. The two-field based model of De Santis et al. 10 is compared to the four-field based model of Meller et al. 9 Furthermore, results obtained by Hysing et al. 20 with a level-set method and the ones obtained by Klostermann et al. 21 with an algebraic VOF method are shown as references here. It becomes clear that the four-field approach (Meller et al. 9 ) that utilizes a numerical closure formulation for interfacial drag 22 delivers a clearly defined resolved interface. Those results are very similar to the reference data, while the hybrid model based on the two-fluid approach (De Santis et al. 10 ) predicts the rising velocity of the gas bubble to be too low. In addition to the interfacial momentum coupling, the four-field approach allows for interaction of large-scale interfaces with dispersed phases of an arbitrary number. A comprehensive verification of a four-field model is carried out by Frederix et al., 8 who apply the method to dispersed bubbly flows, fully resolved single bubble flows, a rising bubble including morphology transfer as well as to two-phase horizontal pipe flows of several flow regimes. They mainly use interfacial closures for F I G U R E 1 Comparison to different hybrid multiphase models 9,10,20,21 in two-dimensional benchmark case of a rising gas bubble according to case 1 of Hysing et al. 20  dispersed flows and for surface tension at large-scale interfaces similar to De Santis et al. 10 However, Frederix et al. 8 replace the interfacial drag formulation for large-scale interfaces with the correlation of Schiller and Naumann 23 and define an artificial value for the diameter in the correlation. It is stated, that an interfacial drag coefficient larger than zero is generally required to stabilize the numerical method in case of vanishing void fraction and, additionally, needs to be applied to resolved large-scale interface. 13 Such a drag closure formulation does not account for the structure of the relative velocity field between the two phases in the vicinity of the interface. Furthermore, the work of Frederix et al. 8 focuses on inter-scale transfers, which are also investigated by Yin et al. 3 However, inter-scale transfers are beyond the scope of the present work. Meller et al. 9 verify the improvement of the four-field approach in accuracy by strong interfacial drag coupling for high spatial resolutions via application to the single rising gas bubble as presented in Figure 1. For fine grids bubble diameter and shape are in accordance to reference data in the literature. However, on very coarse grids, that is, four computational grid cell per bubble diameter or less, the approach fails to preserve the bubble as a compact gas structure in the two-dimensional case described above and, hence, fails to predict bubble rising velocity. A similar dependence of bubble rising velocity on the grid resolution is reported by Frederix et al. 8 Hence, the prediction of rising velocity and shape of gas bubbles generally requires sufficiently high mesh resolution. Typically, a number of 20 cells per equivalent bubble diameter or higher is necessary to correctly predict the dynamics of the gas bubble. 24,25 Nevertheless, a well-defined behavior of resolved interfacial structures is desired on all levels of grid resolution. Such a property of a multiphase model is important for the formulation of morphology adaptive models as situations of insufficient grid resolution are inevitable when coupling scale resolving multiphase models to models based on phase-averaging. It is the belief of the authors, that the shape prediction with coarse grids plays a minor role in case purely overall phase distributions and statistical averaged behavior are of interest as long as interfacial structure (bubbles, droplets) are staying compact and are kept together, that is, unintended numerically induced breakup is successfully prevented.
The present work focuses on the control of interfacial slip on coarse grids in order to improve the prediction of the momentum transfer across an interface, for example, to improve prediction of bubble rising velocity. The four-field model of Meller et al. 9 forms the conceptual basis of the present work and is referred to as MultiMorph model hereafter. The main design criteria are presented by Schlegel et al. 26 The Euler-Euler (E-E) model based on the two-fluid formulation 27 is adopted for small-scale interfacial structures and a VOF-like approach is used to track large-scale interfacial structures. For both edge cases of very low and very high spatial resolution the individual base models are validated and proven to work as intended by Rzehak et al. 28 and Meller et al., 9 respectively. The model performance for the prediction of resolved stratified flows with large-scale horizontal interfaces is assessed by Tekavčič et al. 29 Contrary to Frederix et al., 8 the present work attempts the formulation of a well-defined interfacial drag for large-scale interfaces tracked with the computational grid of given resolution by accounting for the local velocity field. This is achieved by imposing an interfacial no-slip condition, that is, individual phases stick together and, hence, move with an identical velocity. That capability is indicated by the simulation result of a three-dimensional gas bubble in a skirted regime in Figure 2 according to regime II as defined by Frederix et al. 8 Even with a coarse grid resolution of eight grid cells per initial bubble diameter the bubble is preserved as a compact gas structure. Apart from that in the Euler-Euler mode of the hybrid model phase-specific velocity fields generally differ from each other, that is, interfacial slip takes place. This duality of the hybrid model suggests that the prediction of interface dynamics in gas-liquid flow with intermediate spatial resolutions requires the amount of interfacial slip to be locally controlled depending on the local degree of resolution. While this should be achieved in flows directed tangential to the interface surface via adjustment of interfacial momentum exchange, phase specific velocities still need to be identical in flows directed orthogonal to the interface in order to fulfil the kinematic interface condition. This idea is illustrated in the example of a poorly resolved rising gas bubble, which is sketched in Figure 3. On the left hand side, the bubble dimension is small compared to the grid spacing, which is described with the ensemble-average E-E model and, consequently, gas and liquid phases slip along each other across the whole volume of the gas structure. Conversely, on the right hand side grid cells are much smaller compared to the gas bubble. Thus, the VOF-like mode is applied across the whole bubble volume imposing the interfacial no-slip condition everywhere. For the range of intermediate grid resolution a mixed approach is suggested. The lateral parts of the bubble interface are associated with an interfacial shear layer and should therefore experience interfacial slip, that is, locally phase-averaged gas and liquid velocities should differ from each other. On the top and bottom locations of the bubble interfacial flow is expected to be of stagnant type. In those regions a no-slip condition is maintained between gas and liquid in order to achieve a clear interface definition.
After presenting the basic numerical method in Section 2, the present work is focused on the classification of the interface-surrounding flow. Based on that a measure for the degree of spatial under-resolution of interface-parallel shear flow is proposed in Section 3. This quantity is considered analogue to the well-known dimensionless wall-normal coordinate in shear boundary layers of wall-bounded flows. 30 Subsequently, this measure is utilized to control the amount of interfacial slip by means of a dedicated drag model for the interfacial momentum exchange, which is presented in Section 4. In Section 5, the MultiMorph model is then validated in the case of a three-dimensional gas bubble rising in stagnant liquid according to Adelsberger et al. 31 Furthermore, in Section 6 the hybrid model is assessed in the case of a two-dimensional stratified co-current gas-liquid flow with a flat interface, which is experimentally investigated by Fabre et al. 32 In that context, interfacial slip allows for strong velocity oscillations in vertical direction in the interface region, which is handled with a modified formulation of the buoyancy term in the two-fluid model. Finally, in Section 7 the findings of this work are concluded and perspectives for future investigations are presented.

NUMERICAL METHOD
The hybrid multiphase model is based on the morphology-adaptive multifield two-fluid model MultiMorph. 9 The term two-fluid underlines the fact, that this work is based on the assumption of two (or more) interpenetrating continua. Individual physical phases might be described as a set of different numerical phases, which is expressed by the term multifield. More precisely, physical phases are characterized by certain fluid properties, while numerical phases additionally have a certain flow morphology as a property, such as continuous or dispersed flow. For instance, gas is either described as a continuous gas or as a dispersed gas, that is, gas bubbles. The phase-averaged Navier-Stokes-Equations 33 consist of a continuity and a momentum balance equation for each individual phase , The Einstein summation convention applies for Latin indices. A partial derivative with respect to time t or spatial direction x i is denoted ∕ t or ∕ x i , respectively. The phase-specific volume fraction, density, and viscosity are r , = const, and = const, respectively. The phase-averaged velocity vector is u ,i , p is the pressure, which is shared between all phases, and g i denotes the vector of gravity. The symbols , , and n ,i denote the surface tension coefficient of the phase pair and , the curvature, and the normal vector of an interface, respectively. The phase-specific shear-rate tensor The vector of interfacial momentum exchange is referred to as f ,i . The equations are spatially discretized with a second order finite volume method for unstructured grids and integrated in time with a semi-implicit Euler-scheme of first order accuracy. An interface compression term 34 is applied to the phase-fraction transport equations. Pressure-velocity coupling is realized via the compact momentum interpolation method. 35 The resolution of the stiff system of equations due to potentially large interfacial coupling forces is realized by means of an approximate formulation of the partial elimination algorithm 36 according to Meller et al. 9 The hybrid model is implemented in the multiphaseEulerFoam framework of the C++ library OpenFOAM Foundation release. 11 The hybrid multiphase model MultiMorph is available under GPL license. 37

DETECTION OF UNDER-RESOLVED INTERFACE REGIONS
The basic idea is to relate the near interface shear flow to the local thickness of the interface representation in order to obtain a local measure for grid resolution. It is inherent to the underlying algebraic VOF method, that the interface thickness strongly correlates with the local size of the computational grid cells. The local shear is determined in terms of the relative phase-specific velocity and of the interface orientation. The length scale of the interface thickness is then related to a dimensionless shear-based length scale in the interface region, similarly to the dimensionless wall distance y + = y∕L in wall-bounded flows with wall-normal coordinate y and shear length scale L . 30 Such a measure has been proposed by Coste 38 before along with a corresponding procedure to calculate the interfacial drag. Within that approach it is assumed, that the interface consists of exactly three cells and the interface thickness is separately calculated for both sides of the interface. Furthermore, Coste 38 applied the interfacial drag coupling in a tensorial manner by imposing different drag coefficients in interface normal direction and in the direction tangential to the interface, which is not compatible with the compact momentum interpolation method used in the present work. Hence, the procedure presented in the following differs from the approach of Coste 38 by the following: • no explicit distinction is made between gas and liquid sides of the interface, • it is not limited to an interface thickness of exactly three cells, which is likely to be exceeded in case of phase transitions in the hybrid model, • the drag formulation adapts to the local flow morphology rather than adapting to the phase fraction distribution and • the drag formulation is a scalar formulation rather than a tensorial one.
In order to formulate the interface thickness measure, three different pieces of information are required: (1.) the phase-specific velocities u and (2.) u as well as (3.) the interface vector, i = r ∇r − r ∇r . A schematic of the interface including the locations of all three vector quantities mentioned before is shown in Figure 4. Note that the interface vector is not normalized and, hence, has the unit of inverse length.
Both regions that contain the pure phase or are denoted A and B, respectively. They are separated by the interface region with 0 < {r , r } < 1 referred to as region C. As typical for an algebraic volume-of-fluid(-like) method, that is, one edge case of the utilized hybrid model, the interface is smeared out across several grid cells. Therefore, region C has a finite thickness. Due to the underlying phase-averaging procedure, 33 the three pieces of information mentioned above are meaningful at different individual locations. Namely, phase velocity u is of significance inside phase (region A), indicated by phase fraction r , and phase velocity u inside phase (region B), indicated by phase fraction r . The interface vector i is of significance in the interface region C. At the same time its magnitude, |i |, serves as an indicator for the interface region. In order to make the required information available all across the interface at a certain interface position including adjacent locations, a transport of the values of all three quantities is required. This is achieved via an iterative procedure, by which the information are transported step-wise across neighboring grid cells along the negative gradient of each individual distribution of meaningfulness, that is, −∇r (for u ), −∇r (for u ), and −∇|i | (for i ). A constant number of N I iterations is carried out, which defines the maximum accepted number of grid cells across the interface due to numerical smearing. In that way, all three pieces of information are transported form the location of their origin, that is, the region of highest meaningfulness, to the locations of the interface, where the information is missing, that is, the region of low meaningfulness. The procedure is visualized in Figure 5.
In order to achieve a robust algorithm in cases of two-or three-dimensional space, special care has to be taken. Prior to the assessment of the transport direction of the individual quantities via the gradients of individual indicators, namely r , r , and |i |, these three field variables are smoothed by simply interpolating the field data from cell centers to cell face centers and back again. Furthermore, if values of an interface quantity are transported into a grid cell from multiple cell neighbors, the information originating from the location with the higher respective meaningfulness is preferred.
The quantities resulting from the presented algorithm-either directly or via derived quantities-are denoted with superscript I.
Finally, the same meaningful value for each quantity of interest is available at each cell in the interface region. Note, that the procedure strictly avoids to gather velocity information from within the interface region. Instead, it clearly determines values related to positions very close to both sides of the interface. The region within the interface is the region to be controlled by appropriate slip formulations.
First of all, the flow surrounding an interface is classified by its direction relative to the interface normal vector i . One of the two extreme cases illustrated in Figure 6 is stagnation flow, in which the relative velocity between both fluids is directed perpendicular to the interface. The other extreme case is shear flow, where both liquids have a relative flow velocity parallel to the interface. The degree of shear is expressed as the shear flow indicator: In case either factor in the denominator takes the value zero, the shear flow indicator is set to I = 0. Besides pure stagnation flow such a situation is faced in regions, which are not located in the vicinity of the gas-liquid interface and, in consequence, no interfacial shear flow is predicted there. The indicator takes values between zero and one for flow types other than stagnation or shear, depending on the angle between relative velocity and interface normal vector. At this point it is stressed, that the degree of shear, I , is based on the velocity information transported with the algorithm presented before, while the information on the interface direction is local. In this way, a smooth distribution of I is achieved.
In order to estimate the degree of spatial under-resolution, the dimensionless interface thickness is defined analogously to the theory of boundary layers in wall-bounded flows: 30 This quantity represents the ratio of the local interface thickness Δ I to the local shear length scale L * I . Assuming that the phase volume fractions are distributed linearly across the interface and that the interface region reaches all the way from r = 0 to r = 1, the inverse of |i I | gives an estimate for the interface thickness: The distribution of shear stress in a linear shear flow along coordinate y can be expressed as (y) = u∕ y with dynamic viscosity and velocity component in flow direction u. Based on Equation (5), the shear stress across the whole interface thickness can be estimated as regardless of the distribution of individual phase-specific velocities inside the interface region. Using the shear stress approximation, an interfacial shear velocity is defined as F I G U R E 7 Form function for under-resolution in case of a pure shear interface, that is, I = 1.
Analogous to the viscous length scale in wall bounded flows, the shear length scale L * I is calculated as The material properties of the phase pair , namely densitỹ, kinematic and dynamic viscosity,̃and̃, are calculated as arithmetic means of the individual phase-specific values as̃= ( + )∕2. The arithmetic mean value of fluid properties is used instead of a phase-weighted mean value for instance, because the mixture quantities are required to be represented by constant values in order to achieve an integral shear length scale for the whole interface, while the composition of the fluid mixture changes across the interface region. As for wall-bounded flows, it is assumed, that a dimensionless interfacial thickness y +I < 1 indicates a complete resolution of the shear boundary layer. As mentioned before, the expected interface width is not larger than N I grid cells, meaning, that the interfacial velocity gradient is captured on the computational grid in the aforementioned case. In the opposite case, if the quantity y +I is larger than unity, under-resolution is assumed. As for wall-bounded flows, y +I is a measure for the ratio of grid spacing to length scale of velocity gradients.
Based on both the shear flow indicator I and dimensionless interfacial cell size y +I , the under-resolution indicator is proposed as The function is illustrated in Figure 7. This quantity becomes effectively one in case of a fully under-resolved shear flow, that is, y +I ≫ 1, at the respective interface position. If the flow is either of stagnation type or is considered fully resolved, that is, y +I < 1, the indicator takes the value zero. For other cases UR takes values between zero and one.

ADAPTIVE DRAG MODELING
Based on the proposed criterion for detecting under-resolved interface regions, an adaptive drag formulation is derived. The interfacial drag force is generally formulated as In order to achieve tight coupling of phase-specific velocity fields for resolved interfaces (volume-of-fluid-like mode of the hybrid model), the interfacial drag formulation of Štrubelj and Tiselj 22 is adopted. This drag model is referred to as resolving drag model, which is indicated by the superscript R for the drag coefficient K D,R , which is defined as The phase-fraction weighted density of the phase pair is formulated such that it is applicable to situations with more than two phases present. The relaxation time is chosen as r = 10 −8 Δ t with physical time step Δ t , which generally leads to very high values of the resolving drag coefficient K D,R . Hence, the relative velocity between phases and becomes negligibly small.
Considering a gas bubble, which is much smaller than the computational grid spacing, the Euler-Euler model is applied in the frame of the used hybrid approach. However, if the gas structure is of similar dimension compared to the grid spacing, situations arise in which the whole bubble is still depicted on the computational grid, but the dynamics in the region around the gas-liquid interface cannot be captured with the given grid resolution. In that case, the local interface is treated under-resolved by allowing for local interfacial slip. That is achieved with the under-resolving (UR) drag model formulation K D,UR considering a drag coefficient C d : 39 Gauss et al. 40 propose to set the model constant to C d = 0.22 in order to reproduce consistent bubble rising velocities on coarse to very coarse grids. With that value for the model parameter in combination with the hybrid model formulation under investigation, this turns out to cause disintegration of a gas bubble on coarse computational grids. A value of C d = 0.8 has been determined in preliminary investigations by manual calibration and has shown to be better suited for this kind of problem and, hence, is used hereafter. From that drag coefficient follows, that K D,UR takes moderate values, which in turn allows for interfacial slip, that is, a relative velocity between phases in the interface region.
As the values of both drag formulations named above differ by several magnitudes for K D , a simple linear interpolation appears to be inappropriate. The harmonic average is used to obtain the resulting drag coefficient: This formulation effectively blends between the resolving drag formulation for UR = 0 and the under-resolving one for UR = 1. The distribution of K D over UR for exemplary values of K D,UR = 1 and K D,R = 10 4 is shown in Figure 8. The proposed formulation is referred to as resolution adaptive drag model.

THREE-DIMENSIONAL RISING GAS BUBBLE
In order to assess the framework of interface classification and resolution adaptive drag modeling as described before, a three-dimensional test case is selected, which is proposed by Adelsberger et al. 31 It features a three-dimensional single gas (G) bubble initialized as sphere with diameter D b = 0.5 m rising in stagnant liquid (L) under the influence of gravity. The computational domain has a cuboid shape with dimensions of 2D b × 4D b × 2D b in the three respective spatial directions x × y × z with gravity being directed in negative y direction. The origin of the coordinate system is located in the lower left rear corner of the domain. Initially, gas and liquid phases are at rest and the gas bubble is located at (1D b , 1D b , 1D b ), that is, in the centre of the domain one bubble diameter above the lower boundary. All boundaries are considered to be of Dirichlet type for phase-specific velocity fields with value zero. For pressure and phase fraction fields all boundaries are Neumann conditions imposing zero derivative. The selected test case is originally denoted with Case 2 31 and is characterized by the following dimensionless numbers: The ratios of density and dynamic viscosity between liquid and gas are 1000 and 100, respectively. The Eötvös number The domain is spatially discretized with a regular computational grid with cubical cells ( Table 1). The numerical setup is reported by Hänsch et al. 41 In order to assess the performance of the hybrid model, five different grids are used, namely G1 to G5. Detailed information about the different grids are listed in Table 1, which includes the ratio of cell size Δ x to sphere-equivalent bubble diameter D b besides the number of grid cells in each direction.
Assuming the proposed algorithm for the detection of local under-resolution delivers reasonable results, the following model behavior is expected: The reference solution can be recovered with the resolving drag model on the finest grid resolution (G5), which demonstrates the equivalence of the hybrid model approach to recover the homogeneous model behavior in case of high spatial resolution. The results with the resolving drag model obtained on coarser grids (G1 to G4) will reveal wrong bubble dynamics, that is, the gas bubble is rising too slow on coarse grids. Application of the under-resolving drag model will lead to a faster bubble rising velocity exceeding the reference values. The resolution adaptive drag model will lead to more consistent bubble dynamics for all grid resolutions.

Recovering homogeneous model behavior
Adelsberger et al. 31 obtained the reference data on grid G5. With the current approach, the reference data are reproduced with the resolving drag model on the same grid. The results in terms of vertical bubble velocity over time as well as interface location at time t = 3.5 s are presented in Figure 9. The reference lacks information regarding the scale and absolute position of the gas-liquid interface shown in Figure 9B. Hence, those data of interface location are scaled equally in all directions and are then translated such that they match the current results as close as possible. Therefore, the reference data for the interface position purely serves for comparison of the shape of the bubble. It turns out, that reference results are reproduced with the resolving drag model and with an identical grid resolution considering the bubble shape and rising velocity. That demonstrates the equivalence of the hybrid model approach to the homogeneous model in case of high spatial resolutions with the according drag model in this particular case. This observation is in line with the findings of Yan and Che, 42 who have mathematically proven that two-fluid model and homogeneous model are generally equivalent in the boundary case of an infinitely large interfacial drag coupling.

Resolving drag model
The results obtained with the resolving drag model according to Equation (11) are presented in Figure 10. Contrary to the previous section, the simulation is carried out on several coarser grids (G1 to G4 in Table 1) in order to assess the influence of the spatial resolution. As evident from the reference data, the bubble initially accelerates from rest and reaches a temporary maximum rising velocity before maintaining a nearly constant velocity with a slightly lower value. Finally the bubble starts to decelerate for t > 3 s, as it approaches the upper wall. The bubble rising velocity predicted with grid G4 is nearly identical to the reference data, while a strong underestimation is observed on the coarse grids, especially with G1. The temporary velocity maximum is reached too late and the value is 5.7% and 18.9% lower with G2 and G1, respectively. Subsequently, the gas bubble continues decelerating further on these two grids instead of rising with nearly constant velocity.
Considering the bubble shape in Figure 10B, the result obtained with G4 is very close to the reference results, except for the tips of both ligaments at the bottom of the bubble being slightly less sharp. Grid G3 reveals a more narrow bubble shape with more elongated ligaments and the bubble reaches a slightly lower vertical position compared to G4. A lower vertical position is reached at t = 3.5 s on the coarse grids G1 and G2 compared to G4, which results directly from the reduced rising velocity as discussed before. On grid G2, the bubble is deformed in such a way, that its shape is narrower and more bend, resulting in a curved shape of the interface at the centre bottom of the gas bubble. The bubble shape obtained with G1 differs even more from the reference data as it shows a peak at the top of the gas structure. The smaller apparent size of the bubble on the coarse grid G1 is explained with the strong smearing of the interface, that is, gas diffuses into the liquid, due to the low spatial resolution, while the total volume of the gas is conserved. Therefore, the resolving drag model is assessed to be a reasonable approach for high spatial resolution, which allows to reproduce simulation results of the homogeneous model. 9 On coarse computational grids, this drag formulation delivers a bubble rising velocity, which is too slow, while the gas bubble is deformed in a nonphysical way.

5.3
Under-resolving drag model As pointed out in Section 4, the drag formulation according to Equation (12) with a constant coefficient C d with a value of 0.8 is utilised for cases, in which an interfacial slip velocity is allowed across the volume of the whole gas bubble. A grid refinement study with grids according to Table 1 is carried out with the under-resolving drag formulation. The results are shown in Figure 11. Considering the bubble rising velocity in Figure 11A, it is evident, that the bubble is predicted to rise too fast on coarse grids, especially G1, compared to the reference data. With G1, the maximum vertical velocity is overestimated by 40.8%. In the following period, the bubble rising velocity is oscillating and is higher than the reference value until t = 3.1 s. Subsequently, the velocity drops rapidly, because the gas bubble reaches the top boundary of the computational domain too early. With grid G2, the bubble decelerates too much after reaching its maximum velocity. In contrast to the resolving drag model, on the fine grid G4, the rising velocity is too high compared to the reference data. The velocity is overestimated by 4.2% at its maximum and does not reach the reference value until t = 3.5 s.
On grid G1 the gas bubble reaches the top wall (y = 2 m) before t = 3.5 s, which is in line with the strong over-prediction of bubble rising velocity shown in Figure 11A. Hence, in that case the interface is located outside the visible range in Figure 11B. On all other computational grids under investigation, the bubble is too flat compared to the reference data and no sharp ligaments are observed at all. Hence, it is assumed, that with the under-resolving drag formulation the convergence rate of the results towards the reference values is extremely low, if the latter are even achieved at any level of spatial resolution.

Resolution adaptive drag model
The interface classification approach presented in Section 3 and the resolving drag formulation described in Section 4 are assessed in combination with the same grid refinements study. The resulting under-resolution indicator UR GL is shown in Figure 12.
It is evident, that the largest degree of under-resolution is detected on grid G1 at the lateral positions of the gas bubble with a maximum value of approximately 0.7. This is expected as shear flow exists in this region of the interface, which is shown schematically in Figure 6B. Contrarily stagnation flow (see Figure 6A) is observed on the centre top and bottom locations of the interface, which results in a lower under-resolution indicator UR GL . With increasing spatial resolution, the region of detected under-resolution and the maximum predicted value of UR GL become smaller. On G4 full resolution is The grid study results are presented in Figure 13. The reference data in Figure 13B has the same scaling and offset as in Figure 10B.
The rising velocity in Figure 13A as predicted with G1 still overshoots the reference data by up to 21.1% before reaching its temporary maximum. But compared to the velocity over-estimation of up to 40.8% obtained with the under-resolving drag formulation in the previous section the error in the prediction on the coarse grid G1 is nearly reduced by half its value. This is a remarkable improvement, especially considering the extremely low grid resolution of four grid cells per sphere-equivalent bubble diameter. Further improvements are expected if the drag coefficient may be modelled in a more adaptive fashion in the future rather than taking a constant value of C d = 0.8. After reaching the highest rising velocity, with G1 the bubble decelerates and the deviation from the reference data shrinks. With grid G2, a rising velocity is predicted, which is much closer to the reference data than with both the resolving and the under-resolving drag formulations with the same spatial resolution. An improvement in terms of bubble rising velocity, shape and interface position is compared to the basic drag formulations is also observed on G3. The results with G4 are quasi identical to the resolving drag formulation on G4 and, hence, almost no error is observed when compared to the reference data.
Considering the bubble shape in Figure 13B, the solution obtained with G1 reveals, that the gas bubble is located slightly ahead of the position obtained with G4, but the difference is minor compared to result G1 with the under-resolving drag formulation. The bubble shape obtained with this grid is flat and shows no ligaments. It is worth noting, that grid G1 corresponds to a spatial resolution of four grid cells per sphere-equivalent bubble diameter, thus it is impossible to capture more complex shapes on such a grid. Similarly to Figure 10B, the volume bounded by the interface according to definition G = 0.5 appears to be smaller with G1 compared to finer grids due to the smeared interface, while the total gas volume is conserved. For grids G2 to G4, convergent bubble shapes are observed resulting in minor deviations from the reference solution with G4.
In this test case, the resolution adaptive drag modeling framework delivers results, which are reasonable across all levels of grid refinement under investigation improving especially the solution on very coarse computational grids. In contrast, the under-resolving drag formulation with a fixed drag coefficient all across the computational domain turns out to deliver unreliable results for all levels of spatial resolution. For fine grids, the solution of the resolution adaptive drag model converges towards the volume-of-fluid-like solution obtained with the resolving drag formulation and, hence, is quite close to the reference data. This is confirmed by comparing the vertical position of bubble centre of gravity at t = 3.5 s over spatial resolution for different interfacial drag modeling approaches, which is presented in Figure 14.
The under-resolving and the resolving drag model approaches show over and under prediction of the vertical position at the time instance, respectively. The resolution adaptive drag model approach delivers a in-between result, which is less depending on the spatial resolution, hence, contributing to more reliable results on coarse computational grids. At the same time, the results on fine grids are equivalent to the ones obtained with the resolving drag model. That behavior is close to an ideal mesh independent behavior of the resolution adaptive drag model, which is represented by the horizontal line, that is, the final position of the gas bubble would be Y b ∕D b (t = 3.5 s) = 3.17 for all grid resolutions.

TWO-DIMENSIONAL CO-CURRENT STRATIFIED CHANNEL FLOW
The second test case considers a horizontal co-current stratified air-water flow in a rectangular channel according to the experiment of Fabre et al. 32 Both air and water flows are fully turbulent. The case under investigation refers to case 250 of the reference, which is defined by volume flow rates of water and air of 0.003 m 3 s −1 and 0.0454 m 3 s −1 , respectively. In the experiments that results in a non-wavy smooth water surface with a measured mean water depth of 38 mm. Measured bulk velocities for both phases are U water = 0.395 m s −1 and U air = 3.66 m s −1 . That results in channel Reynolds numbers of Re water ≈ 16,800 and Re air ≈ 14,500 based on the respective height of the channel section of each fluid. Fluid properties are presented in Table 2.
The two-dimensional computational domain for the numerical setup features a rectangular channel section with a length of 12 m and a height of 0.1 m. Upstream of the rectangular channel, two inlet regions of 1 m length each for both air and water are attached with the water section being 38 mm high according to the measured water level. Both inlet regions are separated by a thin baffle, which is treated as no-slip wall for both fluids. In order to achieve fully developed turbulent channel flows in each individual inlet section, instantaneous profiles of flow quantities are mapped from the TA B L E 2 Fluid properties of water and air for case of two-dimensional co-current stratified channel flow. end of the respective inlet section back to its beginning. The geometry of the computational domain along with boundary and initial conditions is presented in Figure 15. For each individual phase the k − SST model 43 is applied including asymmetric turbulence damping (in the air only) in the interface region according to Reference 29. Turbulent wall functions 44 are specified as boundary conditions for both, turbulent kinetic energy k and turbulent specific dissipation rate , for the individual phases. Vertical data profiles are extracted at position x = 9.1 m downstream of the tip of the baffle in order to be compared to experimental data. 32 Spatial discretization is realised via an orthogonal grid with equidistant spacing in streamwise direction. In the vertical direction, grid spacing is constant within the individual sections for air and water. The characteristics of different computational grids are presented in Table 3.
In order to compare the results for the stratified air-water flow, profiles of velocity u and of turbulent kinetic energy k are presented as mixture quantities unless specified otherwise. All mixture quantities are calculated from the phase-specific quantities : Similarly to Section 5, the case is initially investigated with the resolving and with the under-resolving drag formulations. As it turns out, a modification of the formulation of the buoyancy term is required, when applying the under-resolving drag formulation either exclusively or in the context of the resolution adaptive drag formulation. This modification will be applied to the under-resolving drag model setup before assessing the performance of the resolution adaptive drag model together with the modified buoyancy term formulation.

Resolving drag model
The results achieved with the resolving drag model (see Section 4, Equation (11)) are presented in Figure 16.
The streamwise mixture velocity profile in Figure 16A shows rather large deviations in the air section among various grid spacings. On the coarse grids G1 and G2, the velocity close to the upper channel wall is over-predicted compared to experimental data. Close to the interface, the velocity gradients are predicted to be too small, most likely as a consequence of the large grid spacing. With an increasing grid resolution, the velocity gradient in the air section just above the interface becomes larger, while the maximum velocity slightly decreases. That leads to an under-prediction of maximum velocity in the air section, especially on fine computational grids. With all grids the vertical position of the maximum velocity is predicted to be too high compared to experimental data.
The streamwise water velocity in the lower section of the channel is presented in Figure 16B. While the water velocity is over-predicted on the coarse grids, the profile is similar on fine grids G4 to G6. With grid G6 the experimental data are reproduced well.
The profile of vertical mixture velocity component is presented in Figure 16C. Compared to negative experimentally measured values of v ∈ (−3.8 × 10 −3 m s −1 , −1.3 × 10 −3 m s −1 ) the vertical water velocity predicted by the numerical model is negligibly small. In the air section vertical velocity contributions are larger compared to most experimentally obtained data but they are very small compared to both experimental outlier values of v = 1 × 10 −2 m s −1 . For grids G1 and G4 the vertical air velocity is negative with small values while it takes slightly larger positive values for all other grids. In terms of vertical velocity no mesh convergent behavior is observed.
The profile of mixture turbulent kinetic energy k is shown in Figure 16D. The qualitative behavior observed in the experimental results is captured with the numerical method across all grids. Generally, k is decreasing with increasing wall distance in both air and water section at the very top and at the very bottom of the channel. The turbulent kinetic energy takes a local minimum in the air section, that is, boundary layers resulting influences of both upper wall and air-side of the interface meet at that point. The vertical position of the extremum coincides with the local maximum of horizontal mixture velocity u and, hence, is also located too high across all grids. Right above the gas-liquid interface quantity k is over-predicted, especially on the coarse grids G1 to G3. Generally, the profile of turbulent kinetic energy reveals a mesh convergent behavior. However, inside the water flow an under-prediction is observed with all grids and the deviation from experimental measurements increases with higher grid resolution.

6.2
Under-resolving drag model The results achieved with the under-resolving drag formulation according to Equation (12) with a constant coefficient C d with a value of 0.8 are presented in Figure 17.
The streamwise velocity profiles in Figure 17A reveal sharper mixture velocity gradients close to the interface compared to results obtained with the resolving drag model in Figure 16A. This is explained with the slip between phase-specific velocities at the interface, which allows less momentum transfer between air and water across the interface due to the chosen interfacial drag formulation. Furthermore, a strong over-prediction of the air velocity is observed for coarse grid G1. The value of the maximum air velocity decreases with finer grid resolution, while its vertical position is again too high compared to experiments. On medium and fine grids the velocity profiles in the air section of the channel show the maximum velocity at a lower vertical position compared to resolving drag model results.
From Figure 17B it becomes clear, that the streamwise water velocity is nearly constant across the channel height for coarse grids. The value of this quantity is heavily under-predicted on the coarse grids, especially with grid G1. This behavior can be explained with the vertical position of the gas-liquid interface, which is indicated by the kink in the streamwise velocity profiles, being located at a much higher position for the coarse grids compared to finer grids. From mass conservation it follows, that such a higher water level results in a lower bulk velocity. For increasing spatial resolution the water velocity is predicted to be larger and more similar to a boundary layer profile. Finally, with finest grid G6 the deviations from experimental data are very small except in the vicinity of the interface.
Profiles of the vertical mixture velocity are presented in Figure 17C. Large negative peaks are observed for coarse grids in the vicinity of the gas-liquid interface. This effect is most pronounced on grid G3 with a minimum of −0.9m s −1 which is two orders of magnitude larger than the extremum of the experimental data. An explanation for this behavior is the gravity acting on both air and water phases via the different phase-specific densities of air and water , respectively. As a consequence, both phases are being accelerated in opposite directions locally in each cell, which is allowed by the two-fluid model. This leads to the strong peaks in vertical velocity. As the under-resolving interfacial drag formulation with constant coefficient allows for interfacial slip, those peaks are not dampened out as in the case of the resolving drag model.
Profiles of mixture turbulent kinetic energy k are presented in Figure 17D. In the air flow k is predicted correctly in a qualitative manner for grids G2 to G6 and the position of the local minimum is consistent for all medium and fine grids, that is, grids G3 to G6. However, in the water flow the predicted profiles are qualitatively wrong for all grids. Especially on the coarse grids, the mixture turbulent kinetic energy is strongly over predicted, while an increase of k towards the interface is contrary to experimental results. This behavior can be explained by the strong peaks in the vertical velocity, which erroneously add to the production of turbulence kinetic energy in the interface region. Additionally, in the interface region a negative peak is observed, which is most pronounced on coarse grids.

6.3
Under-resolving drag model with buoyancy modification In order to overcome the flaw of strong negative vertical velocity peaks, a modification of the buoyancy term is proposed here. In detail, in the buoyancy term r g i , that is, in the third term on R.H.S. of Equation (2), for each phase sharing the large scale interface the phase-specific density is replaced by mixture density , which is defined analogously to Equation (14). With that change, the proposed gravity term in Equation (2) reads: Away from the interface, the proposed formulation is effectively identical to the old one as = , wherever r = 1. This modification is purely applicable to continuous phases, that is, phases, which are never treated dispersed in the context of the Euler-Euler model. The buoyancy still acts on the mixture similar to a usual VOF method. In other words gravity effects on large scale interfaces are still covered, as the mixture density reflects local density changes. The proposed modification purely neglects buoyancy forces inside the interface region separating phases of different phase-specific densities. Both phases are not locally accelerated in opposite directions anymore. The modified buoyancy formulation is applied to the setup with under-resolving drag model with constant coefficient and results are presented in Figure 18.
The streamwise velocity profiles in Figure 18A are generally similar to the ones obtained without modification of the buoyancy term. With grid G1 the over-prediction in the air flow is diminished. With finer grids, G2 to G6, the air velocity is under-predicted similarly to the last section. The prediction of the vertical position of the air velocity maximum is improved on coarse grids, especially with G2 and G3. The velocity gradients right above the interface are slightly smaller compared to the results obtained with the under-resolved drag formulation without buoyancy modification but still larger than predicted with the resolved drag formulation.
The profiles of the streamwise water velocity in Figure 18B are qualitatively improved with the buoyancy modification, that is, in the vicinity of the interface velocity profiles reflect the nature of a wall boundary layer with all grid resolutions. With grid G6 the results are nearly identical to the simulations without buoyancy modification. That is explained by the fact, that the thickness of the interface region is proportional to the grid spacing, which is inherent to the used hybrid multiphase model leading to a comparably small influence of processes inside the interface region. However, on the coarse grids the elimination of erroneous velocity peaks, as will be shown in the following, reveals a large impact improving the prediction of streamwise water velocity.
The vertical mixture velocity presented in Figure 18C is the main focus, when assessing the proposed buoyancy modification. The proposed formulation successfully suppresses vertical velocity peaks and oscillations in the interface region. The results are qualitatively similar to the ones obtained with the resolving drag model. The largest extremum is observed for grid G2 and lies below experimental outlier values. A smooth distribution of vertical mixture velocity is observed for grids G2 to G6.
In Figure 18D it becomes clear, that the prediction mixture turbulent kinetic energy profile is strongly improved compared to the results obtained without buoyancy modification. Especially the profile in the water flow down below the gas-liquid interface shows the same qualitative behavior as observed in the experimental results. That means, that k actually decreases with larger wall distance before jumping to higher values on the gas side of the interface. This is a result of the erroneous peaks in vertical velocity (see Figure 17C) now being effectively dampened by the modified buoyancy formulation. The value of turbulent kinetic energy is still under-predicted in the water section, but results are very similar across all grid resolutions. In the water flow the results are also close to the ones obtained with the resolving drag model. In the air flow, the local minimum of k is recovered with grids G2 to G4. Compared to the resolving drag model the value of the local minimum is consistent across all grids while the vertical location of the extremum is lower on coarse grids G2 and G3 and rises with increasing spatial resolution. Solely with the coarsest grid G1 the steep jump in k cannot be recovered as the quantity rises gradually from the interface location to the upper wall of the channel. It is suspected, that the elimination of vertical velocity peaks leads to an improved prediction of turbulent kinetic energy, which in turn improves also the transport of streamwise momentum across the interface due to turbulent diffusivity.
Overall, the buoyancy modification eliminates the flaw of strong vertical velocity peaks, while also improving streamwise velocity and turbulent kinetic energy profiles.

Resolution adaptive drag model with buoyancy modification
The performance of the resolution adaptive drag formulation is assessed together with the modified buoyancy formulation proposed in the previous section. In that context the procedure for the estimation of the degree of under-resolution from Section 3 shall be further verified. For this purpose the distribution of dimensionless interface thickness y +I across the gas liquid interface on different grids (see Table 3) is visualised in Figure 19. The dimensionless interface thickness y +I only takes positive values in the interface region, which is more narrow on finer grids. Aside from the with of the interface region the quantity y +I itself scales with the grid spacing. For grid G6 a maximum value of y +I ≈ 10 is observed. Hence, the shear flow is considered under-resolved in the interface region, even with the finest grid. However, due to the formulation of the resolution adaptive drag coefficient in Equation (13) in connection with the under-resolution indicator in Equation (9) the effective interfacial drag coupling depends on the value of y +I as long as it is larger than unity.
The vertical profiles of velocity and turbulent kinetic energy resulting from the resolution adaptive drag model together with buoyancy modification are presented in Figure 20. The streamwise mixture velocity profiles presented in Figure 20A reveal, that the over-prediction in the air flow close to the upper wall is reduced or even eliminated for G1 and G2, respectively. The vertical position of that extremum is improved for fine grids compared to the resolving drag model and velocity gradients above the gas-liquid interface are sharper on coarse grids. Value and vertical position of the kink in the velocity profile at the interface location is nearly identical for all grid resolutions under investigation.
The streamwise water velocity profiles are shown in Figure 20B. Compared to the under-resolving drag model together with buoyancy modification the under-prediction of water velocity is further reduced for grids G2 to G5. The result of grid G1 is further improved, while grid G6 produces a water velocity profile, which is comparably good as achieved with the under-resolving drag model with buoyancy modification.
In Figure 20C profiles of vertical mixture velocity are presented. The highest maximum is observed with grid G3 with values less than half of the experimental outlier values. Negative values are much smaller in magnitude compared to experimental data. For grid G2 and finer smooth distributions are observed.
Profiles of turbulent kinetic energy k are shown in Figure 20D and the results are similar to the ones obtained with under-resolving drag model with buoyancy modification. In the water flow, the values are very consistent across grids G2 to G6, while the under-prediction compared to experimental results is maintained. In the air flow, the experimentally obtained distribution of k is qualitatively recovered. While a under-prediction is observed close to the upper channel wall, the values above the gas-liquid interface are rather close to experimental reference data, especially for grids G4 to G6. The value of the local minimum in the centre of the air section is nearly grid independent and its vertical position coincides with the streamwise velocity maximum.
Overall, the predicted velocity profiles obtained with resolution adaptive drag model and buoyancy modification are superior over all other results presented earlier in this section. Compared to both resolving and under-resolving drag models the results are more consistent across all grids and improved especially on coarse and very coarse grids.

CONCLUSION AND PERSPECTIVE
Based on the morphology-adaptive multifield two-fluid model MultiMorph, 9 a new procedure for the modeling of gas-liquid interfaces is proposed, which adapts to varying degrees of grid spacing. For this purpose, two characteristics are determined across the interface region: the type of flow surrounding the interface and the degree of spatial resolution of the interfacial shear layer. The derived indicator function allows detection of regions of under-resolution, which is used to adapt the drag locally, allowing for local interfacial slip. This is realized as a controlled combination of common resolving and under-resolving drag model formulations. The functionality of this framework is demonstrated for two cases: a 3D gas bubble rising in stagnant liquid and a 2D co-current stratified channel flow. In the context of the latter case, a modification for the formulation of the buoyancy is proposed, which prevents erroneous peaks in vertical velocity in case of interfacial slip in horizontal interfaces. It is shown, that solutions converge towards reference data if the grid is refined. It is remarkable that the present approach predicts bubble rising velocity and shape as well as velocity and turbulent kinetic energy profiles reasonably well for coarse grid resolutions. These contributions are integrated in the MultiMorph model allowing more reliable predictions of gas-liquid flows by means of arbitrary spatial resolutions.
In the future, the applicability of the proposed buoyancy modification to flows with interface orientations other than horizontal ones shall be investigated. Moreover, performance and robustness of the proposed method need to be assessed for non-orthogonal, non-uniform, and non-hexahedral computational grids and improved, if necessary. Regarding the dimensionless interfacial grid spacing y +I , a boundary value or a transition range between under-resolved shear flow regimes and resolved ones shall be assessed in more detail to further improve independence of results across a large range of grid spacings. The coefficient C d of the under-resolving drag model might be adjusted after assessing model performance in further test cases. Moreover, a more sophisticated drag modelling approach for under-resolved multiphase flow with smooth or rough interfaces might be adapted to the proposed modelling framework, such as the one proposed by Coste. 38 Furthermore, the information about the estimated degree of under-resolution might be used for other modeling aspects in the context of multiphase flow simulation, for example, turbulence damping in the vicinity of the interface or controlled morphology transitions between continuous and disperse phases. Finally, the application of the presented framework to more complex and industrially relevant test cases might be focus of future endeavors. The overall goal is to achieve consistent results over a large range of spatial resolutions in the context of the hybrid multiphase model MultiMorph in order to flexibly adapt between high precision with fine grids and reasonable results with coarse grids, depending on available computational resources.

ACKNOWLEDGMENTS
This work was supported by the Helmholtz European Partnering Program in the project Crossing borders and scales (Crossing). The authors gratefully acknowledge the financial support provided by the Slovenian Research Agency through the grant P2-0026. Open Access funding enabled and organized by Projekt DEAL.

CONFLICT OF INTEREST STATEMENT
The authors declare that they have no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in RODARE at https://rodare.hzdr.de/record/1877, reference number 1877.