Modeling multiphase flow with a hybrid model based on the Pore-network and the lattice Boltzmann method

Lattice Boltzmann method (LBM) simulations provide an excellent description of two-phase flow through porous media. However, such simulations require a significant computation time. In order to optimize the computation resources, we propose a hybrid model that combines the efficiency of the pore-network approach and the accuracy of the lattice Boltzmann method at the pore scale. The hybrid model is based on the decomposition of the granular assembly into small subsets, in which LBM simulations are performed to determine the main hydrostatic properties (entry capillary pressure, capillary pressure liquid content relationship and liquid morphology for each pore throat). A primary drainage of a random packing of spheres is presented and contrasted to the results of the same problem fully resolved by the LBM. Liquid morphology and invasion paths are correctly reproduced by the hybrid method.


Introduction
The mechanical behavior of partially saturated granular systems is strongly influenced by the fluid flow. Several works [1,2] include hydro-mechanical couplings in large grain-fluid systems by means of the Discrete Element Method (DEM). However, the hydromechanical response of such studies was restricted to low liquid content to ensure the pendular regime, where liquid is retained in the form of bridges between particles. When the liquid content is increased, pendular bridges coalesce forming complex liquid morphologies. At this point, the pendular regime is replaced by the funicular regime. It is challenging for the hydro-mechanical couplings to capture the geometry of phases and interfaces for the whole range of saturation regimes. Even though 3D images from X-ray tomography are a useful and common tool to characterize the fluid morphology in the porous medium during the funicular regime [3], very few attempts have been carried out to obtain a geometrical description of the liquid distribution and morphology within the porous media [4,5].
Multiphase flow within porous media can be described solving the governing equations. However, the complex geometry and non-linearity of the problem hinders the analysis without crude approximations. Accordingly, most of the approaches rely on geometrical simplifications and transport properties are taken from experiments in laboratories. Alternatively, several numerical models can be employed to reproduce the fluid displacement within the porous media. Lattice Boltzmann method [6], volume of fluid (VOF) method [7] or Lagrangian mesh-free methods [8] are examples of micro-scale continuum models widely adopted to predict multiphase flow through porous media. Despite the high computational cost, these models have been proven successful to reproduce the physical phenomena at different scales. Besides, the nature of these models make them suitable for complex geometries and benefit from parallel computing. Finally, pore-scale models arise as an alternative to macro-scale and micro-scale continuum models. Pore-Network models manage to reduce the high complexity of porous morphology by a discretization of the pore bodies and pore throats. The hybrid model introduced in this work combines the efficiency of the Pore-Network approach and the accuracy of the LBM to evaluate the local rules at the pore scale [9].
The hybrid method presented in this article follows the pore-scale approach referred as "two-phase pore-scale finite volume-discrete element method" (2PFV-DEM) [10]. In this method, the invasion of the non-wetting phase is controlled by the entry capillary pressure (p e c ). When the capillary pressure (p c ) exceeds the entry capillary pressure (p c > p e c ), the non-wetting phase passes through a narrow throat invading a pore body. p e c can be predicted by the Incircle method [11], the Mayer-Stowe-Princen (MS-P) method [10], or direct fluid simulations. The hybrid model follows the decomposition scheme proposed by [12]. Such technique leads to a discretization of the pore space as a set of connected throats. LBM simulations are employed to determine the main hydrostatic properties of each throat (p e c , liquid morphology, capillary pressure -saturation relationship).

Lattice Boltzmann Method
The multicomponent multiphase Shan-Chen LB model [13] is implemented using the open source library of Palabos. The fluid motion is governed by the lattice Boltzmann equation: (1) where the superscript ω denotes the ωth fluid, f ω k is the particle distribution function and the main variable of the LBM equation, f ω k gives the probability to find a particle in a certain node (x k ), time (t) and specific velocity (e k ). τ ω is the rate of relaxation towards local equilibrium, f ω,eq k is the equilibrium distribution function, ∆t is the time increment, e k are the discrete velocities which depend on the velocity model, in this work, D3Q19 (three-dimensional space and 19 velocities) model is adopted, and k varies from 0 to Q − 1 representing the directions in the lattice. The left-hand side of Eq. 1 corresponds to the streaming step, which describes the motion of the fluid particles between nodes, whereas the right-hand side stands for the collision step based on the Bhatnagar-Gross-Krook (BGK) approximation. The collision operator represents the evolution of the system towards the equilibrium.
The macroscopic density and momentum variables are recovered from the distribution functions: In order to ensure a system of immiscible fluids, a repulsive force between the fluid phases is introduced. According to [13], the non-local force responsible for the fluid-fluid interaction can be expressed as: where Ψ k is the interparticle potential that induces phase separation and G ωω is the interaction strength between components ω,ω.
Finally, the non-ideal equation of state can be expressed as: where ρ ω and ρω are the densities of the fluids at a certain position and c s is the speed of sound.

Hybrid model 3.1 Pore space decomposition
Following the 2PFV-DEM scheme developed by [10] based on a three-dimensional triangulation method and DEM to study the hydro-mechanical behavior of unsaturated deformable granular materials, the hybrid model relies on the same triangulation to generate the pore-scale network. This decomposition technique leads to a tetrahedral mesh whose vertices coincide with the centers of the spheres. Each tetrahedral element defines a pore body and four pore throats. Therefore, the pore space is decomposed so that pore throats can be solved independently. The LBM computation domain for each pore throat is a triangular-shaped prism defined by three solid walls orthogonal to the pore throat. Each of these solid boundaries passes through two of the spheres centers (two vertex of the triangle defined by the 3 spheres centers). We refer the reader to figure 1 for a comprehensive review.

Fluid displacement
The fluid-fluid interface displacement in the porous media is mainly controlled by the entry capillary pressure (p e c ). The interface curvature increases as the capillary pressure builds up following the Laplace-Young equation.
When the local capillary pressure is larger than p e c , the non-wetting phase passes through the pore throat invading the pore body.
In the hybrid model, the pore space is decomposed into a series of throat-domains defined by a regular triangulation (see figure 1). LBM simulations are performed for each pore throat to predict the primary drainage curve, p e c and liquid morphology. The computational domain displayed in figures 1b and 1c is enclosed by two triangles at the top and bottom of the prism representing the inlet and outlet sections of the LBM simulations. Initially, both wetting and non-wetting phases are in equilibrium for each elementary problem. Figure 2 shows the evolution of a typical subset. In this case, the pore throat is formed by three equal-sized spheres in contact. The mass source and sink at the nodes located in the inlet and outlet sections are imposed to gradually increase the capillary pressure. Such increment triggers the displacement of the fluid-fluid interface (see figure 2). The evolution of normalized capillary pressure (p * c = 2p c R γ , where γ is the surface tension, R is the radius of the spheres and p c is the capillary pressure) during the invasion is related to the change of volume of the wetting phase in figure 2d. Dimensionless volume is defined as: The entry capillary pressure (p e c ) is defined as the maximum value attained by p c during the invasion process. When the capillary pressure reaches the entry capillary pressure (p c = p e c ), the non-wetting phase penetrates into the pore body (see figure 2). Immediately after the invasion, the interface meniscus expands leading to a reduction of the capillary pressure. This procedure is repeated for all the pore throats of the granular assembly. Thus, p e c is determined for all the subdomains and incorporated into the global problem handled by the pore-network.

Numerical setup
A random sphere pack is created by the open-source code YADE [14]. A cubic box of 10 mm x 10 mm x 10 mm is defined in which 40 polydisperse spheres are packed. The mean sphere radius is 1.26 mm. We assume: − Negligible gravitational forces.  A fully saturated granular assembly is considered as initial state (see figure 4a -Full LBM). In order to reproduce a fully resolved LBM drainage simulation of the sample, a porous membrane is located at the bottom of the sample in order to prevent the non-wetting phase to reach the outlet and ensure a complete drainage. Additionally, the displacement of the interface is driven by a fixed flow controlled by a pressure difference between the non-wetting reservoir (top of the sample) and the wetting reservoir (bottom of the granular assembly). The interface is displaced towards the bottom driven by the pressure gradient invading the largest pores at first. Then, the non-wetting phase progressively occupies narrower pore throats of the sample. Capillary pressure is increased and recorded until all the nodes above the porous plate are filled with the non-wetting phase (see figure 4d-Full LBM). At this moment, the remaining wetting phase is trapped in the granular assembly in the form of liquid clusters. Figure 3 evidences the presence of pendular bridges and trimers (liquid cluster formed between three grains) in the sphere packing after the drainage.

Liquid morphology and preferential paths
This section illustrates the invasion paths of the nonwetting phase through the porous media and the evolution of liquid morphology during drainage. Figure 4 captures the non-wetting pathways for different time steps for each method. The full LBM interface is shown as translucent blue isosurface in figure 4 while the wetting phase is depicted in opaque blue for the hybrid model. Both methods (fully resolved LBM and hybrid model) presented here follow the same first preferential path (identical intrusions are observed in figures 4b and 4c). As the drainage proceeds, the non-wetting phase keeps occupying the sample leaving behind some disconnected liquid clusters (see figure  4d). Even though wetting phase trapped in the assembly is well represented with the hybrid method, figure 4d suggests that there are some discrepancies between the models in some regions. Overall, we observe that the non-wetting phase does not happen as a stable piston-like front. Instead, both methods demonstrate a fingering mechanism that leaves pendular bridges and other cluster behind after the main invasion event. In terms of liquid morphology, figure 5 compares details of the fluid-fluid interfaces obtained from the full LBM (translucent interface) and the throat-scale simulations (opaque interfaces) for a given p c . Meniscii from the two methods are nearly identical, which validates the choice of regular triangulation to decompose the pore space.
It is worth mentioning that the CPU time for a complete drainage corresponds to 29.6 days for the full LBM. On the contrary, by using the hybrid model, the computational cost is reduced down to 11.2 days.

Conclusions
Compared to a fully resolved method (in this work: full LBM), the hybrid model presented in this article provides more efficient results in terms of computation time. Results provided by the hybrid method are able to mimic the interface displacement during the drainage of pore space with excellent accuracy. This work evidences that tracking the p c − V relationships for each pore during drainage is crucial to determine the main hydrostatic properties accurately. The total liquid content for a certain capillary pressure is recovered as E3S Web of Conferences 195, 02009 (2020) E-UNSAT 2020 https://doi.org/10.1051/e3sconf/202019502009 the sum of filled tetrahedra plus the meniscii volumes. The later volumes are frequently neglected or approximated. In the hybrid model, however, the volume over each tetrahedron facet is easily accessible by means of the p c −V curves for the corresponding pore throat. Additionally, the strategy adopted by the hybrid method can shed some light on the analysis of complex phenomena such as the pore refilling. Finally, the hybrid model allows knowing the geometry of interfaces at any moment of the drainage process. This feature is crucial to determine the forces on the solid particles, which may be incorporated in hydro-mechanical couplings.
This technique would probably have a stronger impact on a lager problem with thousands of pore throats. In parallel simulations with relatively small computational domains, the execution speed is multiplied by a factor two doubling the number of cores. Large-scale problems that require the exploitation of massively parallel systems are frequently accompanied by a drop in efficiency. In the hybrid method, however, the scalability is optimal. Due to the fact that each pore throat are assigned to a single core, efficiency does not depend on the size of the granular assembly. Besides, some regions of the sample can be excluded from the simulation (i.e. empty pores and isolated cluster with no flux).