Benchmark cases for a multi-component Lattice–Boltzmann method in hydrostatic conditions

Graphical abstract


Introduction
Fully resolved numerical solutions to multiphase pore-scale problems are used increasingly in simulation domains extracted from 3D imaging. There is, simultaneously, a growing interest in the development of simplified methods based on pore-network idealizations since simulating spatiotemporal evolutions at the REV scale would require tremendous computational resources when intricate couplings are at play. This was done primarily for low-porosity materials (typically rock materials) [1,2] . Extensions of the pore network approach to granular media appeared more recently and they still rely on strong assumptions and simplifications [3][4][5] . In the pore-network approach the movements of phases and interfaces are governed by local rules such as the entry capillary pressure, the capillary pressure -saturation curve and the capillary forces. When the local capillary pressure is larger than the entry capillary pressure of a pore throat the non-wetting phase penetrates it invading the pore body. Several approaches can be considered to compute the entry capillary pressure. The most common approximations are the Haines incircle method and the Mayer-Stowe-Princen(MS-P) method [3][4][5] . Unfortunately, these approximations predict just a single pressure value missing crucial information before and after the invasion that could be provided with an accurate local capillary pressure -saturation relationship. Establishing those local rules is another use-case for fully resolved solutions -for elementary microstructures in that case. The lattice Boltzmann method (LBM) is frequently used for producing well resolved solutions. In this study we assess the accuracy of a multiphase LBM scheme for the solution to hydrostatic problems. A background motivation of this work is the extension of pore-network methods to deformable granular media, following the strategy employed previously for saturated flow [6] . We therefore focus on elementary microstructures. Nevertheless the conclusions in terms accuracy and mesh dependency apply equally well to simulations of REVs. It is, thus, worth mentioning that this benchmark is intended to serve as validation of the numerical simulation method to be applied in a practical situation. More specifically, this benchmark is used to justify the mesh resolution and flow conditions employed in [7] , where the pore space is decomposed into small subsets of three spheres (pore throats) that are solved independently to determine the main hydrostatic properties.
The LBM is a mesoscopic model capable of simulating fluid dynamics in complex geometries [8] . Many works using the LBM have focus on a single saturating fluid phase and proven to be successful [9][10][11] . However, multiphase LB models in partial saturation have less satisfactory results due to the complexity of phases interactions. Several multiphase LB models have been proposed in the literature: the color model [12] , the pseudopotential (Shan-Chen) model [13,14] or the free-energy model [15] . The so-called Shan-Chen model has single-and multi-component variants which both apply to the problem of immiscible phases. The single-component method is simpler. It has been used to simulate, for instance, flow in porous media with realistic rock geometries [16,17] or the hysteretic response of idealized sphere-pack systems in drainage-imbibition [18] . More recently, [19] investigated with this method the meniscus profile and the effect of contact angle on fluid displacement through polygonal capillary tubes. According to [20] however, the gas-liquid interfaces tend to be more diffused in single component simulations, which may hinder the approach of strongly immiscible situations. Fewer studies have applied the multicomponent method [21,22] although it is supposed to reflect the fluid composition more accurately [23] . Very few authors -if any -examined the accuracy of the multicomponent scheme for hydrostatic solutions. In this paper, the multicomponent Shan and Chen model is employed using the open-source lattice Boltzmann library Palabos [24] to complement the results and conclusions of previous studies and benchmarks [23,25] .
The paper is organized as follows: in section Numerical method , the lattice Boltzmann method and the Shan-Chen model are briefly described; section Model calibration explains the way that surface tension and contact angle can be computed and tuned; in the section Validation LBM results are compared to analytical solutions for capillary tubes and pendular bridges between spheres; finally, conclusions are drawn in last section.

Lattice Boltzmann method
In this section we provide a brief explanation of the LB method. The LBM has its origin in the lattice gas automata (LGA) [26] , a kinetic model based on discrete space-time field. While LGA method described the evolution of individual particles on a lattice, the LBM solves a discrete kinetic equation (Boltzmanns equation) for a particle distribution function f σ ( x , t ) . Where the superscript σ indicates the fluid component, x refers to the lattice node and tis the time. In the LBM, the motion of fluid is described by the lattice Boltzmann equation. Based on the simple and popular Bhatnagar-Gross-Krook (BGK) collision operator [27] , the standard LB equation can be expressed as follows: where τ σ is the rate of relaxation towards local equilibrium, f σ,eq k is the equilibrium distribution function, tis the time increment, e k are the discrete velocities which depend on the particular velocity model, in this work, D3Q19 (three-dimensional space and 19 velocities) model is used, and k varies from 0 to Q − 1 representing the directions in the lattice. The left-hand side of Eq. (1) describes the streaming step (particles move to the nearest node following its velocity direction) whereas the right-hand side stands for the collision operator (particles arriving to the nearest node modify their velocity towards a local equilibrium). The collision operator correspond to the viscous term in the Navier-Stokes equation. For the D3Q19 model, the discrete velocity set e k is written as: where w k are the weight factors.
The local equilibrium f σ,eq k depends on the lattice type and the macroscopic variables ρ σ = k f σ k (density) and ρ σ u σ = k f σ k e k (momentum) [28] . The equilibrium distribution can be seen as an expansion of the Maxwell-Boltzmanns distribution function for low Mach numbers: is the speed of sound and u σ,eq is the equilibrium velocity defined as [13,14] : where u = σ ρ σ u σ τ σ σ ρ σ τ σ is an effective velocity and F σ is the total force (including body forces and the fluid-fluid interactions that will be presented in Pseudopotential model section) acting on each component.

Pseudopotential model
The interactions between components (or phases) in the Shan and Chen model are defined by pairwise interaction forces. These forces modify the collision operator through an equilibrium velocity and produce a repulsive effect between the phases. We focus on biphasic mixtures (i.e., σ = 1,2), Hereafter, ρ w and ρ nw will refer to the wetting and non-wetting phases. ρ o is defined as the reference density which is kept at ρ o = 1 . The non-local force responsible for the fluid-fluid interaction is expressed as: (5) where k is the interparticle potential that induces phase separation and G σσ is the interaction strength between components σ, σ .
Previous works [14,[29][30][31] have employed several interparticle potentials. For simplicity, we consider k = ρ k , as done in other papers [18] . The interactions within each component, G 11 and G 22 , are set equal to zero for biphasic mixtures. On the other hand, the interactions between components, G 12 = G 21 , are set positive in order to induce a repulsive force between the phases. Low values of G 12 lead to dissolution processes seen in typical miscible mixtures. On the contrary, significantly high values of G 12 result into almost immiscibile binary mixtures with sharp interface prone to numerical instability. Thus, special attention must be paid when choosing the interaction strength as it controls the surface tension and immiscibility of the mixture. The interaction force given by Eq. (5) leads to a non-spherical pressure tensor P deduced from the condition: [32,33] . The components of the pressure tensor can be computed as: Following Eq. (6) , the non-ideal equation of state (EOS) can be determined as: Model calibration

Contact angle
The fluid-solid interaction is implemented in the Shan-Chen model by a mid-grid bounce back scheme applied on the boundaries [34] . This scheme assigns fluid properties to the solid wall. Among them, the pseudo wall density ρ wall (non-real density assigned to the nodes of the solid boundary) controls wettability [19,35,36] . The interparticle potential at the wall in Eq. (5) is = ρ wall ). We perform simulations of static droplets on a flat solid surface and we analyze the dependence of ρ wall on the contact angle. Simulations are performed in a 150 ×150 ×150 lattice domain. Once the simulation is stable and converged, the base length ( b) and the height ( h ) are measured. Knowing the geometrical characteristics of the droplet allows us to determine the contact angle θ [37] (see Fig. 1 (b)). Some error is introduced during the base measurement due to the thickness of the interface layer in the vicinity of the solid wall. In order to overcome the problem, the base and height of the droplet are determined from a reference point located 2 lattice units away from the wall ( Fig. 1 (a)). Moreover, as further discussed in Numerical method section, ρ w /ρ o = 0 . 7 is the density threshold used for positioning the interface (dark line in Fig. 1 (a)).

Surface tension
Surface tension is adjusted by tuning the interaction between different fluid species. The typical numerical set-up to investigate the surface tension consist of a series of spherical drops with different  radii inside a domain with periodic boundary conditions. The droplet and the surrounding fluid are at rest and the pressure difference inside and outside the droplet is balanced by the surface tension according to the Young-Laplace law ( p c = 2 γ R ). Fig. 2 (a) depicts the pressure along a line passing through the center of the droplet. There are two significant drops in pressure when the line crosses the interface [38] , it denotes to surface tension. The pressure difference pcorresponds to capillary pressure. Fig. 2  Surface tension can be also be determined based on a two-phase system with a flat interface having a constant pressure in both phases far from the interface [39] . This technique has been adopted in many works relying on the single-component Shan-Chen model [14,32,38,40] . Literature on the multicomponent model is more scarce yet the flat interface has also been used in that case [35] . We reproduced it for comparison with the droplet test. The logic of the analysis is as follows. The pressure inside the bulk phases corresponds to the scalar quantity p. However, near the interface, due to the surface tension contribution, the pressure is defined as a tensor incorporating different pressure components. Moreover, in order to ensure the mechanical stability, the gradient of the pressure tensor must be zero everywhere in the fluid [41] . The symmetry of the surface requires that pis a diagonal tensor p(x ) = p xx e x e x + p yy e y e y + p zz e z e z with p xx ( x ) = p zz ( x ) , where x and zcorrespond to horizontal directions parallel to the flat interface, y refers to the axis orthogonal to the planar interface and e j is a unit vector in the j-direction. Furthermore, p xx and p zz are function of y only, while p yy is a constant: where p T and p N are the transverse and normal components of the pressure. Both p T and p N can be computed using Eq. (6) .
Surface tension is obtained by integrating the difference between p T and p N along a line crossing the interface [39] : The results from droplet test and the flat interface test are compared in Fig. 3 , they are in good agreement.

Note on interface thickness
The numerical thickness of the interfaces, as seen in Fig. 1 is often considered an issue in the multicomponent Shan-Chen model. Physically inter-molecular interactions lead to a fluid-fluid interface thickness, i.e. a region where the two phases coexist even though they are considered immiscible from a macroscopic point of view. On this basis the fact that the multicomponent Shan Chen model produces diffused interfaces is not strictly unphysical (see Fig. 1 (a)). In many applications however the real interface thickness is well below all characteristic lengths of the problem (such as pore size or radius of curvature), hence negligibly small, and then the interface is considered a single surface. In LBM however the thickness of the simulated interface does not correspond to the physical thickness in general. Previous works [42,43] have evidenced that a fluid-fluid interface of 4-6 lu is required for numerical stability, which could be neglected only at the price of extreme mesh refinement and tremendous computational effort. Some works [23,25] have attempted to increase the accuracy at fluid-solid interface by introducing new boundary models. Despite the effort s and the better results obtained near the solid region, numerical artifacts are still found to decrease the global accuracy. In order to overcome this issue we propose to redefine the solid boundaries based on a wall retraction logic, including a part of the fluid-solid interface in the region normally occupied by the solid phase in the physical problem (as shown in Fig. 4 ). This is tested in the next section in the context of capillary tubes.

Validation
Simple numerical simulations are performed and compared with analytical solutions in order to validate the model. Detailed results are presented for quasi-static displacement of interfaces inside cylindrical tubes and fluid bridges between two spherical bodies.

Invasion of capillary tubes
In order to gain better understanding of multiphase flow at the pore scale, it is common to idealize the pores throats as cylindrical capillary tubes [44] . Immiscible flow in such capillary tubes has been simulated with various cross-sectional shapes ( Fig. 5 ). The dimensionless capillary pressure p * c = p c L c γ is defined with reference to the following characteristics lengths: L c is the radius for the circular cross-section, the side length for the square, the distance between two vertices for the triangle and curved triangle. The fluid displacement corresponds to drainage (invasion by the nw -phase) and it is imposed by including mass sink terms in the time integration: wetting phase density is decreased while non-wetting phase density is increased [18] . In order to keep the flow quasistatic the density is only modified when its fluctuation on one time iteration, at interface nodes, is less than a fixed . Otherwise the solution is considered out-of-equilibrium and the mass sink is delayed.

MS-P method
The Mayer and Stowe-Princen (MS-P) model predicts the capillary pressure and the curvature of the arc meniscus of a fluid droplet of infinite length inside a cylindrical tube [45][46][47] . The assumptions of the MS-P method are that that capillary pressure is uniform and that there is no longitudinal curvature away from the main terminal meniscus. Under these assumptions the cross-sectional radius of curvature R (see Fig. 11 ) defines the total curvature and, after Young-Laplace equation, Furthermore, the balance of forces at equilibrium implies a relationship between capillary pressure and surface tension. The force due to the pressure difference on the cross-sectional area must balance the force from surface tension at the interfaces. Thus, p c A nw = γ (P s cosθ + P ns ) (12) where P s is the length of the line between the non-wetting phase and the solid, P ns is the perimeter of the interface between the wetting phase and the non-wetting phase, and A nw is the area filled with the non-wetting phase. The MS-P method consist in deducing R by combining Eqs. (11) and (12) :

R =
A nw P s cosθ + P ns (13) From now on the MS-P is considered exact for cylindrical throats and used as a reference for comparisons. The errors in LBM solutions will be evaluated using two possible approaches: where p LBM e is the entry pressure obtained in the saturation curves ( Fig. 7 ).
where k MSP is the curvature defined by the MS-P (the inverse of the radius of Eq. (13) ) compared with the curvature of the main meniscus after achieving the entry pressure. k LBM is defined in Appendix A .

Results
The entry capillary pressure p LBM e in the LBM simulations is deduced from drainage curves similar to the plots in Fig. 7 , where V is the volume occupied by the wetting phase within the tube. The  Fig. 6 . Note that the difference is not expected to vanish even with very small tolerance since geometrical discretization errors adds to the error relatively independently of dynamic effects. In the sequel of this study we set the tolerance value to 10 −5 , as it leads to marginal dynamic errors.
Several mesh discretizations have been tested: 40 ×40 ×160, 70 ×70 ×256, 90 ×90 ×320 and 110 ×110 ×384 (last value along the axis of the tube). From now on they are referred to as L c = 40 l u, L c = 70 l u,L c = 90 l u, and L c = 110 l u, respectively. The pressure-volume evolution for each mesh size are compared in Fig. 7 (for the square-shaped tube). The errors with respect to the MS-P prediction are given by Fig. 8 . When the numerical solid wall coincides with the physical wall (no wall retraction) the convergence is superlinear, with an exponent of approximately 1.4. When the interpretation includes wall retraction by two lattice units, the error is smaller and the convergence Fig. 8. Convergence of the LBM result with mesh refinement, with regard to error defined in Eq. (15) . Each simulation is ran in parallel using 8 cores. L c is defined as the distance between the numerical walls (unchanged by wall retraction). becomes quadratic, which is a substantial improvement. This technique was used systematically for all simulations presented in the next sections.
A justification of the optimal retraction length is possible by selecting different iso-density surfaces in the result to represent the interface. A consistent definition of the interface should satisfy Eq. (12) .
Selecting a value of ρ w /ρ o to define the interface enables the determination of the geometrical parameter A nw ,P s , and P ns in that equation. The optimal contour is the one which minimizes the deviation from Eq. (12) . Based on Fig. 9 the optimum is ρ w /ρ o = 0 . 7 , which corresponds approximately to the average density between both phases. In our results this specific value of density was generally reached approximately two nodes away from the solid nodes, which led to the decision to retract the walls by two lattice units. This value is only valid for Gρ o = 1 . 25 . Different interaction strength parameters (i.e. other surface tensions) would result in thicker or thinner interfaces, in such case, the same procedure should be repeated to determine the position of the new retraction wall.
The various cross-sectional shapes have been simulated with domain sizes 80 ×80 ×256 lu 3 . The results are compared to the MS-P solution in Fig. 10 . We find a reasonable agreement between the  simulations and the analytical solution overall. However, larger errors are observed for triangular and curved cross-sections. This can be partly attributed to the artificial roughness introduced by the staircased surfaces. These cross-sections are not aligned to the regular lattice grid. Furthermore, due to the bounce-back boundary condition, these cases lead to mesh-dependent results. In fact, an asymmetry is evidenced in Fig. 11 , where the remaining liquid retained in the corners of the equilateral triangle is different in some parts. Nonetheless, Fig. 11 shows relatively similar numerical and analytical profiles.
This mesh dependency is frame dependent: it depends on the orientation of the throat with respect to the axis of the grid. The evolution of the errors with rotation is shown in Fig. 12 , which reveals that the frame-dependent effects are actually small (of the order of 1%, dominated by other errors).  To conclude this section, we review the hypothesis stating that MS-P solution is valid for cylinders of infinite extension. Due to computation limitations, short domains had to be considered. In order to test the accuracy of the numerical results under these conditions, the error on pressure has been plotted along the cylinder. In other words, capillary pressure was computed using Eq. (12) for various positions of the cross-section in the final, nearly fully invaded, configuration. On the left part of Fig. 13 we observe that the remaining fluid in the corners is parallel to the cylinder walls (no longitudinal curvature). It is concluded that H/L > 1 is sufficient to approach the situation assumed for the MS-P method, i.e. the cross-section must be behind the main meniscus by a distance approximately equivalent to the throat aperture.

Pendular bridge
The shape and volume of a pendular bridge between two spheres have been obtained from the LBM and compared to the theoretical solution for a range of capillary pressure.
The simulation setup was as follows: a droplet of the wetting phase was inserted between two identical spheres of radius R with a gap equal to 0 . 14 × R . Once a stable state was reached the volume of the liquid bridge was reduced slowly, by an imposed mass sink, until p * c = p c R γ = 0 . 3 . The shape of the pendular bridge when p * c = 0 . 3 is compared to the direct solution of Young-Laplace equation [48] in Fig. 15 . They show strong similarity. After reaching p * c = 0 . 3 the LBM simulation was continued by further reducing the amount of wetting phase and recording the volume of the simulated bridge for quantitative comparison with Young-Laplace solution. This was continued until breakage of the bridge. Fig. 14 shows the volume-pressure dependency until breakage. The LBM simulation and the Laplace-Young solution follow a very similar trend, with the relative error generally less than 10 −2 .
Likewise, the critical distance S c (sphere separation that leads to breakage of the bridge) can be compared. S c can be obtained on a theoretical basis: it is the distance beyond which the Laplace-Young problem degenerates into a solutionless problem (practically approached by the upper bound of the actual solutions). Previous works [48] have shown that S c is approximately proportional to the cubic root of the volume of the bridge. This empirical relation is also compared to the results. Fig. 16 shows the rupture distance obtained by the different methods. The LBM follows a correct trend yet the distance is systematically underestimated, by 4% approximately. It is less accurate than the cubic approximation. The systematic underestimation can be explained by the difficulty to approach a mechanically unstable solution numerically.

Conclusions
The hydrostatic properties and pore-scale morphology of immiscible phases have been obtained by the multicomponent Shan-Chen LBM for systematic comparisons with other methods. This article provides estimates of discretization errors and guidelines to calibrate the method and minimize errors.
Two-fluid-phase flow through capillary tubes has been analyzed and compared to the solution given by the MS-P method. Entry pressure, curvature and interface profile obtained from LB simulations converge to the analytical solution with mesh refinement. The capillary bridges simulated between 2 spheres also converge to the solution obtained directly from Laplace-Young equation, in terms of both shape and rupture distance. Discretization errors are introduced in part because of the solid boundaries: curved surfaces are modeled as stair-cased lines, which may not approximate the  curved wall properly if the lattice resolution is not fine enough. In addition the numerical thickness of the fluid interfaces around the solids is also a source of error. These discretization errors were found to scale nearly linearly with mesh size, and relatively independently of rotations of the grid frame. For the error due to interfacial thickness we showed (section Results ) that a significant reduction was possible with appropriate geometrical corrections of the solid boundaries. This correction leads to shrink the size of all solid objects by a mesh-dependent length to minimize the mesh-dependency of the result. This technique has been used systematically throughout this study and proved to give satisfactory results.
The aim is to progressively improve the local rules introduced in pore-network approaches from the analysis of elementary subsets, following [49] . Indeed, this article is meant to be a validation of the multicomponent Shan-Chen model to simulate multiphase flow in porous media and justify the mesh resolution and flow conditions used in [7] where an sphere packing is decomposed into a series of subsets that are solved separately using the LBM.

Declaration of Competing Interest
The Authors confirm that there are no conflicts of interest.

Physical and LBM units
Correlating physical properties to lattice units is an essential task in order to simulate physical problems. Moreover, choosing the right conversion will avoid stability problems and help us to have accurate results. As suggested in [50][51][52] , physical units can be related to lattice units through unit conversion or dimensionless numbers such as the Reynold, the Froude or the Bond number. The parameters involved in the physical and the LB systems are summarized in Table A.1 . Conversion factors for length, time, velocity and density are: . Similarly, we can find expressions for the kinematic viscosity C ν = C 2 x /C t , and pressure C p = C 2 Kinematic viscosity is also related to the relaxation time τ as The method presented above is consistent and can be applied to find other quantities [51] . Nevertheless, one important constraint must be kept in mind. LBM is limited to low Mach numbers due to compressibility effects that lead to numerical instabilities [50,53] . In order to conduct numerical simulations of quasi-compressible flows and reduce the numerical error, lattice Boltzmann velocities should be significantly smaller than the speed of sound ( v lb << c s ). Dimensionless numbers are extensively used to overcome this limitation. The first step consist of converting the physical system into a dimensionless system. After that, dimensionless units are transformed into lattice units. For the sake of clarity, let us use the Bond number to illustrate the unit conversion in terms of

Curvature analysis
In order to analyze the multiphase flow it is crucial to study the shape of the fluid-fluid interface. Thus, in this section we introduce a method to determine the interface curvature following [54] .
Given a fluid-fluid interfacial surface Senclosed by an arbitrary volume element V such as the one displayed in Fig. A.1 , we can perform a force balance on V : Inertial force = Body force + Hydrodynamic force exerted on S + Surface tension force exerted along C , equivalently where, • dlrefers to a length increment along the closed curve Cthat forms its boundary (see Fig. A.1 ), • ρis the fluid density, • γ is the surface tension,