Coupled flow and deformations in granular systems beyond the pendular regime

A pore-scale numerical model is proposed for simulating the quasi-static primary drainage and the hydro-mechanical couplings in multiphase granular systems. The solid skeleton is idealized to a dense random packing of polydisperse spheres by DEM. The fluids (nonwetting and wetting phases) space is decomposed to a network of tetrahedral pores based on the Regular Triangulation method. The local drainage rules and invasion logic are defined. The fluid forces acting on solid grains are formulated. The model can simulate the hydraulic evolution from a fully saturated state to a low level of saturation but beyond the pendular regime. The features of wetting phase entrapments and capillary fingering can also be reproduced. Finally, a primary drainage test is performed on a 40,000 spheres of sample. The water retention curve is obtained. The solid skeleton first shrinks then swells.


Introduction
The situation of two immiscible fluids -say, water and air -going through a deformable granular material is widely encountered in nature and in many areas of engineering and science.Modeling such multi-phase systems has long been restricted to the so called pendular regime, in which the wetting phase (i.e., water in air/water systems) is only present in the form of bridges connecting grain pairs [1].The volume content of wetting phase is constrained to a very low level of saturation (usually less than 5%).The hydro-mechanical couplings for highly saturated states, capillary regime and funicular regime, for instance, have been rarely addressed.
In this work, we present a pore-scale numerical model to simulate the quasi-static primary drainage and the hydro-mechanical couplings, where the wetting phase is drained from a fully saturated state to a residual state but beyond the pendular regime.The basic idea is to combine the discrete element method (DEM) and a pore-network model for, respectively, the granular solid phase and two fluids.The drainage criterion and procedure can be established based on the decomposition of the pore space.The capillary effects, i.e., the fluids forces on grains, can be formulated as well.Then motion of the grains can be predicated by the DEM solver and consequently a macroscale deformation of the skeleton is obtained.In the end of this paper, we will apply the model to simulate the primary of a random packing in order to obtain the water retention curve and the mechanical response.

Network
We consider that the solid skeleton is represented by a dense random cuboid packing of poly-disperse spheres, which is created by DEM [2].The void space is decomposed into large pores connected by narrow throats using Regular Triangulation (RT) method.RT method generalizes the Classical Delaunay Triangulation (CDT) to weight points, where the weight accounts for the size of each solid sphere [3].The difference between RT method and CDT method is shown in Fig. 1.The merit of RT method is that its dual Voronoi graph is entirely contained in the void space, and therefore is appropriately used to describe the flow path.For the boundaries, we represent the rigid confining walls of the cuboid sample by fictitious spheres with infinite radii (i.e., assuming R → +∞, P 1 , P 1 , P 2 and P 2 become aligned) in order to keep the consistency of the geometrical algorithm.
A typical network generated by RT method is shown in Fig. 2 (a).Based on the decomposition, we define the void volume in a tetrahedron as a pore body and its four related facets as four pore throats (see Fig. 2 (b)).The throat does not enclose any volume, but it will play a key role when defining the criterion of the drainage.

Assumptions
We propose the model to simulate the quasi-static hydromechanical behaviors of multiphase granular materials.The system contains of a pair of immiscible wetting and non-wetting fluid phases and a solid phase, which are  noted W-phase, NW-phase and S-phase, respectively.The S-phase is perfectly wetting by the W-phase.The "quasistatic" is restricted by defining the capillary number C a approaching 0. (C a = μv γ , where μ is the viscosity of receding phase, i.e., W-phase; v is the average velocity of the receding phase; and γ is the interfacial tension between the two fluids.)In such case, the viscous effects are negligible.
We consider the state of a single pore is completely either drained or undrained.The saturation is simply binary (0 or 1 depending on which phase is present).Although practically a pendular bridge is formed between two neighboring solid grains, we neglect its volume and capillary effects currently.
For the solid phase, we assume that the motion (mainly sliding) of the rigid solid grains induced by the capillary forces is reversible.Consequently, the macro-scale deformation of the skeleton is of small amplitude and elastic.We assume such deformation will not cause large changes on the RT network connectivity and the fluid phases pressure.To conclude, we present a one-way coupling.

Fluid phases
Locally, the displacement of W-NW interface is controlled by the relation between local capillary pressure p c (the pressure difference between the NW-phase and the Wphase: p c = p n −p w ) and a local threshold value of pressure p c e (termed "entry capillary pressure").If p c > p c e , the lo- ( But C is difficult to define for an interface near a pore throat of complex geometry when p c approaching p c e .Approximations are necessary. We propose to calculate p c e by following MS-P (Mayer-Stowe-Princen) method, which employs the balance of forces on the NW-W interface ( [4,5]).The balance can be written as (see Fig. 3), where, F c is the capillary force acting on pore throat section domain (i.e., A e f f ); and F t is the total tension force along contact lines (i.e., L e f f ).p c e is the value of p c such that F(p c ) = 0.It can be numerically solved.More details for determining p c e can be found in [6].Overall, the drainage follows the invasion percolation model by [7] and [8].Initially, the pores are fully saturated by W-phase.The sample connects to NW and W reservoirs by two opposite boundaries.Drainage starts by increasing p c .The NW-phase invades through the pore throats which separate the phases connected to the two reservoirs.
The first displacement of the W-NW interface occurs through the throat with lowest p c e , leading to the evacuation of the first pore.As soon as this pore is drained the NW-phase reaches new throats, possibly triggering a recursive cascade of displacements and invading more than one pore for one single value of p c , until no more throats satisfy p c e < p c .Then a state of equilibrium is obtained and ready for the invasion of next p c .As the NW-phase is invading, the W-phase may form clusters of pores which are disconnected from the W reservoir.We assume that there are no flux exchanges between individual separated clusters, thus the disconnected regions will remain saturated by a fixed amount of the W-phase throughout subsequent drainage.

Solid phase
During the drainage, the motion of the fluids and interface will cause fluid forces acting on the solid grains, leading to their movement (mainly sliding motion between contact grains).Consequently, a macro-scale deformation is formed.Considering a single solid grain k, the fluids forces can be written by, where ∂Φ k denotes the contact surface between NW-phase and particle k, ∂Θ k denotes the contact surface between Wphase and particle k, and ∂Γ k denotes the W-NW-S contact lines.
Based the RT topology, the domains of ∂Φ k , ∂Θ k and ∂Γ k can be the decomposed and expressed by the geometrical terms of the incident tetrahedral (see Fig. 4).Therefore, Eq.3 can be written in terms of individual fluid pore attributions, where k j represents the tetrahedron j incidental with particle k.F f k j and F t k j represent the corresponding forces induced by fluid pressure and surface tension incidental with tetrahedron j.The algorithm follows [9].After solving the total force on each particle, the DEM solver is employed to predicate its motion.Currently, we are using the traditional Cundall's linear elastic-plastic law [10] to define the contact relations.

Implementation
The model has been implemented in Yade platform, where it is called 2PFV (two-phase pore-scale finite volumes) Weight(%) Diameter(mm) 30 1-1.4 35 0.850 35 0.600 package.The RT network generation has been implemented in C++ [11] using the CGAL library [12].A set of functions has been implemented for solving the hydraulicmechanical couplings (eq.2 and eq.4).As a whole, the model is freely available as part of the open-source discrete element code Yade-DEM.
3 Simulation and results

Numerical setup
We apply the model to simulate the primary drainage of a cuboid packing in order to obtain the water retention curve (WRC) and to observe the mechanical evolution of the solid skeleton.The sample contains 40,000 particles with random positioning, which is generated by using the DEM software Yade [2].The particle size distribution (PSD) is shown in Table .1 and the porosity is 0.34 (replicating the parameters of [6]).The simulation is performed the oedometer test conditons, where only the one-dimensional deformation induced by the drainage can be obtained.The rigid confining walls are set to prevent the lateral displacement of the skeleton.The initial pressure state of the fluid phases is: p c = p n − p w = 0.The sample is fully saturated.Drainage is controlled by decreasing p w and keeping p n constant.The results are put in dimensionless forms.p c is represented by: p c = p c D γ nw , where p c is termed normalized capillary pressure, and D is the average sphere size.

Results and discussion
We observe the drainage process by cutting a slice of the sample, as shown in Fig. 5.The top boundary connects to NW-reservoir and the bottom boundary connects to Wreservoir.The invasion starts from the pores with larger throat, where the entry capillary pressure is smaller.In the beginning, the fluid phase distribution show strong gradients of saturation.The capillary fingering can be observed.After the NW-phase penetrates the sample, the distribution tends to homogeneous.At the end, the remaining W-phase is in the form of disconnected clusters in short-term entrapped by the NW-phase.
The results of the normalized capillary pressure -Wphase saturation -one dimensional strain (p c − s w − 11 ) relationship are reported in Fig. 6.From the WRC (p c − s w ), we can find that the sample is fully saturated when p c < 7.0; the NW-phase invades the pores gradually during 7.0 < p c < 11.0, leading to a small amplitude change of   W-phase content; then s w drops quickly from 0.9 to 0.4 with narrow range of increase of p c , which can be explained by Haines jumps going through large cluster of pores; finally a residual saturation of 0.2 is achieved due to the entrapments.
From the 11 − p c curve, it can be found that the relation is strictly linear when p c < 7.0 (i.e., the sample is fully saturated).It can be explained by Terzaghi's effective stress principle for saturated conditions [13]: σ = σ− p w I, where I the identity tensor.The decrease of p w leads to the increase of effective stress σ .Then the sample keeps shrinking until reaching the peak during 7.0 < p c < 11.0, but the relation loses linearity due to the invasion of NWphase.
With the receding of W-phase, the sample becomes dry.The effective stress for dry scenario is: σ = σ − p n I, where p n = 0 and keeps constant.It is the same as in the initial state (i.e., p n = p w ).Thus, the sample starts swelling.Since the residual W-phase is entrapped in the pore, the deformation of skeleton induced by the capillary effects will not be completely eliminated even the p c keeps increasing.

Conclusions
A pore-scale numerical model for the simulation of primary drainage and the hydro-mechanical couplings has been developed.The underlying idea is to build an explicit and efficient link between fluid space and solid skeleton by a network topology.For the fluids, the local invasion criterion and global logic of drainage has been defined.For the solid phase, the capillary forces acting on grains has been formulated.The model has been implemented in open-source code Yade-DEM.A drainage tests has been performed under the oedometer conditions.The water retention curve has been obtained.The p c − s w − 11 curve illustrated the evolution of the deformation during drainage, in which the solid skeleton showed first shrinking and then swelling.

Figure 1 .
Figure 1.Comparison of 2D triangulation and dual Voronoi graph obtained from classical Delaunay and regular Delaunay: classical Voronoi graph has branches inside discs, while regular Voronoi gives all branches in the void space.

Figure 2 .
Figure 2. (a) Definition of pore network for packing of polydisperse spheres, generated by regular triangulation.(b) A pore geometry defined by tetrahedral element of the finite volume decomposition.

Figure 3 .
Figure3.Geometry of a meniscus.r c is the curvature of meniscus; L nw is the length of contact line between NW and W phases; L ns is the length of contact line between NW and S phases.L e f f = (L nw + L ns ), in which L e f f is termed effective entry pore throat perimeter.A e f f is the planar project area of NW-W interface, and is termed effective entry pore throat area.

Figure 4 .
Figure 4. (a) Decomposition of the fluid space based on regular triangulation method.Particle k is incident with a NW pore and a W pore; (b) Fluids pressure and interfacial acting on particle k. (in 2D for clarity).

Figure 5 .
Figure 5.The process of drainage (40,000 particles), NW-phase invade from top. Brown (gray) is solid phase, blue (black) is Wphase, and light cyan (white) is NW-phase, see colour version of this figure in the HTML.

Figure 6 .
Figure 6.The relationship between normalized capillary pressure-saturation-strain (p c − s w − 11 ) during drainage simulation under the oedometer condition.

Table 1 .
Particle size distribution of porous medium