Role of blockages in particle transport through homogeneous granular assemblies

The transport of passive particles in porous media is a widely studied problem with applications in many different contexts [1]. Particle transport is determined by the interaction of particles with each other and with porous medium surfaces. This letter focuses on the transport of solid particles of diameter larger than 5 · 10⁻⁵ m which are moderately large in comparison with the constrictions and which interact exclusively via excluded volume repulsion. Other forces, such as those acting between colloidal particles in suspensions, are not considered. These conditions are met, for instance, during the internal erosion of silty or sandy soils. Even when the porous medium and initial boundary conditions (for both the fluid flow and the particle transport) are supposed to be macroscopically homogeneous, the experience has confirmed the presence of preferential fluid and particle flows [2,3]. The reasons behind them not only have to do with local variations in the permeability of the clean porous medium but also with the relocation of transported particles itself. Indeed, in some cases, particles can be temporarily immobilized, adding complexity to their transport. For example, particles may get trapped because they stick at surfaces, crevices, caverns, etc. [4], or because they clog constrictions (individually or collectively).

Anomalous transport processes [5,6] are referred to those cases in which the advection-Fickian dispersion paradigm fails to reproduce the evolution of the concentration of particles in space and time. The complexity of natural porous media [7] may involve multiscale heterogeneities, which are known to be a reason for anomalous transport [8–12], but even in quite homogeneous systems, preferential flow paths may lead to anomalous transport [13,14]. For instance, in random fractured networks subtle features of the steady flow-field distribution through the network are behind time-dependent anomalies [15].

Our objective is to investigate the role played by the intermittent formation of blockages in the constrictions of the porous medium (fig. 1). We want to understand the length and time scales affected by these events and their consequences for particle transport. We focus on purely mechanical clogging: one or several inert particles form an arch in a constriction via mechanical forces that locally block particle flows. The probability of forming an arch in these conditions depends on the size and shape of constrictions and particles, on the friction, and on the number of them passing simultaneously through the constriction. Indeed, sometimes the size of the crowd is known to be more significant than the driving force to explain the ability of a constriction to clog [16], but this ability decreases with the square root of the kinetic energy [17]. Once an arch is formed, other particles can be retained upstream, increasing the blockage. Dynamically changing fuid flows or collisions from other particles can break up the blockage and cause the collapse of the clogging arch. Whether or not these events have implications for particle transport is something strongly related with the statistics of their duration and frequency. Exponential and power laws have been found in other cases of intermittent flow to explain, for example, avalanche sizes during silo discharge [18], burst sizes in escape panic experiments [19], clogging duration in the discharge of granular hoppers under external vibration of low amplitude [20] or time lapses between the passage of consecutive particles through a bottleneck [21]. In our case, if the statistical distributions of waiting and traveling times have finite second moments, on a time scale much longer than the time scale of blockages, the classical advection-dispersion paradigm will be valid to explain the evolution of the concentration of particles. Similarly, on a time scale comparable to that of blockages, the distinction between retained and transported particles together with appropriate mass transfer laws, is required to properly model particle transport. These laws have been often proposed either phenomenologically or after micromechanical approaches (e.g., [4,22–24] for colloids). In contrast to these approaches, we deduce them for the specific case of mechanical clogging due to physical crowding, after a statistical analysis of the frequency and duration of blockages.

Zoom In Zoom Out Reset image size

In particular, we study the transport of passive particles in macroscopically homogeneous media. We consider unchangeable granular assemblies and viscous and incompressible fluids that flow through the assemblies due to stationary hydraulic pressure gradients $\boldsymbol{\nabla}p$. Gravity is not included. The fluid flow can be described at the scale of the interstitial space by the Stokes equations (valid at small Reynolds number). For a given flow field the hydraulic tortuosity $\mathcal{T}$ accounts for the effective elongation of fluid paths. Equivalently, on the scale of the granular sample, Darcy's law can be applied to obtain the average fluid flux $\bm{q}_\text{w}$ (discharge per unit area) as $\bm{q}_\text{w} = -\frac{\bm{\kappa}}{\mu} \boldsymbol{\nabla}p$, where μ is the dynamic viscosity and κ is the permeability tensor. κ depends on the pores structure, the porosity ϕ (the fraction of the volume of the voids over the total volume) and the constriction size distribution CSD [25]. The fluid velocity in a clean porous medium of permeability $\bm{\kappa}_\text{0}$ is $\bm{v}_\text{0,w} = \bm{q}_\text{0,w} / \phi_\text{0}$.

The introduction of fine particles in the voids of the structure results in a new granular assembly of longer and thinner tortuous passages of varying cross-sections. The permeability is reduced since $\phi < \phi_\text{0}$ and the constrictions are smaller, too. If fine particles are not attached to the granular assembly, the Stokes flow is able to convey them through the porous medium via pressure and viscous forces. When transported particles are of size comparable to that of pores, the relative motion between them and the carrier fluid is very significant. Moreover, the average tortuosity of particle paths $\mathcal{T}_\text{p}$ is different from $\mathcal{T}$ because fine particles do not follow closely the streamlines. As particles move they change the domain for the fluid flow.

If $C = V_\text{s}/V$ is the concentration of transported particles (being $V_\text{s}$ the volume of the transported particles and V the total volume of an element), the advection-dispersion equation ADE, based on the macroscopic conservation of particles, reads

$$\begin{equation} \frac{\partial C}{\partial t} + \boldsymbol{\nabla}\cdot(\bm{v}_\text{p} C ) - \boldsymbol{\nabla}\cdot\left( \bm{D} \boldsymbol{\nabla} C \right) = \bm{0}, \end{equation} \tag{ 1 }$$

where $\bm{v}_\text{p}$ is the advective velocity of the particles and D is the hydrodynamic dispersion tensor. Although the ADE is a continuum balance law, it can be also derived from grain level considerations. In effect, the advective flux of particles $\bm{j}_\text{a} = \bm{v}_\text{p} C$ accounts for the fact that particles tend to follow, on average, the prevailing directions marked by the fluid flow. Then $\bm{v}_\text{p}$ is parallel to $\bm{v}_\text{w}$ and $v_\text{w}$ is always higher, and often much higher, than $v_\text{p}$. For a given porous medium, $v_\text{p} / v_\text{w}$ depends on the size of the particles with respect to that of pores, their shape, the interparticle friction and the value of C, as has been shown in interparticle percolation by gravity [26–29]. The dispersion mechanism is exhibited in a frame of reference that moves with the center of mass of the cloud of particles. Dispersion is random in nature. Fick's law $\bm{j}_\text{d} = - \bm{D} \boldsymbol{\nabla}{C}$ takes up the idea that dispersion occurs from high- to low-C regions. The dispersion tensor depends on $\bm{v}_\text{w}$ (e.g., as expressed for isotropic media in [30]).

Normal dispersive transport can be conceived as the macroscopic manifestation of an underlying microscopic uncorrelated and Gaussian stochastic process. Normal and anomalous transport can be identified by tracking particles. In effect, having an ensemble of trajectories over large times and distances, the mean displacement gives the advective velocity $(\langle \bm{x} \rangle = \bm{v}_\text{p} t)$ while the growth of the centered mean square displacement $\sigma^2 = \langle \vert \bm{x} - \bm{v}_\text{p} t \vert^2 \rangle$ with time unveils the kind of transport. $\sigma^2 \propto t^\alpha$ with $\alpha = 1$ is normal dispersion while anomalous dispersion means that $\alpha > 1$ (superdispersion) or $\alpha < 1$ (subdispersion).

Anomalous transport occurs when the assumptions made to explain the Fickian law do not hold. This is the case of very short time or length scales, when the motion involves correlations (e.g., typically below the free path or the pore size lengths). On large time scales, memory effects are sometimes facilitated by the complexity of porous media. Although homogeneous media do not foster any hierarchy of pipes, the effect of trapping events on the motion deserves investigation.

To study the effects of non-Gaussian statistics, the simplest model of Brownian random walk (with always the same walk lengths and times) is not enough. Deviations from normal dispersion can be interpreted in terms of a stochastic continuous time random walks model CTRW [31,32]. In a CTRW the length $\Delta x$ and duration of walks $\Delta t$, as well as the times elapsed between consecutive walks, are drawn from a probability distribution function $f_{(\Delta \bm{x},\Delta t, \Delta t_\text{p})}$ [5]. Walk lengths and/or pauses drawn from PDFs with no finite second moment lead to anomalous transport processes. In these cases, the behavior of particles is dominated by the largest displacements and/or the longest pauses, because the system memory about such rare events is never erased. What is a large displacement or a long pause depends on the observation scale.

We carry out numerical experiments (fig. 1) to check the aspect of the PDFs of pause and walking times in cases in which blockages are possible. An originality of our approach is that we simulate on a microscopic level the interplay between fluid, transported particles and porous media. Our experiments capture events like drag forces, interparticle collisions, bouncing and rolling, formation and destruction of clogging arches and blockages.

We use a micro-hydromechanical numerical method, the DEM-PFV [33,34] implemented in YADE [35] to model particle transport. This method applies for Stokes flows of incompressible and Newtonian fluids. The DEM-PFV couples the discrete element method [36] with a pore-scale finite-volume formulation that solves the flow problem.

In the PFV framework, finite volumes are defined from a weighted Delaunay tessellation of the granular assembly. Therefore, at a given time, the geometry of the packing defines an open network of connected voids: finite volumes are tetrahedra whose apexes coincide with particle centers (fig. 2). Each finite volume represents a pore and deforms as the solid particles move. In our case, not only the coarse particles of the granular assembly but also the transported ones define the network of connected voids and hence this network changes dynamically with particle transport.

Zoom In Zoom Out Reset image size

The deformation of a tetrahedron leads, in general, to a change in its volume that must be balanced by the fluid fluxes through its four faces, so the continuity equation applied to the fluid phase in each pore reads [34]

$$\begin{equation} \dot{V}_{\text{f}i} = \sum_{j=1}^4 q_{ij} = \sum_{j=1}^4 g_{ij} \frac{\left( p_i - p_j \right)}{l_{ij}}, \end{equation} \tag{ 2 }$$

where $\dot{V}_{\text{f}i}$ is the volumetric strain rate of pore i deduced from the particle velocities at the four apexes of the tetrahedron i, and q_ij is the flux between the pore i and one of its four adjacent pores j. As explained in eq. (2), inter-pore fluxes q_ij are related to the hydraulic pressure gradients between the corresponding pores. p_i, p_j is the hydraulic pressure in pore i, j, respectively, and l_ij is the distance between i and j. A Stokes flow leads to linear relationships between fluid fluxes and pressure gradients that are characterized by the hydraulic conductivities g_ij. For a given fluid dynamic viscosity, g_ij only depends on the hydraulic radius and the cross-section areas of constrictions [33].

For a full coupling, the forces applied by the fluid on solid particles have to be taken into account in the DEM in addition to inter-particle contact forces. The fluid forces result from the integration of fluid pressure and viscous stresses over the surface of each solid particle. The development of these surface integrals shows that the pressure force $\bm{F}_{ij}^{\text{p},k}$ and the viscous force $\bm{F}_{ij}^{\text{v},k}$ acting on a particle k due to a fluid flow between pores i and j are proportional to the product between the pressure drop $(p_i - p_j)$ and the surface areas $A_{ij}^k$, $A_{ij}^f$ and $\gamma_{ij}^k$ (fig. 2) that describe the geometrical features of the interface plane and the constriction between two tetrahedra:

$$\begin{eqnarray} \bm{F}_{ij}^{\text{p},k} & = A_{ij}^k \left( p_i - p_j \right) \bm{n}_{ij} \text{,} \nonumber\\ \bm{F}_{ij}^{\text{p},k} & = \frac{\gamma_{ij}^k}{\sum_{k=1}^3 \gamma_{ij}^k} A_{ij}^f \left( p_j - p_i \right) \bm{n}_{ij} \text{,} \end{eqnarray} \tag{ 3 }$$

with $\bm{n}_{ij}$ the unit vector pointing from pores i to j. $A_{ij}^k$, $A_{ij}^f$ and $\gamma_{ij}^k$ are relatively easily computed for spheres.

We use granular assemblies of monodispersed spheres of size $d_\text{c}=5 \cdot 10^{-3}\ \text{m}$ as porous media. We report results from random, quasi-crystalline and perfect crystalline structures included in parallelepiped periodic simulation boxes. Quasi-crystalline structures (face-centered cubic (FCC*) and simple cubic (SC*)) are generated from the corresponding crystalline ones, by enlarging unit cells and slightly displacing particles from crystalline positions. $\phi_\text{FCC*}=0.30$ and $\phi_\text{SC*}=0.53$ are slightly higher than in perfect crystals but the tortuosities are almost the same. However, as the CSD is modified, blockages are only possible in some of the constrictions. On the contrary, in crystalline structures (face-centered cubic (FCC), $\phi_\text{FCC}=0.26$) the formation of blockages is not conditioned by the porous structure since identical constrictions are found all over the medium.

Fine particles of size $d_\text{f} = d_\text{c}/8$ are randomly and uniformly introduced in the voids of the assembly up to an average concentration C. Hookean contact is used to model interactions between solid particles and all the particles are frictionless. $d_\text{f}/d_\text{c}$ is below the maximum value that a frictionless particle can have in order to pass through the smallest possible constriction found in FCC media $(2 \sqrt{3}-1)$. Therefore, fine particles cannot get trapped alone. However, they can form a collective clogging arch when they locate in a constriction as shown in the inset of fig. 1(a). Obviously for a given CSD, blockages are more likely at high concentrations since clogging is basically a matter of the probability of the presence of particles [37,38]. We have checked this by changing the total amount of fine particles (i.e., the average C). We report exclusively the longitudinal velocity and dispersion (in the direction of $-\boldsymbol{\nabla} \bm{p}$). In normal dispersion, $\sigma^2 = 2 D t$ and the Péclet number $Pe=v_\text{p} d_\text{c} / D$ gives an idea of the importance of advective with respect to dispersive transport. We use $v_\text{0,w}$ as a reference for velocities and $t_\text{0} = d_\text{c} / v_\text{0}$ for time.

In SC* media, fine particles cannot form blockages and after some time the displacements of particles are uncorrelated because of collisions. Then we observe that the transport is normal, $\sigma^2 \propto t$ (fig. 4(a)). In fig. 3 we show how $v_\text{w}/v_\text{0,w}$ decreases with C, whereas $v_\text{p}/v_\text{w}$ grows because a higher number of transported particles in the voids increases $\mathcal{T}$ and reduces $\mathcal{T}_\text{p}$. These results are similar to those of particle percolation by gravity [27–29].

Zoom In Zoom Out Reset image size

In FCC* and RND media at low C, blockages are very unlikely but they may last a long time. Sometimes this time is much longer than that of observation, so that we find that some particles never move, while the remaining ones are normally transported. In effect, the data shown in figs. 3 and 4 correspond exclusively to transported particles. Regarding the growth of $\sigma^2$ with t (fig. 4(b), (c)), if the whole sample of fine particles is considered and a few of them are trapped, the center of mass moves more slowly than that of transported particles. Therefore, the contribution of the dispersion to $\sigma^2$ is hidden by a "ballistic" effect $(\alpha = 2)$. In effect, if a fraction r of particles is transported as the ADE states (with $v_\text{p}$ and D) while the others are stopped, after a time t, we find that $\langle x \rangle = r v t$ and $\sigma^2 = 2 r D t + r (1-r) v^2 t^2$. When $t \gg 2D/((1-r) v^2)$ then $\sigma^2 \propto t^2$. If the trapped particles are discarded, then only the advection and dispersion mechanisms are observed. In summary, in this case the transport is normal but we cannot forget that particles have a small chance of being trapped and that if they do, they will rest for a long time.

Zoom In Zoom Out Reset image size

In FCC media there are many traps so that the probability for blockage formation grows very quickly with C but, conversely, blockages are easily destroyed, too, since there are many particles in the voids and the fluid flow changes very quickly. As a result, on the time scale of observation we find that blockages appear and disappear and the growth of $\sigma^2$ with time is not that of normal dispersion, as shown in fig. 4(d). In the framework of CTRW this motion can be seen as a succession of walks (displacements between consecutive blockages) and pauses (times that particles remain in a blockage). During walks particles move on average at an effective velocity (fig. 3) $v_\text{p}^, > v_\text{p}$ (since $v_\text{p}$ implicitly includes the pauses) and disperse according to the normal transport paradigm.

Then the kind of transport on a longer time scale depends on whether or not the PDFs of pause and walking times have finite second moment. Some examples of power laws with non-finite second moment have been found in other situations [39], but in this case there is a reason to consider exponential decays, i.e., $f_{(t)} = (1/\tau_\text{p,w}) \exp{(-t/\tau_\text{p,w})}$ being $\tau_\text{p,w}$ the time scale factors of pauses and walks, respectively. In effect, if particles pass through all the constrictions in similar conditions (number of neighbors, velocities, constant probability of finding other particles blocking the constriction, etc.) and trapping events occur independently of each other, this corresponds to a homogeneous Poisson process, and the time between each pair of consecutive events follows an exponential distribution. Similar assumptions have been made in silo discharge problems, where a constant probability of clogging during the whole avalanche duration is assumed [40]. Regarding blockage destructions, an exponential decay also means a constant probability of blockage failure and non-dependence between consecutive events.

The PDFs of walking and pause times computed from our experiments fit quite well with exponential decays (except for the pause time of FCC $(C=0.009)$, probably as a consequence of the insufficient size of the statistical sample). In fig. 5 we plot a cumulative distribution function $F_{(t)} = \int_t^\infty f_{(t)}\,\mathrm{d} t$ since this helps to identify power law and exponential decays [39]. For the latter we expect $F_{(t)} = \exp{(-t/\tau)}$. The corresponding least squares approximation is also shown in fig. 5. Exponential PDFs lead to normal transport at time $t \gg \tau_\text{p}, \tau_\text{w}$. Therefore, $\tau_\text{p}$ and $\tau_\text{w}$ establish the time scales on which the memory effect of blockages formation is erased. These scales are, in the end, determined by the number of traps in the medium and the probability of finding several particles in the same pore (i.e., on C). We have found (fig. 6) that the higher C, the shorter $\tau_\text{w}$, with $\lim_{C \to 0} \tau_{\text{w}(C)} = \infty$ and $\tau_{\text{w}_{(C_\text{max})}} = \tau_{\text{w,min}}$. $\tau_\text{p}$ decreases slowly with C but this relationship deserves further investigation. At high C, the local changes of fluid flows are presumably more frequent and intense, but with more particles in the blockages, those that are very embedded may need more time to be released. For $C=0.012$, and especially for $C=0.009$, the time scale is much longer than our observation time, and this anticipates that media having wider CSDs blockages (i.e., fewer traps) lead to higher $\tau_\text{w}$, indicating how to approach FCC* cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With these results in mind, some already known transport models can be used to study particle transport when mechanical clogging of constrictions is expected.

In particular, a continuum transport model with 1st-order kinetic mass transfer laws can be used. This needs the inclusion on the right-hand side of eq. (1) of laws of type $f_\text{f} = - C/\tau_\text{w}$ and $f_\text{e} = + S/\tau_\text{p}$, being $S_{(x,t)}$ the concentration of trapped particles, and $\partial S/\partial t = - f_\text{f}-f_\text{e}$. We have compared the predictions of this model to the results of DEM-PFV simulations (fig. 7(a)). A steady-state solution of the corresponding 1D problem with periodic boundary conditions is $C_{(x,t)} = \tau_\text{w}/(\tau_\text{w} + \tau_\text{p})$ which is in the end the expected value of the average fraction of traveling particles $\langle r \rangle$. $\langle r \rangle$, directly measured in the simulations, is compared to $\tau_\text{w}/(\tau_\text{w} + \tau_\text{p})$ in fig. 8 (where $\tau_\text{w,p}$ come from the statistical fitting). Furthermore, $v_p / v_\text{p}^{\prime}$, also shown in fig. 8 should be equal to $\langle r \rangle$, since, on average, particles are expected to move at $v_\text{p}^{\prime}$ for a fraction of time $\tau_\text{w} /(\tau_\text{w} + \tau_\text{p})$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, CTRW models can also be used. We have also compared the DEM-PFV simulations to a CTRW model (fig. 7(b)) in which we consider that at a given time, each particle is either traveling or trapped. In the first case, as the motion of particles between consecutive blockages is supposed to be normal, its path can be described (until the end of the corresponding traveling period) as a succession of advective $+v_\text{p}^{\prime} \Delta t$ and random steps $\pm \Delta x$ that occur at fixed $\Delta t = \Delta x^2 / (2D) \ll \tau_\text{w,p}$. During a trapping period, the particle does not move. The duration of traveling and waiting periods (i.e., the number of steps that particles move or do not move) is drawn from exponential laws of the parameters $\tau_\text{w}$ and $\tau_\text{p}$, respectively. The inclusion of these pauses is what makes the difference with respect to an usual random walk model, and this explains the behavior of the ensemble of particles on the short time scale. Moreover, the time scale for normal transport (i.e., $\sigma^2 \propto t$) and the effective transport parameters can be obtained in turn from this model.

Thanks to the DEM-PFV numerical method, "ab initio" simulations of particle transport in dense granular assemblies have been performed. Fluid-solid interactions were considered on a particle scale (including transported particles). We find reasonably good fits with exponential decays for the duration of trapping events and of displacements between consecutive trapping events. This establishes the scale for normal transport. On shorter scales we show that an ADE with 1st-order kinetic mass transfer laws (a model originally introduced for other systems) is valid to account for the intermittent formation of blockages during the transport of moderately large, passive, solid particles that interact only mechanically with both the fluid and the porous medium. We also show, how, for some porous structures, exponential factors and transport parameters depend on the concentration.

This work has been partially supported by the LabEx Tec 21 (Investissements d'Avenir - grant agreement nANR-11-LABX-0030).