Multiscale modeling of transport of grains through granular assemblies

We investigate the transport of moderately large passive particles through granular assemblies caused by seeping flows. This process can only be described by highly nonlinear continuum models, since the local permeability, the advection and dispersion mechanisms are strongly determined by the concentration of transported particles. Particles may sometimes get temporally trapped and thus proper kinetic mass transfer models are required. The mass transfer depends on the complexity of the porous medium, the kind of interaction forces and the concentration of transported particles. We study these two issues by means of numerical and laboratory experiments. In the laboratory we use an oedo-permeameter to force sand grains to move through a gravel bed under conditions of constant hydraulic pressure drop. To understand the process, numerical experiments were performed to approach particle transport at the grain scale with a fully coupled method. The DEM-PFV combines the discrete element method with a pore scale finite volume formulation to solve the interstitial fluid flow and particle transport problems. These experiments help us to set up a continuum transport model that can be used in a boundary value problem.


Introduction
Particle transport in porous media is a widely studied problem with applications in many different contexts [1]. Even when the macroscopic conditions for the flow and transport problem are homogeneous, preferential paths for particle and fluid flows as well as localized filtration regions are often observed [2,6]. We investigate how seeping flows convey passive particles through granular assemblies made of coarser particles (size ratio, transported to coarse particles, around 1/10). In such circumstances there is a significant relative motion between the transported particles and the carrier fluid, as shown by their average velocities or tortuosities (the length of a path connecting two points divided by the distance between them).

Continuum transport model
Under a continuum approach [3], the system can be viewed as a three-phase medium with a representative volume element V consisting of (i) a solid fraction of coarse and fixed particles (the coarse skeleton) V c , (ii) a set of solid fine particles not attached to the coarse skeleton V p and (iii) a fluid that fills the space left by solid particles V w . The porosity of the clean coarse medium (i.e. without fine particles in the voids) is Φ 0 = (V − V c ) /V while the actual porosity is Φ = V w /V. The concentration of transported particles is defined as In the case of moderately large particles, the particle flux, e-mail: ignacio.gtejada@upm.es q p , is much smaller than the water flux, q w = Φv w (where v w is the velocity of water), so that the flow problem can be separately solved: Darcy's law states q w = 1 μ κγ w ∇ (h), where μ the viscosity of the fluid, κ the intrinsic permeability and h = p/γ w + z the hydraulic head (with z the elevation, γ w the specific weight of water and p the fluid pressure). Some theoretical approaches aim at explaining the influence of the features of the coarse skeleton and the concentration on κ. For example, in the Kozeny-Carman's [4,5] model, when the porous medium consists of a collection of curving pipes through which the water flows according to Poiseuille's law. The average tortuosity of the pipes is T . We refer to κ 0 as the permeability of the clean structure, which can be expressed in terms of the average tortuosity of the clean structure, T 0 . The inclusion of fine particles changes T and Φ and makes κ dependent on C.
Regarding the transport of passive particles, the advection -(normal) dispersion paradigm states that the evolution of C is given by: where v p is the advective velocity of the particles and D is the hydrodynamic dispersion tensor. The advective and dispersive fluxes of particles are, respectively q p,a = v p C and q p,d = − D∇C (Fick's law). The advective velocity is parallel to that of the fluid. The relative advective velocity of the particles is defined as λ = v p /v w and depends on the coarse skeleton, the transported particles and C (as observed in interparticle percolation by gravity [9,10]). The hydrodynamic dispersion depends on v w (e.g., [7]) and can be expressed as D = |v w |δ, where δ is the dispersivity, to emphasize the fact that the dispersion is scaled by the water velocity. Eq. 3 is valid when the number of transported particles is constant (they do not join the coarse skeleton and no new particles are eroded from the later). However, sometimes particles can be temporarily immobilized. For example, when they get trapped due to sticking at surfaces, crevices, caverns, etc. [8], or when they clog constrictions (individually or collectively). Then, kinetic mass transfer laws must be included in the right-hand side of Eq. 3. These laws depend on the complexity of the porous medium and interaction forces, but also often do on C. For example, in this case, the intermittent, purely mechanical, formation of blockages in the constrictions due to physical crowding can be properly reproduced by 1 st -order kinetic mass transfer laws, and Eq. 3 becomes: and C f(x,t) are, respectively, the concentrations of trapped and free particles (C (x,t) = C t(x,t) + C f(x,t) ), τ t is the mean lifetime of trapping events and τ f is that of travelling periods. These mean lifetimes may change considerably with C. As κ (and hence v w ), λ, τ f and τ s , and, to a lower extent, δ, are strongly determined by C, the transport of moderately large particles can only be described by a highly nonlinear continuum model. A micromechanical approach may help to get some insight on the effect of C on these parameters.

Micro-mechanical transport model
Seeping water conveys fine particles through the coarse skeleton by means of viscous and hydraulic pressure forces. The particles interact mechanically with each other and with the coarse skeleton. Under a micromechanical approach, their individual paths are considered.
From grain level considerations, Eqs. 1 and 4 can be derived. The computation of the q w over a representative volume element is straightforward, so that if the hydraulic pressure drop is known, then κ is known for each given value of C. Regarding q p , individual particles can be tracked to generate an ensemble of trajectories over large times and distances. The center of mass of the particles is expected to move linearly with time (t) at the advective velocity: The Fick's law requires a linear growth with time of the centered mean square displacement of transported par- However on very short time scales (e.g. below the typical time for interpore transit) the motion of particles involves time correlations and Fick's law does not applies. On timescales comparable to the lifetime of blockages or to the typical time between consecutive blockages [11], the total C can be splitted into trapped C t and free fractions C f and kinetic mass transfer laws can be used (as in Eq. 4). 1 st -order kinetic mass transfer laws are justified by exponential statistical distributions of the duration and frequency of trapping events [11]. For a given C the probability of a particle to get trapped depends on the number of constrictions that it crosses per unit of time, so that τ f should decrease almost linearly with the particle velocity.
In other words, τ f |v p | = l f , where l f is the filtration length.
In contrast, the duration of trapping periods, represented by τ t , should be related to the hydraulic pressure gradient, so that τ t = f (C). As a first approximation, τ t |v p | = l t , where l t is the length that the particle would have moved if it had not been trapped.

Methodology
We performed a laboratory experiment and many simulations of particle transport. The laboratory experiment allowed us to understand the transport of moderately large particles caused by water seepage whereas the numerical simulations helped us to set up a continuum transport model. At this stage of the research the simulations did not aim at reproducing the laboratory experiments. The experimental device used at the laboratory was a large oedo-permeameter (see a detailed description in [12]). This device comprises of a cylindrical cell made of a rigid transparent acrylic tube of diameter 0.279 m. The granular specimen, depicted in Fig. 1, consisted of a binary mixture of sand and gravel (1:1) placed over a 0.200 m high bed of gravel. We used Fontainebleau sand (median diameter, D 50 = 0.207 · 10 −3 m, i.e. 50 % of the grains, by mass, is finer than 0.207 · 10 −3 m) and a gravel from Sablière Palvadeau (D 50 = 3.309 · 10 −3 m). The binary mixture was prepared with a water content of 3 % (by mass) and placed to fill the gap between the acrylic tube and an impervious internal plastic bucket with an average diameter of 0.209 m. The bucket left a 0.150 m high ring (see Fig. 1). This configuration allowed us to observe how the sand spread in the gravel bed as the water flowed. The whole sample was saturated first with CO 2 and then with water and subjected to a downward water seepage caused by a fixed hydraulic pressure drop. After some time, the water flow was stopped and the cell was gently drained. Some samples were collected at different sites and their sand content was measured (the ratio of the dry weight of sand to the dry weight of the sample). The fact that the specimen remained partially saturated, helped to keep sand particles where they were at the end of the seepage stage.
The numerical method we use is the DEM-PFV, which applies for incompressible and Newtonian fluids at low Reynolds numbers. The discrete element method (DEM) [16] is combined with a pore-scale finite volume formulation (PFV) [13][14][15]. We use linear elastic contact laws between the frictionless particles. The finite volumes are geometrically defined after a tessellation of the granular assembly 1 : finite volumes are the tetrahedra whose apexes coincide with particle centers. When solid particles move, tetrahedra change, and the variations of volume must be equal to the net water fluxes over their four faces (mass balance and incompressibility conditions). Interpore fluxes are supposed to be proportional to the difference of hydraulic pressure between consecutive tetrahedra, the distance of their centers of mass and the hydraulic section of the shared face. This is congruent with Stokes flows [13]. The forces that the fluid flow exerts on solid particles are computed after the integration of fluid pressure and viscous stresses over the wet surfaces of solid particles and the hydraulic sections of tetrahedron faces. This ensures a full coupling.
In the simulations we used random, crystalline and quasi-crystalline 2 fixed granular assemblies of monodispersed spheres of size d c = 5 · 10 −3 m. These structures are 1 Transported particles are also considered for the tessellation. 2 Randomly generated from the corresponding crystalline structures by enlarging unit cells and slightly displacing particles from crystalline positions. included in parallelepiped periodic simulation boxes. Fine particles of size d f = d c /8 are randomly introduced within the voids (all over the simulation box) up to an average concentration C . The size ratio ensures that fine particles cannot get trapped alone but allows them to collectively form clogging arches in the constrictions that may cause larger blockages. For any coarse skeleton, blockages are more likely at high C. The simulations provide us a statistical sample of particle trajectories. The longitudinal velocity, the dispersion (in the direction of −∇h), τ f and τ t can be measured from the statistical sample.

Results and discussion
In the laboratory experiment, the ring containing the mixture was eroded in a few seconds under the action of the water seepage. The process was localized. The main water flow occurred through a way opened in the ring. The sand particles that had been eroded from the ring moved down and settled. No significant amount of sand was collected at the outlet of the cell after the experiment and the sand content at different sites within the specimen is shown in Fig. 2. Although the corresponding boundary value problem was supposed to be axisymmetric, the transport was not symmetric at all. This instability unveils a highly nonlinear behavior: the water flows where there is less sand and the particle transport is stopped when water finds continuous paths to flow and deposits the sand aside. This emphasizes the importance of the relationship between κ and C.
The numerical experiments show (Fig. 3a) how κ changes with C. These values have been fitted with Kozeny-Carman models (Eq. 2). For each coarse skeleton, an expression of the form T = a √ C, for C > C 0 , and T = T 0 otherwise, was used. Each coarse skeleton had a porosity and tortuosity that determined its permeability. In all the cases, v p decreased with C, whereas λ increased with C (Fig. 3b). This is because with a higher number of fine particles makes them to move more evenly with the fluid. D is scaled by v w , but δ remains considerably constant at different C levels (Fig. 3c) (i.e. the dispersivity mainly varies from one coarse skeleton to another). Finally, with respect to the formation of collective clogging arches and blockages in the constrictions (the only trapping mechanism that our model is able to reproduce; only for the FCC structure), at low C, the expected walking length can be very long while it decreases greatly with C (Fig. 3e). The lifetime of trapping events is less variable (Fig. 3d).

Concluding remarks
Laboratory and numerical micromechanical experiments show that the transport of moderately large passive particles caused by seeping flows requires nonlinear transport models which may lead to preferential paths, where water flows and particles are deposited aside. In these models the macroscopic permeability is strongly determined by the presence of particles. In consequence, the concentration scales advection, dispersion, filtration and erosion processes. Furthermore, the relative advective velocity of the particles, the filtration and erosion times and the dispersivity also depend on the concentration. The next challenge is to embed such a macroscopic model, including all these dependencies, in a continuum numerical method to tackle boundary value problems.