Toward multiscale modelings of grain-fluid systems

Computationally efficient methods have been developed for simulating partially saturated granular materials in the pendular regime. In contrast, one hardly avoid expensive direct resolutions of 2-phase fluid dynamics problem for mixed pendular-funicular situations or even saturated regimes. Following previous developments for single-phase flow, a pore-network approach of the coupling problems is described. The geometry and movements of phases and interfaces are described on the basis of a tetrahedrization of the pore space, introducing elementary objects such as bridge, meniscus, pore body and pore throat, together with local rules of evolution. As firmly established local rules are still missing on some aspects (entry capillary pressure and pore-scale pressure-saturation relations, forces on the grains, or kinetics of transfers in mixed situations) a multi-scale numerical framework is introduced, enhancing the pore-network approach with the help of direct simulations. Small subsets of a granular system are extracted, in which multiphase scenario are solved using the Lattice-Boltzman method (LBM). In turns, a global problem is assembled and solved at the network scale, as illustrated by a simulated primary drainage.


Introduction
The simulation of wet granular materials with pendular bridges between pairs of solid grains has been a popular application of discrete element models (DEM). The study of mixed situations in which pendular bridges and larger clusters of wetting phase co-exist (the typical situation), eventually reaching one single cluster in the saturated regime, has received much less attention. This is certainly due to the modeling and computational difficulties more than the lack of scientific an practical interest. An hybrid modeling framework taking advantage of both the efficiency of the pore-network approach and the accuracy of direct fluid dynamics solutions at the microscale is proposed in this paper. In the first section physical processes are associated to the basic geometrical entities of a triangulation and the main governing equations of fluid-grain mixtures are listed. The saturated case is included since one-phase flow in large single-phase clusters is intrinsically linked to the more general multiphase flow problem. The last section describes how improved accuracy is gained by embedding direct simulations of multiphase events into the pore-network model. The network representation of the pore space in three dimensions is obtained by Regular Triangulation [1] of the sphere assembly. The algorithm follows [2] and the implementation is based on the open source library CGAL [3,4]. e-mail: bruno.chareyre_a_3sr-grenoble.fr In the triangulated geometry, a pore body is defined as the irregular cavity within a (3-D) tetrahedron, surrounded by four solid spheres (whose centers are the vertices of the corresponding tetrahedron, see Fig.1(a)). A pore throat is defined by the cross sectional area extending within a tetrahedral (2-D) facet ( Fig.1(b)). Pair interactions (including solid-solid contacts) act along the (1-D) edges of the triangulation. Finally the 0-D points are attached to the solid objects to which total forces are assigned (Fig.  2). This Lagrangian mesh follows the solid grains.

One-phase flow
The pressure of a fluid flowing within a granular medium may not immediately appear as a discontinuous problem yet the changes of pore pressure are strongly localized near the pore throats (i.e. the facets of the triangulation), as shown in Fig. 3. This feature has been exploited by the  DEM-PFV scheme [2,5]. In this method the fluid mass balance in one pore links interface fluxes with the rate of changeV i of the volume of pore i (a function of particle velocities). In addition the linearity of the flux-pressure relations (Stokesian regime) leads to the discrete form of a Poisson problem with theV i as source terms, for i iterating over the poresV (1) K i j is the local hydraulic conductivity reflecting the local pore geometry [6]. Solving this linear system gives the p i , from which drag forces on the grains are evaluated and added to the total forces (strong two-way coupling). The lubrication forces, short-range hydrodynamic interactions between immersed particles, are not captured by the pore-scale mass balance. They should be added independently. The normal and shear components of the lubrication force read [7,8] (the rotational terms are omitted for brevity, see [9]) where v n and v t are the normal and tangential relative velocity, μ is the fluid viscosity, and a is the radius of the grains.

Immiscible fluids and interfaces
We now consider the case of two fluid phases in the pore space: a wetting phase and a non-wetting phase (hereafter w-phase and n-phase respectively) such as water and air, typically. For the sake of simplicity we admit that the particles are sufficiently small to make the gravitational effects negligible compared to surface tension effects (small values of the Bond number B o ). The geometry of the nw-interfaces follows the the Laplace-Young relation between the trace C of the curvature tensor and the capillary pressure p c (pressure difference between the n-phase and the w-phase locally): p c = Cσ wn with σ wn is the interfacial tension. The simplest shape formed by immiscible fluids is the axisymmetric wetting bridge connecting two spheres (edges of the triangulation), for which Laplace-Young equation leads to a one dimensional differential equation tractable with conventional methods [10]. Incorporating pendular bridges in DEM boils down to calculating the traction force resulting from a particular state of the bridge where u n is the distance between spheres of radii R 1 and R 2 . The bridges exist only in a bounded region D of the parameter space. Alternatively, many empirical relations can be found [11], determining f cap , volume V, and the bounds of D. Beyond pendular bridges, the interfaces in a porous medium frequently move in sudden piston-like steps (the so-called Haines jumps) from one equilibrium state to another. The typical equilibrium states have interfaces in the form of meniscii near the pore throats (Fig.4). A meniscus is stable as long as  where p c is the local capillary pressure and p e c is the throat entry capillary pressure. p e c may be approximated by the inscribed sphere method [12] or the Mayer-Stowe-Princen (MSP) method [13,14] (Fig. 1.b), or by direct simulations as in the last section. The MSP method has been found to predict rather accurately measured drainage curves (see Fig. 6, based on [14]). The imbibition events are more complex. Piston-type invasion by the w−phase is possible depending on a new threshold value p f c (the previous drainage event reversed). The interface stability inside a pore needs p c ≥ p f c .
Generally p f c < p e c , which is a major cause of hysteretic responses in multiphase systems. Alternatively, imbibition can proceed by progressive coalescence of pendular bridges, i.e. by first filling edges and facets instead of pore bodies; possibly leading to snap-off : the disconnection of a bubble of n-phase trapped in a pore after all the edges/faces have coalesced. Local rules driving such complex events can be found in [16]. Note that the accurate geometry of interfaces is only known for the pendular bridges [17] (unless direct simulation of the microstructure is performed -see last section). Due to this limitation the interfacial area in the funicular regions is only roughly approximated, which is detrimental in view of thermodynamic analysis [18,19].
In addition to the movements of interfaces each phase can flow within a connected domain just like it does in fully saturated cases -this is handled by the single-phase PFV scheme. Transfers may also occur by film flow in the thin layer of water at the surface of the grains [20,21]. A third type of transfer occurs by evaporation/condensation and vapor advection-diffusion [22]. The last two mechanisms tend to equilibrate capillary pressure all over the pore space after enough time, even between distant, apparently disconnected, clusters. Depending on flow rate and viscosity of the fluids, interface movements can be coupled to viscous dissipation (that is, when the capillary number C is too small for the flow to be quasi-static), leading to secondary pore imbibition in drainage as well as various patterns of fingering [23]. Beyond the pendular regime, very few authors derived analytical expressions of capillary forces in partial saturation [22,24]. A general expression of the force is where p is the pressure of the fluid phase present at point x and δ L (x) is a line Dirac delta reflecting the triple wnsline. Obviously the evaluation of this integral needs accurate knowledge of the geometry of the phases around the particle, which is usually lacking.

Multiscale coupling
Nearly every equation in the previous section can be calibrated with the help of direct fluid dynamics (CFD) simulations of the actual microstructure (Fig. 5). For this purpose we employ a multicomponent Shan & Chen [25] Lattice Boltzmann method (LBM) implemented in Palabos [26]. Here primary drainage is shown for a proof of concept (the reader is referred to the oral or poster presentation for more examples). Primary drainage of an initially saturated material depends, mainly, on the distribution of the entry capillary pressure p e c (Eq. 6) in the pore network. The elementary problem that is solved with the LBM to determine p e c is shown in Fig.7 (note the pendular bridges staying behind the advancing main meniscus). The geometry of this throat-scale problem is defined by the positions of three spheres taken from a large (40,000 sphere) packing. The number of LBM problems to be solved is virtually the same as the number of facets in the triangulation (approximately 85,000). A dynamic selection of the potentially invaded throats is done at runtime to eliminate some of them. The task fits very well in massively parallel architectures where many throats can be simulated simultaneously. After determining p e c for every relevant throat the master (pore-network) process assign the p e c values to the network. In turns, invasion by the n-phase in the boundary value problem is simulated at the network scale. It couples capillary effects (Eq. 6 with updated p e c ) and viscous flow of the receding w−phase (Eq. 1).
This multiscale coupling may seem to be merely one very large LBM problem decomposed in a series a subproblems with nearly the same number of LBM nodes in total -more or less what blind parallelization would achieve. However we argue that the computation time is decreased by order of magnitudes practically. Besides the fact not every throats have to be analyzed, this is mainly because the flow of the receding phase is solved as a porenetwork problem. In fully resolved explicit methods like the LBM, the long range interactions associated to incompressible flow within large w-clusters introduce additional constrained and cause the dominant part of the increased computational cost. In contrast, multiscale coupling enable flexible combinations where the Pore-Network and CFD methods are used for what they do best.