Toward a fully resolved volume of fluid simulation of the phase inversion problem

This paper presents an enstrophy-resolved simulation of the phase inversion problem using the volume of fluid (VOF) method. This well-known benchmark for modeling multiphase flows features a buoyancy-driven unsteady motion of a light fluid into a heavy one followed by several large- and small-scale interfacial processes such as deformation, ligament formation, interface breakup, and coalescence. A fully resolved description of such flow is advantageous for a priori and a posteriori evaluations when developing new subgrid-scale closure models for large eddy simulation of two-phase flows. However, most of the previous attempts in performing the direct numerical simulation of this problem have been unsuccessful to reach grid-independent high-order flow statistics such as enstrophy. The key contribution of this paper lies in proposing a new converging configuration for this problem by reducing the Reynolds and Weber numbers. The new setup reaches grid convergence for all the flow characteristics on a 5123\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$512^3$$\end{document} grid. Particularly, the enstrophy which has always revealed a grid-dependent behavior in all the previous studies converges for the proposed setup. Also, we analyze the temporal evolution of interfacial structures including the statistics of the total interfacial area during the process on different grid resolutions. First, no convergence on the interfacial area is observed and the possible reasons for lack of convergence are discussed. The potential remedies are investigated through a comprehensive parameter study. The findings highlight that (i) the enstrophy always converges for these moderate Re and We numbers, and (ii) the convergence of the total interfacial area is sensitive to the choice of initial and wall boundary conditions. Then, a new setup based on this sensitivity analysis is proposed that succeeded in full convergence for enstrophy and a partial convergence for the total interfacial area. The numerical simulations were carried out using the VOF solvers of OpenFOAM with a comparison between the algebraic and geometric schemes. Besides, the convergence of size distribution of dispersed structures is investigated. The present study provides insight into the possible directions toward a DNS of phase inversion problem with all the flow and interfacial structures resolved, which is essential for the future development of multiphase flow models.


Introduction
Modeling and simulation of multiscale two-phase flows in the presence of interfaces has been a prominent research topic in fluid mechanics and heat transfer. A wide spectrum of real-life applications from controlling the fuel injection, engine pollution, and light material design in the aeronautical industry to various engineering processes such as oil transport, steel manufacturing, and nuclear safety benefit from model development for twophase interfacial flows, in particular when the turbulence is involved. In all these applications, the physical core phenomenon consists of multiphase flow with interfaces between immiscible fluids which in turn are subject to unsteady or turbulent interactions and impart a wide range of spatial and temporal scales from micron to meter and micro-second to minute. This multiscale nature entails a crucial challenge to numerical simulation. The continuous growth of computational capacities brought by massively parallel clusters has nowadays allowed to significantly improve the understanding of highly unsteady and turbulent multiphase flows by performing direct numerical simulation (DNS). In the case of interfacial flows, all the flow characteristics, as well as all the interfacial scales, are intended to be resolved using the one-fluid formulation where the Navier-Stokes equations are coupled with an interface tracking/capturing method. Examples include atomization of liquid jets [1][2][3][4], droplet-laden turbulent flows [5,6], bubbly flows [7,8], nano-coating by plasma projection [9,10], and submerged gas injection in metallurgical flows [11]. Nevertheless, despite the use of extremely high grid resolution and very small time steps as well as implementation of sophisticated numerical methods such as Adaptive Mesh Refinement (AMR) techniques [12,13], coupling the volume of fluid (VOF) [14] and Level Set (LS) [15] methods for interface tracking between fluids [2,16] and employing multi-grid parallel solvers [13,17]), the perspective of obtaining a real DNS where all the flow scales and interfacial features are resolved is still imprecise. In other words, extracting unique and realistic interfacial characteristics such as droplet size distribution and representative interfacial areas in most of turbulent interfacial flows is still impractical when the ratio between large and small resolvable scales exceeds 10 4 . Indeed, in this flow configuration at least a grid resolution of 10 12 is required for resolving all the scales in a three-dimensional domain. Whereas currently, the largest simulations on massively parallel computers cannot handle grids larger than 10 10 cells [18]. In addition to the memory size difficulties, the time steps have to become very small (smaller than 10 −8 s) for numerical or physical reasons. This entails another limitation in particular for those physical phenomena that have to be analyzed for long time intervals. Even with high spatial and temporal resolutions, the simulation of turbulent interfacial flows may still be affected by unrealistic results such as numerical coalescence as well as non-converged droplet size distributions to a unique solution. Therefore, the multiscale nature of such turbulent multiphase flows that involve a spectrum of physical length and time scales would make the DNS still a challenge [7,19], and the capability of the current numerical tools to perform DNS for realistic problems involving multiscale interfacial flows is disputable. Nevertheless, in spite of the limitations, direct numerical simulations are still necessary as the reference solutions to assess the accuracy and performance of new model developments for other affordable approaches such as large eddy simulations (LES). The lower level of maturity of the two-phase LES compared to single-phase one, and the continuous research on modeling new subgridscale (SGS) models in this context make this necessity of proper reference problems even more compelling. Thus, such reference solutions have to be defined with simple initial and forcing conditions in a way that a DNS becomes practical on a reasonable grid resolution (less than 10 9 grid points). Also, they have to exhibit unsteady (or turbulent) behavior for both flow and two-phase characteristics, including deformation, rupture, and coalescence of interfacial structures and represent the interactions between flow turbulence and multiphase physics. One such problem is the interfacial homogeneous isotropic turbulence (HIT) which consists of the evolution of an initially flat interface in interaction with an isotropic turbulent flow. This problem has been investigated to acquire a more in-depth understanding of turbulence-interface interaction [20,21]. Nevertheless, this benchmark was restricted to capillary forces in interfacial flow with no density and viscosity contrasts. Thus, it does not lie within the scope of the present study.
A more popular benchmark in this context is the phase inversion problem. This problem features a buoyancydriven unsteady interfacial flow which begins with an upward motion of a light fluid which is initially placed at rest while surrounded by a heavier fluid. This unsteady motion is then followed by several large and smallscale interfacial processes such as deformation, ligament formation, interface rupture, and coalescence. The interesting feature of this problem is that the two-phase flow unsteadiness, as well as the interface evolution toward fragmentation and coalescence, can be controlled by the corresponding dimensionless numbers, e.g., Reynolds (Re = ρU g L/μ) and Weber (We = ρU 2 g L/σ ) numbers. These numbers are usually defined according to the gravity-based characteristic velocity U g = ρ 2 −ρ 1 ρ 1 gL 2 [18], where the subscripts 1 and 2 denote the lighter and heavier fluid, respectively. During the initial stages of this problem, the inertia is generated by buoyancy forces and becomes dominating, whereas during the latest stages of the motion, inertia decreases similarly as in a homogeneous turbulent flow (even though the flow is not isotropic). We have been working on the phase inversion benchmark for the past two decades and so far several 2D and 3D configurations for this problem have been reported [18,[22][23][24][25][26][27]. Also, this problem has been subject to many a priori and a posteriori LES studies [27][28][29][30]. Nevertheless, so far no real DNS solution for this problem has been reported in which all the flow and interfacial quantities reach the grid convergence. In almost all previous studies, the macroscopic flow quantities such as potential and kinetic energies reveal grid independence on very fine grids, but the high-order moments such as enstrophy as well as the evolution of the interfacial area have never achieved a unique grid-independent solution. This may seem as a generally accepted shortcoming by the community; however, Sayadi et al. [31] have recently investigated the plausibility of one-fluid formulation to achieve gridconvergence for this problem. For this purpose, they reduced the Re and We numbers from O(10 3 ) in previous works [24,26] to O(10 2 ) and performed VOF simulations using the AMR technique. Even for such moderate Re and We numbers, their simulation results underline the loss of convergence for enstrophy which might be attributed to the sheet breakup at the top of the box. Nevertheless, we believe the shortcomings in the previous works originate from a couple of improper physical and numerical assumptions. Indeed, in all previous DNS or highly resolved simulations of phase inversion slip walls were considered which seem problematic as no vorticity length scale is imposed by this kind of wall condition. In addition, employing the AMR technique (as used in [31]) for the capturing of thin film or ligament structures near the walls might be also a reason for the unconverged enstrophy.
The major objective of this study is to revisit the abovementioned limitations and modify the phase inversion problem to allow for a fully resolved solution and try to resolve all the flow scales. Thus, we propose a new configuration in which the Re and We numbers are reduced to O(10 2 ) to become comparable to [31]. We consider the no-slip condition for all the surrounding walls, as it is more realistic due to imposing a vorticity scale following a Re 1/2 law. The numerical simulations were carried out for various static grids starting from 64 3 to 512 3 where the grid convergence is achieved for all the flow quantities including enstrophy. Furthermore, the capability of the VOF approach in providing grid convergence for the interfacial area is investigated, and a lack of convergence is observed. Then, a comprehensive parameter study is carried out to investigate the sensitivity of enstrophy and interfacial area convergence to the various parameters such as the numerical scheme, the initial position of the lighter fluid, and the type of wall boundary condition. First, a comparison between two different VOF methods, namely algebraic and geometric, reveals that the proposed setup is enstrophy-converged independent of the numerical scheme. Then, the lighter fluid is initially placed at a certain distance from the surrounding walls in order to facilitate the free dynamics of the lighter fluid. While the enstrophy is still converged, it prevents the formation of sticky thin films near the walls at the early stage of the process, resolving which requires higher grid resolutions and may induce convergence issues for interfacial structures. Finally, the slip boundary conditions are tested at the surrounding walls. The results indicate that the slip boundary condition results in a different flow topology, in particular, by preventing the formation of a thin layer of heavier fluid near the upper wall. This thin layer is identified as the source of the non-convergent interfacial area. According to the findings of this parameter study, a new setup with a combination of these initial and boundary conditions is proposed which not only converges on the enstrophy but also succeeds in a partial convergence on the interfacial area.
The paper is structured as follows. In the next section, we present the governing equations as well as the numerical details for the volume of fluid simulations. In Sect. 3, we discuss the results of the enstrophy-resolved simulations and present a detailed parameter study. Then, the paper ends with conclusions in Sect. 4.

Governing equations
The unsteady motion of an incompressible, immiscible, multiphase flow can be mathematically described with the continuity and Navier-Stokes equations together with the transport equation of the marker function of the interface capturing technique, i.e., here the volume of fluid method. These equations are defined in the frame of one-fluid formulation as ∂α ∂t + ∇ · (αU) = 0. (3) In this system of equations, U is the mixture velocity vector shared with all phases. p is the pressure and D is the strain rate tensor in the form of D = 1 2 [(∇U) + (∇U) T ]. The scalar function α is the volume fraction field which is used to determine the mixture density and viscosity of the fluid as ρ = αρ 1 + (1 − α)ρ 2 and μ = αμ 1 + (1 − α)μ 2 , respectively. The surface tension force, S σ , is commonly computed by the Continuous Surface Force (CSF) method [32]. In this approach, a basic definition of interface normal vector is applied based on the gradient of VOF function, then the unit interface normal, and its curvature are determined aŝ n = ∇α |∇α| and κ = −∇ ·n, respectively. Using these two quantities, the surface tension force reads where σ is the surface tension coefficient and δ s ≡ |∇α| is the mathematical delta function that equals infinity at the interface and zero elsewhere. Therefore, the accuracy in calculation of surface tension is directly influenced by the accuracy of VOF approach in predicting the interface characteristics.

Volume of fluid method
The choice of interface capturing technique is significantly important in modeling interfacial flows. The estimation of the interface normal vector and curvature out of discretized marker function of α in VOF is challenging and may induce numerical errors in the regions involved with interface and flow interactions. All the interface capturing techniques including VOF consist of two major steps that have to be fulfilled accurately: advection of the interface marker function, and estimation of interface curvature. There are two families of advection algorithms, namely, geometric and algebraic methods. Geometric methods are based on reconstruction-advection of the interface. In these methods, first, the interface is explicitly reconstructed inside the computational cell, then advected in the Lagrangian sense to accurately compute volume fluxes and reconstruct the sharp interface. A well-known example of this approach is the Piecewise Linear Interface Calculation (PLIC) scheme which is used by many computational codes. In contrast, algebraic methods do not involve complex interface reconstruction procedure, and hence computational numerical sharpening treatments are required to prevent interface smearing. One such effort involves adding an artificial compression term to Eq. (3) to make the interface sharp. This is the essence of the Multidimensional Universal Limiter of Explicit Solution (MULES) algorithm implemented in the VOF solver package of OpenFOAM [33] called interFoam. This solver has been used in our previous studies for different industrial applications of interfacial flows [34,35]. For more details on the computational algorithms and parameters, we refer to [36,37], where this algorithm is also tested for pertinent benchmarks such as splashing droplet. As an efficient alternative to geometric VOF, OpenFOAM has been known for its algebraic VOF scheme for a decade until the contribution of Roenby et al. [38] who introduced an OpenFOAM-based geometric VOF called isoAdvector. They replaced the MULES algorithm with a two-step reconstruction-propagation operation based on the iso-surface of the volume fraction in each computational cell. First, during the reconstruction step, they use the iso-surface calculation to find the interface position and orientation in each interfacial cells. Then, the intersections of the reconstructed interface with the cell faces are further considered to obtain the time evolution of the interface within a time step. This approach provides the accuracy of geometric methods while the computational complexity is kept at a minimum. They demonstrate satisfactory results for 2D and 3D benchmark problems, where the algebraic VOF performs poorly. The details of this geometric approach can be found in their paper [38] and thus are not repeated. In the present study, the geometric VOF of isoAdvector is employed for all the simulations and analysis, and the algebraic VOF is only used for one of the items in the parameter study.

Simulation setup
According to the schematic description of the problem in Fig considered to account for the vorticity length at the walls. The rest of the domain is filled with the heavier fluid (hereinafter fluid-2). The cavity is fully closed with no-slip walls. Both phases are assumed initially at rest while the buoyancy-driven motion starts in the y-direction immediately after the process begins. The fluid-1 and fluid-2 densities are ρ 1 = 900 kg/m 3 and ρ 2 =1000 kg/m 3 , respectively. An equal dynamic viscosity is assumed for both fluids with the value of μ 1 = μ 2 = 1 Pa s. The interfacial tension coefficient of σ = 0.45 N/m is considered between the phases.
According to [18], a gravity-based characteristic velocity of U g = ρ 2 −ρ 1 ρ 1 gL 2 = 0.245 m/s can be determined. Based on the fluid-1 material properties and the characteristic velocity, the Reynolds and Weber numbers are 221 and 121, respectively. In the context of the finite volume method, the Navier-Stokes and VOF transport equations are solved numerically. A combination of PISO and SIMPLE algorithms, namely PIMPLE [33], is employed for the pressure-velocity coupling and pressure correction. A first-order implicit discretization scheme is used for the temporal term. The convective terms are discretized using a second-order accurate interpolation scheme (central differencing). The gradient terms are discretized by using a Gaussian second-order linear scheme. For the discretization of the convective term in the VOF equation, a Van Leer scheme [39] is utilized for both algebraic and geometric approaches. The computational domain is discretized uniformly with the size of L/2 n in each direction, where n ranges from 6 to 9 to provide four various cases with 64 3 , 128 3 , 256 3 and 512 3 number of computational cells. Having applied both geometric and algebraic VOF schemes for each case, in total 8 cases were simulated. Based on a time step dependency analysis, a simulation time step of Δt = 1e−4 s was chosen for all the cases. An in-house computational cluster was used for the simulation with 64 processors for the lowest grid resolution cases and 512 processors for the highest grid resolution case. The maximum computation time was three weeks for the highest grid resolution of 512 3 until the t = 10 s of the real time in phase inversion process. It has to be noted that the kinetic energy of both phases goes down to zero before this instant, therefore we simulated the highest grid resolution cases until this time, and the corresponding flow sequences will be presented.

Results and discussion
The numerical simulations using the volume of fluid method are performed and different results are obtained. For comparison purposes, the quantities obtained from the simulations are all made dimensionless. Similar to [26], a dimensionless time t * = t t c is defined based on the characteristic time scale for the fluid-1 to move alongside the domain which reads t c = L U g . The macroscopic evolutionary behavior of the liquid-liquid interface for the proposed configuration is presented in Fig. 2. The unsteady motion initiates due to the buoyancy forces which move the fluid-1 upward. Fluid-1 remains in one major piece until around t * = 1. At this instant, the interface exhibits several stretching followed by rippling from different sides, and generates several thin liquid layers prone to breakup. At t * = 1.5, those thin films disintegrate into many ligaments that further break up into small droplets at t * = 2.5. The process continues with several small-scale interfacial sub-processes such as fragmentation of droplets into smaller ones that move upward until the major part of the fluid-1 reaches the top of the cavity. At t * = 3.7, many droplets are still dispersed in the domain, but the process approaches its equilibrium state where the kinetic energy level is almost zero and the dispersed droplets are gradually moving upward. If the simulations further continue, the fluid-1 will eventually reach equilibrium position at the upper L/8 of the cavity. It has to be noted that the flow sequences in Fig. 2 are obtained from the 256 3 grid as this case was simulated until t * = 5. Also, unless otherwise stated, all the flow visualizations in the remainder of this study are based on the results of the geometric VOF scheme. We will further discuss all the simulation cases from two different aspects: flow characteristics and interfacial structures.

Flow characteristics analysis
For the direct numerical simulation, the global physics of the flow is fully resolved if all the flow characteristics of the flow converge to a unique physical solution independent of the computational grid resolution. To evaluate the convergence of low-order moments, the temporal evolution of the kinetic and potential energies in both fluid phases is considered. Figure 3 presents the temporal evolution of the domain-integrated kinetic energy (KE = Ω 1 2 α i ρU 2 dV ) and potential energy (PE = Ω α i ρ g ydV ) of fluid-i on different grid resolutions and VOF schemes. The kinetic energies in fluid-1 and fluid-2 are made dimensionless by the initial value corresponding to the characteristic velocity as respectively. The potential energies in fluid-1 and fluid-2 are normalized by the potential energy of each phase at the final equilibrium condition (i.e., t * → ∞) as P E 1 = 15 128 ρ 1 gL 4 and P E 2 = 49 128 ρ 2 gL 4 , respectively [26]. Immediately after the process begins, there is a drastic increase in the KE in both phases. Fluid-1 reaches its maximum KE at t * = 0.4, while this takes a bit longer for the fluid-2. Then, the KE begins to decrease at both phases until it reaches zero around t * = 2 and t * = 3 for fluid-1 and fluid-2, respectively.
The total potential energy has a clear evolution during phase inversion at both phases. As fluid-1 is moving upward, its PE increases continuously until t * = 1.75. At this instant, the major content of fluid-1 has reached the equilibrium position and the upward motion of the dispersed liquid structures increases the PE with a much slower rate until around t * = 2.25 where it becomes nearly constant for the rest of the process. Conversely, fluid-2 loses potential energy until t * = 1.75 where its major content reaches the equilibrium. From this moment, the PE still decreases but with a very slight slope until t * = 2.25 where it reaches its equilibrium value. The grid sensitivity analysis reveals that both energies converge at the grid resolution of L/256 and higher.
As the major objective of this study, the convergence of second-order statistics of the flow is examined. Figure 4 presents the temporal evolution of the domain-integrated enstrophy (E = Ω 1 2 α i (∇ × U) 2 dV ) for each phase. The enstrophy in each phase is made dimensionless based on the maximum enstrophy resolved by the highest grid resolution. At the early stage, the enstrophy increases drastically as the buoyancy force derives the fluid-1 to move upward. This abrupt motion entraps a portion of fluid-2 between the fluid-1 and top wall which generates strong vortices in the upper region. This is reflected as the peak of enstrophy in fluid-1 at around t * = 0.5. Afterward, the maximum magnitudes in fluid-2 happen around t * = 0.85, which corresponds to the highest shear rate of the two-phase flow and the interaction between the vortices in both fluids. As time elapses, the enstrophy decreases gradually due to the damping of momentum under viscous force action. Both phases reach their zero levels for t * > 2.25. It is evident from Fig. 4 that the enstrophy of both phases has reached the grid convergence.
The convergence of enstrophy is a key finding for further physical investigations. The smallest flow structures are proportional to the Kolmogorov microscale which can be computed by the kinematic viscosity and turbulence dissipation rate as The dissipation can be approximated by the maximum enstrophy as ≈ 2νE max [23]. Using the maximum domain-integrated total enstrophy for the highest grid resolution of 512 3 , the smallest flow structure can be  computed concerning the material properties of each phase as presented in Table 1. It is evident that for the computation based on either of the liquids, the grid spacing is sufficiently smaller than the Kolmogorov length scale. It has to be noted that at such a moderate Re number we do not expect the flow to be turbulent. Therefore, we computed the Kolmogorov length scale purely based on the definition to indicate the size of the smallest flow structures and the adequacy of the grid size in this study.

Interfacial structures analysis
The evolution of fluid-fluid interfaces is of great importance in real-life problems involving interfacial energy, chemical reactions and mass transfer. A macroscopic characteristic of interfacial flows is the total interfacial area during the process. The interfacial area can be approximated using the mathematical proportionality to the gradient of volume of fluid (i.e., Σ ∝ ∇α) [27]. Therefore, the total interfacial area in the domain is computed by where S cell = (Δx) 2 and Ω is the computational domain. It should be noted that Eq. (6) is simple to implement but not the most accurate way to compute the interfacial area. Vincent et al. [18] show that this approximation could result in a maximum 20 % difference in interfacial area compared to the proper integration of the interface in each computational cell. Nevertheless, it is still plausible to study the convergence of Σ for the phase inversion problem. Initially, the fluid-1 is at rest in the corner of the cavity in touch with three walls. Thus, the total interfacial area equals 3(L/2) 2 , i.e., 0.75 m 2 . This value is used to make Σ dimensionless. The temporal evolution of the domain-integrated Σ for different grid resolutions is presented in Fig. 5. After a short while that it almost keeps its initial magnitude, the upward motion of fluid-1 causes drastic stretching and large deformation of the interface that yields a rapid increase in the interfacial area. Although fluid-1 remains in a major piece together with few ligaments, the interface reaches a maximum total area at t * = 1.1. As 3 (a) grid 256 (b) grid 512 3

Fig. 6
The iso-surfaces of liquid-liquid interface at the instant of t * = 1 for 256 3 and 512 3 grids. The pointed under-resolved zones A and B are symmetrical with respect to diagonal plane time elapses, the corrugated interface is fragmented into many ligaments, which further break up into small droplets. The total interfacial area decreases by time as several droplets move upward and coalesce with the major portion of the fluid-1 accumulated at the top of the domain. It is expected that the interfacial area reaches its equilibrium value during much later stages of the process which is 1 m 2 corresponding to dimensionless area of 1.33.
In contrast to the potential and kinetic energies as well as the enstrophy, the total interfacial area seems not to converge to a unique physical solution, and higher grid resolution results in resolving more interfacial structures. But as it is evident in Fig. 5, the Σ seems to almost converge until around t * = 0.7 at the highest grid resolution. Shortly after this instant, a drastic change between the total interfacial area on the 256 3 and 512 3 grids arises that remains nearly the same for the rest of the process. For a better understanding of the source of this change, the iso-surfaces of the liquid-liquid interface at t * = 1 are shown from the side view for the 256 3 and 512 3 grids in Fig. 6.
The upper left regions look almost alike but the thin rippling sheets on the top-right regions (zone A) as well as the vertical sheet (zone B) reveal the major difference between the cases. Besides, a couple of symmetrical slim ligaments are resolved by the 512 3 grid. The thickness of these sheets and ligaments is in the size of the grid resolution. Therefore, they remain under-resolved for the 256 3 grid and result in the lack of convergence for the total interfacial area. It is also observed in Fig. 5 that the difference between the grids in predicting the maximum interfacial area reduces by increasing the grid resolution. Thus, it can be expected that the total interfacial area may converge in the case of higher grid resolutions such as a 1024 3 grid or higher. However, such simulation was beyond the computational resources in the present study.

Parameter study
The present study on the phase inversion problem is distinguished from the previous studies by reducing Re and We numbers and applying no-slip boundary conditions for the velocity at the walls. As a result, the setup succeeded in a full convergence for enstrophy on the 512 3 grid. To evaluate the generality of this conclusion, a parameter study is carried out. As a comparison measure, we only focus on the enstrophy in fluid-2 as well as the interfacial area Σ. For better comparison in the remainder of the paper, the enstrophy and interfacial area are made dimensionless by the maximum resolved enstrophy and initial magnitude of area in the original case, respectively.

Effect of VOF scheme
First, the proposed configuration was simulated by the algebraic VOF method. The temporal evolution of the domain-integrated enstrophy of fluid-2 obtained on different grids is presented in Fig. 7. No specific difference is observed among the VOF schemes indicating that even algebraic VOF is capable of reaching the enstrophy convergence for the proposed configuration of the phase inversion problem. It has to be noted that the results on the 256 3 grid are already conclusive, thus the simulation on the 512 3 grid is skipped to minimize the computational costs.
The effect of VOF schemes on the integral characteristics of interfacial structures is also explored. Figure  8 compares the domain-integrated interfacial area obtained by each VOF method. It is evident that both VOF schemes are capable of predicting the evolution of Σ with the same level of accuracy. But most of the time the geometric VOF yields a higher total interfacial area compared to its algebraic counterpart. This holds true except for t * > 1.7 on the coarsest grid of 64 3 where the trend is vice versa. This is attributed to the fact that the algebraic VOF smears the interface over a few cells. Therefore, the |∇α| is nonzero in more cells compared to the one-cell-occupying, sharp interface computed by the geometric VOF. This leads to an overestimation of Σ in Eq. (6).

Effect of initial condition
As the next parameter, we study the effect of initial conditions on the convergence of enstrophy and interfacial area. We repeated the simulations for the same Re and We numbers with a new initial condition where the lighter fluid is initially placed at a certain distance h from the surrounding walls. This configuration is depicted in Fig. 9. In this parameter study, h = L/10. Figure 10 presents the fluid-2 enstrophy on different grid resolutions for both initial conditions. It is evident that the enstrophy is still converging; however, a different trend is observed which originates in the different dynamics of the fluid-1 during the process. When fluid-1 is initially attached to the walls (original setup), its upward buoyancy-driven motion is affected by the formation of a thin layer at the wall which results in a lower shear rate in fluid-2. But when the initial position is shifted from the walls, the free dynamics of the lighter fluid is facilitated. Therefore, the enstrophy reaches its peak faster with a slightly higher magnitude. This trend is similar for all the grids and reveals the convergence. Although the initial position of fluid-1 has no impact on the enstrophy convergence, we further evaluate the convergence of the interfacial area on both setups as shown in Fig. 11. It shows that by shifting the initial position of fluid-1 far from the walls, it remains more intact during the upward motion with no interface-wall interactions. This strongly reduces the magnitude of the interfacial area on all the grids compared to the original setup. Also, the formation of the elongated thin layer at the walls is prevented, resolving which requires higher grid resolutions and makes the results more grid-dependent. Accordingly, shifting initial position improves the partial convergence of Σ until t * = 0.85 as shown in Fig. 11. After this time, the temporal evolution of Σ reveals a similar trend on different grids with no convergence on the higher grid resolutions. But, the difference in Σ between the higher grid resolutions becomes smaller. Because as shown in Fig. 12 in the absence of the elongated thin film at the walls, the only contributors to the lack of convergence are the unresolved dispersed droplets that may be eventually resolved on much higher grid resolutions. It can be concluded that shifting the initial position of the lighter fluid could be a potential remedy for the lack of convergence in the interfacial quantities of the phase inversion problem.

Effect of the wall boundary condition
It remains to discuss the effect of the wall boundary conditions on the convergence of the flow quantities. As mentioned in the previous sections, the converged enstrophy values presented in Fig. 4 were obtained by imposing no-slip conditions for the velocity at the surrounding walls. We repeated these simulations by changing the no-slip boundary condition to the slip one. The fluid-2 enstrophy is plotted in Fig. 13 for both boundary conditions. Interestingly, the enstrophy still remains converged despite yielding different magnitudes. Similar to the previous cases, the enstrophy results obtained on the 256 3 grid are already conclusive, thus the simulation on the 512 3 grid is skipped to minimize the computational costs. It is evident that employing slip BC at the walls results in a much lower enstrophy in fluid-2 compared to the no-slip counterparts. Because the no-slip condition usually imposes a velocity constraint that enhances the vorticity generation at the wall regions which eventually results in higher enstrophy. This is demonstrated in Fig. 14, where the vorticity magnitudes at the middle-plane in the x-direction are illustrated for the instant of maximum enstrophy in the original case (i.e., t * = 0.85). The high vorticity regions at the top and the left walls clearly indicate the higher vorticity generation due to the no-slip constraint at the walls. Moreover, a shift in the instant of maximum enstrophy is observed in Fig. 14 which is attributed to the different procedures of vorticity generation in each case. The major source of vorticity generation in the case of no-slip wall BC is the shear rate near the walls which happens immediately after the fluid-1 starts moving upward and reaches the upper wall. But when the slip BC is imposed at the walls, the fluid-1 reaches the upper wall, slides over it, and then slides down from the other walls to the bottom of the cavity. The latter downward motion results in the highest vorticity generation which is of lower magnitude compared to the no-slip case.   Finally, the effect of slip boundary condition on the evolution of interfacial structures is studied. Figure 15 compares the temporal evolution of Σ for both slip and no-slip wall BCs. It is evident that the interfacial area is still not converging, and imposing slip BC has no impact on the maximum Σ on each grid, but causes a shift in the instant of maximum Σ, similar to the enstrophy. In fact, in the case of slip walls, Σ increases with a gentler slope and reaches its peak at many later times (t * = 1.6) compared to the no-slip cases. But after this instant, it decreases with a relatively steeper slope in comparison with the no-slip case. Thus, it is expected that Σ reaches its equilibrium value (i.e., 1.33) sooner than the original setup. This can be explained by the difference in liquid-liquid interface topology at the top of the cavity. As the buoyancy-driven dynamics of both systems are similar, the lighter fluid reaches the upper wall at almost the same time (t * = 0.75) regardless of the type of wall BC. In the case of no-slip BC, there is a high shear rate near the upper wall to satisfy the zero velocity condition. This hinders the motion of both fluids near the wall and subsequently causes the entrapment of a thin film of fluid-2 between the fluid-1 and at the upper wall of the cavity. This film is the main source of non-convergence of Σ (as also observed by [31]), and remains until the end of the process. Figure 16 compares the 3D snapshots of the interfacial structures as well as the volume fraction contours on the upper wall plane for both no-slip and slip cases at t * = 2 which is long after the upward motion of fluid-1. Even though fluid-1 is accumulated at the top of the cavity regardless of the wall BC, a thin layer of fluid-2 is trapped near the upper wall in the case of the no-slip case (i.e., the regions with zero volume fraction in Fig. 16a). But the slip wall prevents from formation of such film at the top of the cavity as both fluids slide freely over it. Most of the time, the upper wall is covered by the fluid-1 as shown in Fig. 16b. In this case, the only source of non-convergence of Σ is the dispersed structures such as ligaments and droplets that are mostly formed at t * > 2 and may remain under-resolved on the coarser grids. In Fig. 15, the absence of this thin film in the total interfacial area is reflected as the difference between the slopes of Σ for slip and no-slip cases at t * > 2. It is concluded that while the enstrophy convergence is not affected by the type of velocity boundary conditions at the walls, imposing slip walls could reduce the source of non-convergence for the interfacial area by preventing the thin fluid-2 layer at the top of the cavity.

A new setup for the convergence of the total interfacial area
To recapitulate, the parameter study shows that for such moderate Re and We numbers the enstrophy is always converged on the 512 3 grid no matter what type of VOF scheme, initial and boundary conditions are used. However, the convergence of the total interfacial area reveals a considerable sensitivity to some interfacial topologies that, in turn, are triggered by the initial condition as well as the type of wall boundary condition. Particularly, the results show that shifting the initial position of the lighter fluid away from the walls facilitates the upward motion of lighter fluid and prevents the formation of a sticky film at the wall. Also, imposing slip walls could avert the entrapment of a thin layer of the heavier fluid near the upper wall. Thus, one could control the convergence of Σ by proper changing of these parameters. Following this conclusion, we propose and test a new setup that benefits from a combination of both parameters. Figure 17 depicts the new configuration in which the fluid-1 is initially placed away from the walls with a slip wall boundary condition imposed only on the upper wall plane. The rest of the walls are still no-slip as this is a more realistic boundary condition for many real-scale applications. The domain-integrated enstrophy of fluid-1 is plotted in Fig. 18. It is clear that the new setup is still enstrophy-converged. Also, in Fig. 19 the total interfacial area of the new setup is compared with the one obtained from the original case (with no-slip walls and attached initial position). In the new setup, the overall magnitudes of Σ have decreased significantly as a result of shifting the initial position. Besides, the results show a partial convergence of Σ for t * > 2.25 on the highest grid resolution. It means that the integral of all the interfacial surface areas is converged to a unique value. This is attributed to the fact that after this time the upper wall is fully covered by the fluid-1 and the large-scale interfacial structures remain almost unchanged. Figure 20 presents a 3D snapshot of the interfacial flow obtained on the highest grid resolution of 512 3 at t * = 2.5 where both enstrophy and interfacial area have converged to a unique physical solution. It is expected that by further increasing the grid resolutions, Σ convergence for 0.7 < t * < 2.25 may also be achievable. Nevertheless, this is beyond the computational capacity of the present study and remains for future studies.
Further to this integral viewpoint, it is also worth investigating the size of individual interfacial structures such as droplets and ligaments to examine the convergence of droplet size distributions. It has to be noted that with the VOF approach, the size of the smallest resolved droplets is controlled by the size of computational grid cells, making the convergence of droplet size distribution a difficult task. To analyze the size distribution of droplets, a post-processing tool based on the connected component labeling concept is employed to detect individual dispersed droplets in a three-dimensional snapshot of the flow sequence. The size of droplets is then measured with an equivalent diameter. The Probability Distribution Function (PDF) is defined as the number frequency of droplets. The equivalent diameters d i may belong to any of a number of adjacent bins B i with the width of Δ i (i.e., B i = (d i − Δ i , d i + Δ i )). Further details on this post-processing algorithm are provided in [18]. Figure 21 presents the PDF of droplets of both the original and the new setups at t * = 2.5 where Σ is converged for the latter. In the plots, the bin width varies with grid resolution and is twice the corresponding mesh size for each case [18]. It is evident that the PDF of the droplets in neither of the cases is converged by the grid refinement. In fact, by increasing the grid resolution much smaller droplets are being resolved. This is consistent with the major shortcoming of the interface capturing techniques. The results further imply that the convergence on the integral quantity of Σ does not necessarily constitute the convergence of the size distribution to a unique PDF, because different droplet populations could have the same domain-integrated surface area. This is the case for the PDF of the droplets obtained from the 256 3 and 512 3 grids in Fig. 21b, where the total interfacial area is converged but the number of smaller droplets increases by grid refinement. In addition, we plotted the distribution of area per class of the droplets as the inset plots in Fig. 21. It shows the cumulative area of the droplets in each class of diameter d i normalized by the initial total area in the original case, i.e., A = 1 A 0 n i πd 2 i , where n i is the number of droplets in each class. The integral impact of different ranges of droplet sizes is better revealed in the area per class distribution plots. The smaller classes, despite their large population, have a rather negligible contribution to the total interfacial area, whereas most of the interfacial area in the domain comes with the larger structures. It is also noticeable in Fig. 21 that the cumulative area of most of the droplets with d i > 0.02 m is almost the same for both the 256 3 and 512 3 grids which are reflected in the convergence of Σ for the new case. But such a trend is not evident for the original setup. Furthermore, the converging trend in the area per class plots implies the necessity to define more pertinent measures for the statistical representation of two-phase flow features in the context of interface capturing simulations, since the convergence of the droplet numbers is sensitive to the grid resolution and may not be achieved during the VOF simulations. It has to be noted that the present post-processing tool does not distinguish between different shapes of structures and assigns the size of an equivalent spherical droplet to each detected dispersed structure. Therefore, another feasible measure for the statistical representation of two-phase flows could be the size distribution of different types of interfacial structures and convergence analysis for the PDF of individual shapes such as ligaments of a specific size. However, developing a post-processing tool for such detailed analysis is the subject of ongoing research and is planned for future studies.
To conclude, the new setup with a full convergence on enstrophy and a partial one on the total interfacial area could be considered as the basis for future modifications of the phase inversion problem toward a fully resolved simulation where the full convergence on all integral quantities is expected. It is also essential to

Conclusions
In this study, the VOF simulation of the phase inversion problem is presented for a configuration that provides the grid convergence in different first-and second-order statistics of the flow. Particularly, the enstrophy that has always revealed a grid-dependent behavior in all the previous studies is convergent on a 512 3 structured grid. The numerical simulations were carried out on uniform static grids using a geometric VOF solver in the framework of OpenFOAM. It has to be emphasized that while the previous works could not achieve grid-independent results, this study provides a three-dimensional fully resolved flow where all the energies and enstrophy converge. We also performed a comprehensive parameter study to evaluate the generality of our conclusion. It is concluded that regardless of the type of VOF method as well as the type of initial and boundary conditions, the enstrophy for such moderate Re and We numbers is convergent on the highest grid resolution. We also analyze the temporal evolution of the total interfacial area on different grid resolutions.
Although the total interfacial area shows no convergence first, the potential remedies are investigated through a parameter study. The findings highlight the effect of the liquid-liquid surface topology near the walls (i.e., formation of different film layers of either fluid) on the convergence of the interfacial area. Then, a new setup based on this sensitivity analysis is proposed that succeeded in full convergence for enstrophy and a partial convergence for the interfacial area. The latter shows the possibility of reaching a full-convergence on the grid resolutions higher than 512 3 , which remains for future studies. We also believe this setup can serve as the reference solution for future LES studies such as a posteriori evaluation of the recent functional SGS models [30] as well as the structural ones [26]. In particular, the potential implication of the present results is to evaluate the performance of closure models for the subgrid surface tension and subgrid term in the filtered VOF equation that are essential in capturing subgrid interfacial structures. Future works will focus on the evolution of individual interfacial structures and analyzing the convergence of their size distribution. Also, we will use the herein proposed setup to explore higher Re and We numbers.