Evaluation of the interFoam solver in the prediction of immiscible two-phase flow in imbibition and drainage on the pore-doublet system

. Drainage and imbibition at low to moderate capillary number, are very important phenomena, in which the SCAL data are derived. In this venue, there is an increasing trend to apply computational fluid dynamic models at pore scale to extract the required parameters for multiphase flow study. Before applying these methods, their abilities to capture the real performance of fluids at the pore scale with corresponding capillary numbers should be validated. Moreover, due to the dimensional difference from pore to core scale, rigorous statistical methods are required to successfully extract the corresponding SCAL parameters . Hence for any model, before constructing any SCAL data via digital methods, validation at the pore scale along with the investigation of the possible effect of the dimensional change from pore to core scale are significant priorities. In this study, the performance of interFoam, a solver for immiscible two-phase flow from OpenFOAM, at the pore scale is evaluated. An image of the pore doublet system is digitized and meshed with different resolutions to capture the contact angle movement. The average pore diameters are different, and the pore diameter pathway varies. The experimental results are extracted from an article, which extensively performed the tests for different imbibition and drainage regime on this pore doublet system. Moreover, a sensitivity analysis on mesh geometry, refinement near boundaries, and contact angle is conducted. The drainage simulation and experimental results in the pore doublet system are in good agreement with each other. As expected, the displacement happens in the larger diameter and bypasses the wetting fluid in the narrow constriction. In imbibition, the displacement first takes place in the smaller constriction and then the wetting phase starts invading the larger diameter pore. From the sensitivity analysis, it is concluded that the instability of the front is mainly related to the slip condition at the walls. Moreover, the spreading of wetting film near the wall is prevalent at very low contact angles. The results show that the solver can predict the flow behavior in imbibition and drainage on the pore doublet system and the results are comparable with the experiments.


Introduction
Understanding two-phase flow behavior especially immiscible flow is crucial in different areas of science and industry including the petroleum industry [1].Multiphase fluid flow through porous media, at the macro scale (core and/or field), is often modeled based on the derivation of relative permeability curves which are extracted from drainage and imbibition experiments.Another approach is to extract the pore structure and apply CFD equations at the pore scale to capture the flow behavior at micro-scale.Typically, in hydrocarbon recovery the scale of fluid velocity is on order of cm/s (ft/day) and the flow regime away from the wellbore and even near the wellbore is in the laminar regime.This article aims to investigate the applicability of CFD modeling on pore scale structures and determine if they can reproduce immiscible flow behavior.
OpenFOAM, a leading, free CFD open-source software owned by OpenFOAM Foundation is being used for this purpose [2].OpenFOAM (Open-source Field Operation and Manipulation) is an object-oriented C++ framework capable of performing a variety of computational fluid dynamics in continuum mechanics.The solver uses the finite volume discretization method [3].This software includes solvers for incompressible, compressible flow, and covers multiphase flow modeling.Specifically, the interFoam solver is being used for this study.Deshpande et al. [4] analyzed the two-phase flow behavior at the scale of 10 -6 -10 -3 m under capillary dominant and viscous dominant regimes.In the capillary dominant regime, they investigated the effect of perturbation at the interface of fluids in a standing column.They found interFoam can produce satisfactory results when spurious currents are controlled.Recent modifications in interFoam include an adaptive interface compression scheme to minimize parasitic current [5].However, a study is needed to investigate the pore scale behavior of this solver for drainage and imbibition process.
We use the standard interFoam solver from OpenFOAM version 5.0.In this article, first the formulation of the underlying equations and the solution procedure are presented.Afterward, a meshed pore-doublet is used for the CFD experiment to investigate the capability of the solver to predict flow behavior in imbibition and drainage processes.The sensitivity to the contact angle, and the resolution of mesh near the surface are investigated.

Methodology
The method that is being used for the solver is presented, mostly adapted from [4].In two phase flow, the phase fraction  ,  = 1,2, is the volume of one phase over the total volume of interest: The mass conservation is given as: The density of fluid (which is a mixture of immiscible phases) is calculated as: The result of combining equation 1 and 2 yield the phase fraction equation as: The above equations along with the momentum equation can be used to obtain the phase fraction and the velocity and pressure fields of interest.
The momentum equation is written as: The last two terms on the RHS of the above equation are forces which are related to the gravitational field and interfacial tension between the two phases.These equations can be used at the pore scale for continuous flow.The boundary conditions can be modified to account for the slipflow regime which is prevalent near the surfaces [6].
The interfacial tension is based on the continuous surface force model adapted from [7] and is obtained as: where  is interfacial tension and  is the curvature of the interface defined as: The above equations are coded in interFoam solver with the aid of the finite volume method.The finite volume formulation is presented below (description based on [8]).
Equation 8 is the general balance equation for any physical property being carried in the fluid.
( ) + +.(U) − .  =  () From the LHS of the equation, the first term stands for temporal evolution of property, the second one stands for the convective contribution of property by the fluid and the third one represents the contribution from diffusive distribution of the property.The RHS term represents the contribution of a property from a sink or source within the system.
Integrating over a domain of interest at specified time step yields: Every volume integral can be transformed to an integral over faces using Gauss' theorem: A second order equation discretization in time, known as the Crank-Nicolson method is used: where the F is the mass flux through the faces: This system of equations on every cell will be simplified to: These equations on every cell form a sparse matrix with  on diagonal and ℎ   on the off diagonals: The momentum equation is solved by a combination of SIMPLE and PISO methods.The benefit of the combined procedure is the ability to use a larger time step with the Courant Number higher than unity [9].
The procedure for solving the equation is as follows [10]: 1. Initiate variable based on internal initial (previous time step) and boundary conditions.

Results and Discussion
In Hydrocarbon production practice, the fluid velocity is at very low orders, in a way that capillary forces gain their significant role in immiscible fluid movements.Depending on the relative importance of capillary and viscous forces, the microscopic sweep efficiency may differ [11].Mahmoodi et al. developed a visualization method to capture multiphase fluid displacement in micro-models [12].Saturation evolution were quantified during EOR process at pore scale [13][14][15].The behavior of imbibition and drainage is visualized experimentally based on viscosity ratio and capillary number by Guo and Aryana [16].Based on their findings for the specific microfluidic device used to perform the experiments, stable front movement is guaranteed when log M <0 and log Ca>-4.Here M and Ca are viscosity ratio and capillary number respectively.Raeini et al. [17] used a finite volume method to investigate the controlling mechanisms at pore scale in two-phase flow including the layer flow, snap-off, and the geometry effect on trapping and mobilization of disconnected ganglia.Berg et al. [18] computed the relative permeability curves from pore scale quasi-static method and compared the results with relative permeability obtained from 3-D imaging.
The attempt here is to evaluate the performance of the interFoam solver in the prediction of two-phase immiscible flow at pore scale in drainage and imbibition scheme.A thorough investigation of the problem analytically and experimentally is presented by Chatzis and Dullien [19].In their work, they used the conventional pore-doublet configuration to visualize the microscopic phenomena at narrow constriction on the scale of mm.As displacement may occur in imbibition or drainage schemes, numerical studies on both imbibition and drainage are performed.

Pore Doublet Reconstruction for Numerical Study
To reconstruct the 2-D pore doublet configuration for numerical study, the image of pore-doublet is digitized and imported to Gmsh [20] for meshing.The cells are composed of 54230 unstructured hexahedra grids.The overall domain of the pore doublet structure is 6.94 x 3.19 x 0.5 mm in x, y, z directions.Table 1 shows the overall quality of the mesh where results were obtained by checkMesh, an OpenFoam utility for evaluating mesh quality.A view of meshed structure and representative cell configuration is presented in Figure 1.Right: zoom view of hexahedral grid mesh

Imbibition in a Pore-Doublet
The imbibition process is the displacement of the non-wet phase with the wetting phase due to the contribution of only capillary pressure [21].Generally, for a pore-doublet configuration in either forced or spontaneous imbibition at very low speed when the capillary forces become significant, the wetting phase, first, fill and displace the nonwetting phase at the narrowest pore.Afterward, the process continues to push and displace the non-wet fluid in the larger constriction.The efficiency of the process is reported as 100% and all the non-wet fluid is displaced.
An imbibition experiment performed by Chatzis and Dullien [19] is selected for comparison with numerical results.The imbibition experiment is performed as displacement of n-decane (colored with an oil-soluble dye), by 2% NaCl brine.
The fluid properties for two sets of experiments are extracted from the literature and presented in Table 2. [22,23].The boundary and initial conditions are shown in Table 3.

Numerical setup
In OpenFOAM, the numerical set up to build the discretization matrix for equation 14 is specified in a file called fvScheme.Every term in the differential equation can be discretized based on the specified scheme.The summary of discretization schemes is presented in Table 4.The 2-D solution procedure for the linear system is given in fvSolution file.The method of solving alpha, pressure, and speed are those presented in Table 5.Two iterations for the  loop and one iteration for its sub-cycle are selected.The compression parameter for the saturation equation is set to 1 and the PIMPLE method is employed with 3 inner corrections and one outer correction.
Figure 2 shows the visual representation of the experiment extracted from Chatzis and Dullien [19].The water is white and the n-decane is black.The water imbibes from the right into the pore doublet and first displaces the ndecane in the smaller capillary followed by imbibition into the larger capillary.The results of OpenFOAM modeling are presented in the sequence of plots on the right where water and n-decane are depicted in blue and red, respectively.
As expected from the pore doublet experiment and also observed in the simulation, when the velocity is very low the situation becomes more similar to spontaneous imbibition and capillary force is the main contributing driving force for the system.Consequently, imbibition takes place in the smaller constriction first then displacement occurs in the larger diameter pore.The process continues until 100% displacement of the non-wetting phase by the wetting phase is accomplished, as observed from the experiment (Figure 2).It is evident that the numerical results obtained are similar to the experiment results.
However, it should be mentioned that at higher injection rates, as the viscous forces become comparable to the capillary force, the displacement in the pore with a larger diameter will start even before displacement in the smaller diameter pore is completed.The simulation results at very first stages of the imbibition process are presented in Figure 3.The front advances at both pores, and as the process continues the front in pore with larger diameter stops and it retreats back.Confirmation of the numerical results requires visual experiments to monitor the flow behavior at the early stages of its development.

Drainage Experiment
At the pore-scale, drainage is defined as the forced displacement of wetting fluid by a non-wetting fluid [21].It is expected that the invading fluid occupies and pushes the fluid in the pore with a larger radius because the capillary entry pressure is lower compared to a pore with smaller radius [19].The drainage boundary conditions are shown in table 6.All the settings for discretization and solution of equations are the same as the imbibition process which is explained in the imbibition section.The drainage experimental results in the pore-doublet are presented in Figure 4 (left) with n-decane (black), a nonwetting fluid, injected from the left draining the water (white) to the right (downstream) preferentially through the larger diameter capillary.The numerical simulation results are presented on the right.

Sensitivity Analysis on Wall Boundary Conditions
To evaluate the ability of interFoam on capturing wetting film, some sensitivity analysis is performed by considering different slip conditions and contact angles.Four cases are compared as shown in the following subsections.

Slip Condition on Walls
In this case, all the parameters and boundary conditions are the same, just slip conditions on walls are active.The slip boundary condition erases the normal component of the variable at the patch and keeps the tangential components untouched.The results do not show any obvious difference with no-slip as boundary conditions.
The movement of the front at the juncture of two conduit, however, becomes unstable which leads to a snapoff effect.The overall behavior of front movement remains the same as the previous case.Snapshots of the simulation results are shown in Figure 5.

Slip Condition with Dynamic Contact Angle
Here, a dynamic contact angle was set at the boundaries.The advancing and receding dynamic contact angles are set to 45 and 60 degrees respectively.Results do not show any obvious difference compared to previous case, i.e., slip condition on walls.Here, the receding and advancing contact angles are set to 10 and 15 degrees respectively, a spread of wetting film near the walls can be seen.This is expected because, by setting the contact angle as low as 10 degrees, we assumed the walls are strongly water wet.The sequence of front advancement, in this case, is also presented in Figure 6.The advancing fluid (water) starts moving in both pore conduits, but at some point, the front retreats and stays stagnant in the pore with a larger diameter.At the same time, the front at the pore with small diameter moves and pushes the n-decane out of the pore.The simulation results are comparable to the experiment.Large diameters control imbibition and the wetting fluid imbibes the pore with a smaller diameter first.Using snappyHexMesh the mesh is refined near the walls to better capture two-phase interface curvature near the wall.Consequently, the number of cells increased to 107891.The properties of cells are shown in Table 7 and a picture of the cells is shown in Figure 7.The dynamic contact angle with receding and advancing angles of 10 and 15 o is being used in the simulation.Threedimensional flow simulation, despite the huge increase in run time, is being applied.Some snapshots from simulation on interface movement are being represented.The effect of thin-film flow of water that spreads over the wall is evident.As the slip condition is prevalent near the wall and water is being considered as wetting fluid, thin-film extension near the wall is expected.This phenomenon is being predicted by simulation.As demonstrated, interFoam can capture the advancement of film flow near the walls due to the high wettability of water which is being set at boundaries. Figure 8 shows a thin film of water is being constructed while water imbibes throughout the pore doublet.The results are in agreement with the simulation obtained in case 3.5.3.Comparing the results of slip conditions with no-slip conditions at the rock surface, the progression of film flow cannot be captured in no-slip conditions.However, at low to medium wettabilities, these conditions have no preferences in determining the movement of fluid and interface in the pore-doublet system.For strong wettability, applying slip flow is important to predict the creation and evolution of film flow near the rock surface.The effect of film flow due to the geometry of pore interconnection and their tortuosity will be more pronounced even in fluids with weak to medium wettabilities.Hence, applying slip conditions at the rock surface is important to capture the entrapment of non-wet fluids.

Conclusions
The interFoam solver for immiscible two-phase flow was evaluated by predicting imbibition and drainage in a pore doublet and qualitatively comparing the results to the experiments carried out by Chatzis and Dullien [19].As expected, the wetting fluid first enters the smaller diameter pore, in imbibition, displacing the non-wetting fluid then continues to the larger diameter pore.The results from the experiments and simulation are, qualitatively, in good agreement with each other.The simulation shows initial displacement in the smaller pore, then front movement in the larger diameter pore due to the competition of viscous and capillary forces.The wetting fluid in the larger diameter pore then retreats to the pore inlet before displacement in the smaller diameter pore continues.These phenomena are analytically explained by Chatzis and Dullien [19].
Drainage results from the simulation and experiment are in agreement with each other.As expected, the displacement happens in the larger diameter pore first and bypasses the wetting fluid in the narrow pore.
The sensitivity analysis shows that instability of the front is mainly related to the slip condition at the walls, and spreading of wetting film near the wall is prevalent at very low contact angles or in other words, in advancing fluids with strong wettabilities.Film flow is predicted at low contact angles, as shown in sections 3.5.3and 3.5.4.Moreover, the effect of film flow especially in fluids with strong wettabilities can be predicted by applying slip conditions at the rock surface.
In this study, the ability of interFoam in prediction of imbibition and drainage for simple pore doublet system is evaluated.Numerical analysis must be conducted on a pore scale system, representative of real porous media for validation in the next step.
Overall, interFoam solver has the acceptable capability to follow the behavior of two-phase immiscible flow in either drainage or imbibition for pore doublet system.Other CFD tools can be evaluated and be compared with each other.

2 . 4 .
Adjust the time step based on the Courant Number.3. Use old time volumetric flux to solve the phase fraction equation.Update physical properties; density, viscosity, face densities, unit normal vector, and curvature.5. Solve the momentum and pressure equations using SIMPLE, PISO, or combined SIMPLE/PSIO methods.6. Move to the time step.

Table 4 .
Solution Method for Alpha, Pressure, Generalized geometric-algebraic multi-grid b FDIC: a faster version of the diagonal based incomplete Cholesky preconditioner

Fig. 5 .
Fig. 5. Advancement of the Imbibition Front in Pore-Doublet with Low Dynamic Contact Angle

Fig. 7 .
Fig. 7. Spreading of Wetting Film Near the Walls during Imbibition (right to left)

Table 1 .
Mesh Quality used for CFD Study

Table 2 .
Physical Properties of the Fluids

Table 3 .
Boundary & Initial Conditions in Imbibition Simulation

Table 3 .
Discretization Schemes used for the Set of Equations

Table 6 .
Mesh Properties Generated for Three Phase Flow with snappyHexMesh