Pore-Scale Flow Simulations: Model Predictions Compared with Experiments on Bi-Dispersed Granular Assemblies

A method is presented for the simulation of pore flow in granular materials. The numerical model uses a combination of the discrete element method for the solid phase and a novel finite volume formulation for the fluid phase. The solid is modeled as an assembly of spherical particles, where contact interactions are governed by elasto-plastic relations. Incompressible Stokes flow is considered, assuming that inertial forces are small in comparison with viscous forces. Pore geometry and pore connections are defined locally through regular triangulation of spheres, from which a tetrahedral mesh arises. The definition of pore-scale hydraulic conductivities is a key aspect of this model. In this sense, the model is similar to a pore-network model. Permeability measurements on bi-dispersed glass beads are reported and compared with model predictions, validating the definition of local conductivities.


INTRODUCTION
Microscale measurements and modeling of porous materials are active areas of research in geomechanics and geosciences. They aim at a better understanding of some of the complex phenomena taking place in porous materials, in connection with flow of pore fluids (Mohammadzadeh and Chatzis, 2010), mechanical behavior (Gueguen and Bouteca, 1999), thermo-chemical processes and couplings between these effects (e.g. Renard et al., 2005;De Boever et al., 2012). Such phenomena are strongly linked with some of the challenging problems of energy production.
For modeling the solid phase alone, the Discrete Element Method (DEM) has been developed in the last few decades (Cundall and Strack, 1979). It is now commonly employed for studying aspects of mechanical behavior on the microscale for granular materials (Scholtès et al., 2009) and rocks (Scholtès and Donzé, 2012). It is based on the description of materials as a collection of individual particles interacting at contact points. The DEM has had great successes in the study of stress-strain behavior and fracturing in dry materials. The simulation of problems including a coupling with one or more interstitial fluids with the DEM is also very attractive and it has many potential applications but it is much more challenging. Various strategies have been adopted for such couplings. They essentially differ in the modeling techniques adopted for the fluid part of the problem. They can be classified into three groups: -Microscale models are based on a fine discretization of the void space. Finite Element Methods (FEMs) are often used because of their flexibility in the definition of the numerical mesh (Glowinski et al., 2001) but they involve a large computer memory footprint and long computational times. Lattice-Boltzmann (LB) (McNamara et al., 2000, Mansouri et al., 2011 does not resort to the solution of large systems of nonlinear equations. In general, LB is faster than the FEM even if commonly implemented fixed-size grids in 3D can result in considerably larger computer memory occupancy (Ma et al., 2010); -Continuum-based models define flow and solid-fluid interactions on the mesoscale using empirical relations (e.g. Darcy's law), generally resulting in acceptable computational costs (Nakasa et al., 1999). There is no direct coupling on the local scale: the forces acting on the individual particles are defined as a function of a mesoscaleaveraged fluid velocity obtained from porosity-based estimates of the permeability. Therefore, the individual particle behavior is not accurately reproduced and this limits the application of the model to problems such as strain localization, segregating phenomena, effects of local heterogeneities in porosity and internal erosion by transport of fines, which are all inherently heterogeneous on the microscale; -Pore-network modeling has been commonly developed to predict the permeability of materials but has also been extended to multiphase flow (Bryant and Blunt, 1992). It relies on the representation of the pore space as a network (Jiang et al., 2007). Crucial for its success is an adequate definition of how fluids are exchanged between pores with respect to the local pore geometry. In the context of a coupling with DEM models, the accurate definition of forces applied to individual entities constituting the solid phase is a key point. This is only rarely addressed in porenetwork techniques, which tend to disregard the deformability of the solid matrix (definitions of forces may be found, however, in Jing et al. (2001) for polyhedral solids). Early ideas on coupling DEM models of granular materials and pore-network methods can be found in the works of Hakuno (1995) and later by Bonilla (2004) but these studies were limited to 2D models of disc assemblies. We recently adapted this approach to 3D sphere assemblies using a Pore-scale Finite Volume (PFV) scheme in Chareyre et al. (2012). The three-dimensional aspect of the problem let us define local hydraulic conductivity as a function of actual pore-space geometry, which is described as a network of connected pores (Fig. 1a, b). This in turn opens up the possibility of obtaining macroscale permeability as a result.
In the first section of this study, we give the principles of the DEM and we summarize the governing equations of the DEM-PFV coupling. Permeability measurements on mixtures of glass beads of different sizes are reported in the second section. The results are compared with the predictions of the PFV model and with existing empirical or semi-empirical relations.

Discrete Element Method
The Discrete Element Method (DEM) has been extensively used to study soil and rock mechanics providing, for instance, some insights into shear strength and deformation properties of geomaterials. The DEM is essentially a Lagrangian (meshfree) technique where each particle of the material is a sphere identified by its own mass, radius and moment of inertia. We employ here the DEM in 3D as implemented in the open source platform Yade-DEM. In what follows, we only briefly recall some basic aspects of the method (for more details see Šmilauer et al., 2010). For every time step of the computation, interaction forces between particles and consequently the resulting forces acting on each of them, are deduced from spheres' positions through the interaction law. Newton's second law is then integrated through an explicit second-order finite difference scheme to compute the new positions of spheres.
The linear elastic-plastic model for contact interactions between solid particles is presented in Figure 2a. A fundamental assumption of the DEM is that the deformation of the solid particles is concentrated in a very small zone in the  vicinity of the contact area. This deformation corresponds to an apparent overlap, as shown in Figure 2b, actually reflecting the relative movement of centers of mass along the normal of the contact point. A linear elastic law provides the contact force as a function of this relative displacement between two interacting grains (see Fig. 2c, d).
The normal stiffness K n defines the normal force F n for a given normal displacement U n as: (1) and the tangential stiffness K t defines the increment in the shear force vector F t induced by an increment in tangential relative displacement dU t when in contact: In our simulations, K n and K t are dependent functions of the particle radii R 1 and R 2 and of a characteristic elastic modulus E of the material such as: ( 3) with α being a fixed parameter. This definition results in a constant ratio between E and the effective bulk modulus of the packing, whatever the size of the particles.
Shear and normal forces are finally related by a slip Coulomb model such that F t max = tan(φ).F n , where φ is the contact friction angle.
Although the linear formulation employed here is quite simplistic compared with more realistic Hertz-Mindlin solutions, it has been proven to give accurate results (DiRenzo and DiMaio, 2004), with significant computational efficiency. However, it should be noted that the fluid model presented in the next sections is independent of the contact law, hence not limiting the potential usage of different contact laws in future studies.

The DEM-PFV Coupling Method
We give here the main features of the Pore-scale Finite Volume (PFV) scheme that we developed for the coupling with DEM. For a more detailed derivation, the reader can refer to Chareyre et al. (2012).
A partition of the pore space in a sphere assembly can be obtained by constructing the regular triangulation of the packing, as allowed by the open-source library CGAL (Pion and Teillaud, 2011). The elementary geometrical objects emerging from this procedure are tetrahedra with spheres at each vertex (Fig. 1c), and enclosing what we call hereafter voids. Then, it is possible to derive a finite volume formulation for Stokes flow in such a mesh. The mass balance equation will link the rate of change of an element's volume to the fluxes q ij exchanged between adjacent pores through the facet S ij : (4) where u n and v n are, respectively, the normal component of fluid velocity on the contour of the element and the velocity of the contour. Since the geometry of the mesh follows the motion of the solid particles, (u n -v n ) is the velocity of the fluid relative to the velocity of the solid. Hence, this equation couples the deformation of the solid skeleton and the fluid flow.
The flux q ij can be expressed as a function of the pressure jump between elements i and j. We proposed to define the function by using the hydraulic radius of the throat between the two elements (Chareyre et al., 2012), the throat itself being defined by the branch that is the dual of the facet S ij in the Voronoi graph (Fig. 1d, e). This finally links the deformation to the pressure field {p k }, that we assume to be piecewise constant. If the conductivity of the throat is noted k ij , Equation (4) becomes: This equation gives the linear system that has to be solved at each time step of a simulation, giving the pressure field as a function of the velocity field of the particles. The forces exerted by the fluid on each particle can finally be derived from the pressure field, with the help of a momentum conservation equation. It can be decomposed into three terms which are contour integrals of the hydrostatic pressure ρgz (Archimede's force), of the piezometric pressure p* and of the viscous shear stress τ, respectively: (6) with n the unit normal of the contour. These forces can be easily integrated into the conventional time-stepping algorithm of the DEM presented in Section 2.1 by summing them with the contact forces.  Computation cycle of the coupled DEM-PFV model. of the pores is computed at each time step, then the system defined by Equation (5) is solved to obtain the pressure field, and new forces are computed for the next step according to Equation (6), as shown in Figure 3.

PERMEABILITY OF BI-DISPERSED SPHERE PACKINGS
A series of experiments and numerical tests were performed on granular samples in order to evaluate the ability of the PFV model to predict their permeability. The experimental setup is presented below. Some classical relations giving estimates of the permeability as a function of PSD and porosity are also recalled and later compared with the experiments.

Experiments
The experiments were done on samples made of sphericalshaped glass beads of two different sizes. The small particles have diameters ranging from 0.50 to 0.63 mm, while the coarse particles are between 2.80 mm and 3.15 mm. They were classified by a mechanical sieve analysis. The two types are mixed with various mass ratios: M = m 1 / (m 1 + m 2 ), where m 1 and m 2 are, respectively, the mass of the fine and coarse beads. In order to get homogeneous samples, it is important to pour the beads in thin successive layers, one on top of the other, to minimize segregation effects. For each definite mass ratio, three samples were prepared, two of which were used for the estimation of density, porosity and median grain size, whereas the third one was prepared for the permeability test. It must be noted that the different mass ratios resulted in different sample porosities (Fig. 4), which will have to be reproduced in the numerical simulations (Sect. 2.2).
Water (temperature of 20°C) flowing through samples of section S = 0.01 m 2 and length L of 0.20 m under a controlled water head gradient was used for permeability tests. To see how permeability varied with time and to ensure that saturation was achieved, flow rate (Q) measurements were repeated for a constant hydraulic head Δh = 0.05 m. Once saturation was achieved (no variability of Q in time), the permeability test was carried out. A series of water heads (from 0.02 to 0.20 m) were tested to measure the hydraulic conductivity k. Flow rates are estimated from measured volumes at the outlet. Each measurement was repeated at least twice for each head level.
In what follows, we will compare the values of a normalized permeability, defined as, where is a normalization factor that corresponds to the hydraulic radius, V v is the void volume filled with the fluid, and S w is the wetted surface of grains. K* only depends on the shape and connectivity of the pore space.
The results obtained for the different values of m 1 /(m 1 + m 2 ) are shown in Figure 5 and in Figure 6 for the case of m 1 / (m 1 + m 2 ) = 0.2.

Simulation Procedure
A requirement for meaningful comparison with experiments is that the simulated microstructure must be as close as possible to the experiments. In the absence of advanced experimental techniques, it is not possible, however, to determine the position of each individual particle. In this study, it was decided to at least match global quantities: the porosity and the Particle Size Distribution (PSD).
It is generally difficult for purely geometric space-filling algorithms to generate random sphere packings when both porosity and PSD are prescribed. Instead, the numerical samples were obtained by simulating a mechanical problem corresponding to isotropic compaction, as initially proposed in Chareyre et al. (2002). This method combines growth of particles and lubrication of contacts under constant applied stress: -in the first step, a loose cloud (porosity n = 0.7) of nonoverlapping spheres is generated in a cubic box with random positions. The radii are assigned with respect to a Porosity measured after sample preparation and in the reconstructed numerical samples.
PSD that is homothetically equivalent to the prescribed PSD; -in the second step, a steady growth of the particles is simulated. The particles are free to move and tend to push each other while the growth creates contact with repulsive forces, progressively filling the whole volume of the box. The growth goes on until the packing is stable and a certain level of stress is applied to the boundaries. The interparticle contacts created during this growth phase are assigned a high friction value. As a result, the internal reorganization of the microstructure is limited and the porosity obtained at the end is typically higher than the target value; -in the last step, the friction of contact is decreased progressively. This decrease triggers reorganizations of the contact network and compaction of the packing, while the stress is kept constant by further increasing the sizes. This lubrication process is continued until the porosity reaches the value of porosity that was measured in the experiments. Figure 7 shows the evolution of the system in this algorithm. In the end, a stable and isotropic assembly of spheres is obtained. The unbalanced force is represented as an indicator of the stability of the system (Šmilauer et al., 2010); it tends to zero at static equilibrium. In Figure 7, the values are of the order of 10 -2 , which indicates that the evolution is close to the quasi-static limit during the process.
The resulting porosity is homogeneous inside the box on average, and matches the target values very well (Fig. 4). The regions near the boundaries have generally higher porosity, however. The impact of this boundary effect on simulated permeability will be evaluated below.
Note that the PSD obtained after this procedure is only homothetically equivalent to the experimental one. Although it would have been possible to scale the final packing by a size factor in order to match the absolute PSD of the experiment, it is not necessary for comparisons. The numerical model implies a strict proportionality between the permeability and the squared reference length (Chareyre et al., 2012). Therefore, the permeability given by the simulation can be easily extrapolated to the real PSD using the factor of homothety.  Finally, the permeameter tests are simulated on the numerical samples by imposing different values of the fluid pressure at the top and bottom of the box (Fig. 8). No flow is allowed through the lateral boundaries, where a no-slip condition is assumed. The permeability is deduced from the fluxes through the top and bottom boundaries, using Equation (7).
We recall that the current PFV formulation is assuming Stokes flow of pore fluid. It is usually considered that this regime is a valid approximation as long as the Reynolds number R e is smaller than 1. In our permeameter tests, R e is between 0.04 and 0.81, hence it is relevant to compare the model with such experiments.
The number of spheres in the simulations is much smaller than in the experiments, due to the limitations of computing power. In order to evaluate the possible size and/or boundary effects, we simulated permeability tests on packings of N = 5000 to N = 200 000 spheres. The results are plotted in Figure 9. For N = 5 000 to N = 50 000, twenty samples with the same porosity and same PSD were generated and the standard deviation of the permeability is also provided in Figure 9. Given the small dispersion that was found for N ≤ 50 000 (less than 1%), it was decided to simulate only one test for the higher values of N. These results suggest that boundary effects are sensitive up to N = 50 000 (although the fluctuation of the mean value is not extremely high). For higher N, the normalized permeability is stable. N = 50 000 was taken as a sufficiently large sample, to reduce the computational time.

Empirical Relations
The empirical studies of Schlichter and the semi-empirical approaches of Kozeny-Carman, and Terzaghi and Peck (1964) relate K to particle size and shape, PSD, and porosity. They proposed relations that we will compare with the experimental results.
The formula of Schlichter (1905) reads: The following semi-empirical formula was developed by Terzaghi and Peck (1964): (9) where C is the sorting coefficient, ranging from 6.1 × 10 -3 to 10.7 × 10 -3 ; the average value of C was used in this study.
The Kozeny-Carman (K-C) equation, first proposed by Kozeny (1927), then modified by Carman (1956), is known Evolution of porosity, the contact friction coefficient and unbalanced force in the procedure of sample generation, for m 1 / (m 1 + m 2 ) = 0.70.
for accurate estimation of K (Carrier, 2003;Ishaku et al., 2011;Odong, 2007): where S v is the specific surface area of the grains.

RESULTS AND COMPARISON
The normalized values of permeability that were found in simulations, experiments and calculations from the different formulas are summarized in Figure 5 as functions of the mass ratio M = m 1 / (m 1 + m 2 ). The general trend of the simulated and empirical curves resembles that from the laboratory tests. The maximum K* is obtained when the packing is mono-dispersed (i.e. when M is 0 or 1). The permeability is minimized near M = 0.3 in all cases. The general shape of these curves is consistent with the evolution of porosity (Fig. 4) but is not strictly the same. Namely, K * is almost constant in the range 0.7 < M < 1 while the porosity is not. This trend is reflected correctly by the simulation and by the empirical relations.
Quantitatively, it is found that the empirical relations tend to overestimate the measured values for small M. Surprisingly, the K-C equation gives the worst estimate in all cases and overestimates K* by a factor of 5 on average, while the other relations give relatively good estimates when 0.2 ≤ M ≤ 0.9. The K-C equation is commonly accepted for well-graded materials but our results suggest that it is not well suited for bi-dispersed (0 < M < 1) or mono-dispersed (M = 0 or M = 1) packings. The PFV predictions are satisfying for M ≤ 0.5, with errors of the order of 50%. The worst estimates are for 0.7 ≤ M ≤ 0.9, with an error of a factor of 2.5 in the worst case. Overall, the errors in model predictions are of the same order of magnitude as those of the empirical relations. Qualifying these errors as negligible or significant is a difficult task since the impact of such errors will depend on the application to be treated.
We recall that, ultimately, the DEM-PFV coupling is intended to simulate hydro-mechanical couplings. As such, it is not supposed to compete with numerical models more specifically dedicated to permeability predictions. However, a requirement for the model is that it should be computationally efficient and should reflect changes in permeability when the porosity is modified by deformation or by internal erosion. On this aspect, the results are satisfying, since the model follows the trend of the experimental data over a wide range of PSD (0 < M < 1) and porosity (0.28 < n < 0.415).

CONCLUSION
The coupled DEM-PFV method has been presented. Permeability tests on bi-dispersed packings of glass beads have been reported and a procedure to generate numerical samples matching experimental PSD and porosity has been proposed. The predictions of the model are in good agreement with the experimental data, reflecting the role of PSD and porosity in permeability. This conclusion is in agreement with some comparisons between PFV and FEM models, reported previously by Chareyre et al. (2012). The difference between the experimental results and the simulations is of the Influence of number of grains in the simulation on the results of permeability; m 1 / (m 1 + m 2 ) = 0.5. same order as the error obtained with the empirical relations of Terzaghi and Schlicher and it is much smaller than the error obtained with the Kozeny-Carman relation. It is worth noting that the flow model does not restrict the variety of constitutive behavior that can be assumed at contact between particles. Hence, the coupling allows, for instance, the study of fracturing processes in cohesive materials simulated with the DEM in the presence of (or because of) a fluid.
A restriction of the model is due to the discretization of the pore space via regular triangulation, assuming spherical particles. However, it can be noted that a current trend in the field of DEM modeling is to approximate complex shapes of particles by clusters of spheres. In this situation, regular triangulation holds and the method could be easily generalized.
In terms of accuracy of the predicted permeability, the PFV model will probably not compete with other numerical methods which solve the flow problem on a smaller scale. Since natural materials are rarely composed of perfectly spherical particles, the practical interest of a very precise estimate would be limited anyway. This weakness may be balanced by robustness and ability to solve large problems efficiently. It may be the case especially in situations of strong coupling with mechanical effects and/or other physical effects (e.g. thermo-chemical), as found in reservoir engineering.
The method presented here has been successfully applied to transient regimes in poromechanics problems . Internal erosion and transport of fines (as in Sari et al., 2011), as well as hydro-fracturing, are examples of the challenging problems which may benefit from the microscale insights offered by this type of coupling.