# Particulate transport in porous media at pore-scale.
## Part 1: Unresolved-resolved four-way coupling
CFD-DEM

Laurez Maya Fogouang, Laurent André, Cyprien Soulaine

► To cite this version:

Laurez Maya Fogouang, Laurent André, Cyprien Soulaine. Particulate transport in porous media at pore-scale. Part 1: Unresolved-resolved four-way coupling CFD-DEM. Journal of Computational Physics, 2025, 521, 10.1016/j.jcp.2024.113540. insu-04796295

HAL Id: insu-04796295

https://insu.hal.science/insu-04796295v1

Submitted on 21 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Particulate transport in porous media at pore-scale. Part 1:
Unresolved-resolved four-way coupling CFD-DEM

Laurez Maya Fogouang<sup>a,b,*</sup>, Laurent André<sup>a,b</sup>, Cyprien Soulaine<sup>a,**</sup>

## ARTICLE INFO

Keywords:

CFD-DEM model

Particulate flows

Unresolved coupling

Resolved coupling

Pore-scale modeling

Pore-clogging

## ABSTRACT

Computational Fluid Dynamics - Discrete Element Method (CFD-DEM) is a powerful approach to simulate particulate flow in porous media at the pore-scale, and hence decipher the complex interplay between particle transport and retention. Two separate CFD-DEM approaches are commonly used in the literature: the unresolved (particle smaller than the grid cell size) and the resolved (particle bigger than the grid cell size) approach. In this paper, we propose a novel CFD-DEM coupling approach that combines both unresolved and resolved coupling. Our new modeling technique allows for the simulation of particulate flows in complex pore morphology characteristic of porous materials. It relies on an efficient searching strategy to find grid cells covered by the particles and on an appropriate calculation of the fluid-solid momentum exchange term. The robustness and efficiency of the computational model are demonstrated using cases for which reference solutions – analytical or experimental – exist. The new unresolved-resolved four-way coupling CFD-DEM is used to investigate pore-clogging and permeability reduction due to the sieving and bridging of particles.

## 1. Introduction

There is a growing interest in simulating particulate flow in porous media at pore-scale [1–3]. Applications include diesel particulate filters [4,5], filtration of solid aerosols [6], groundwater treatment [7,8], production and injection of fluids in the subsurface in geothermal operations [9,10], carbon dioxide storage in underground reservoirs [11], and rock-based hydrogen production and storage [12,13]. In porous media, particles are suspended in a carrier fluid or attached at the fluid-solid interface clogging the pores [14–16]. Pore-clogging can be a desired phenomenon in filtration processes as it increases filter efficiency [17], or on the contrary, an unwanted effect during subsurface exploitation because it induces permeability damage near wellbore that reduces drastically the operating times [18].

The heterogeneity of porous media structure made of a network of pores of various sizes and shapes leads to complex feedback between particle transport, deposition, and pore-clogging. The description of the clog formation at the pore-scale – where the pore space is fully resolved – enables to decipher and characterize these processes [19]. Three mechanisms are responsible for pore-clogging: (i) sieving, (ii) bridging, and (iii) aggregation of particles [20]. Sieving refers to the blockage of particles based on size exclusion. Bridging consists in the formation of bridge-like structures composed of a few particles arriving at the same time at the

<footnote><sup>a</sup> Institut des Sciences de la Terre d'Orléans (ISTO), UMR7327, Université d'Orléans, CNRS, BRGM, F-45071, Orléans, France</footnote>

<footnote><sup>b</sup> BRGM, F-45060, Orléans, France</footnote>

<footnote>* Corresponding author at: Institut des Sciences de la Terre d'Orléans (ISTO), UMR7327, Université d'Orléans, CNRS, BRGM, F-45071, Orléans, France.</footnote>

<footnote>** Principal corresponding author.</footnote>

<footnote>E-mail addresses: laurez.maya@cnrs-orleans.fr (L. Maya Fogouang), cyprien.soulaine@cnrs-orleans.fr (C. Soulaine).</footnote>
pore entrance in which the flow is converging. Aggregation of particles results from successive depositions of colloidal particles driven by long-range electrochemical forces between the fluid, particles, and the solid surface [21]. Thus, a good comprehension of those clogging mechanisms passes through an accurate method for modeling particle transport at the pore-scale. In this paper, we focus on the first two mechanisms (sieving and bridging) that depend only on hydromechanical forces.

Several computational approaches exist in the literature for solving the carrier fluid flow combined with particulate transport in porous media. They are divided into two main families: Euler-Lagrange and Euler-Euler. In the first one, the fluid flow is solved on an Eulerian grid, and the particles are tracked individually using a Lagrangian frame [22]. In the second one, the motion of the particles is described by continuum equations solved on the Eulerian grid. Accurate Euler-Euler models for simulating particle retention in porous media up to pore-clogging are still scarce because the description of deposition and plugging relies on constitutive relationships that depend on the physicochemical properties of the suspension [23]. Euler-Lagrange approaches based on the Discrete Element Method (DEM) resolve the hydrodynamic interactions between the fluid and the particles and the mechanical contacts between the particles. In this kind of model, the fluid flow can be solved by the Lattice-Boltzmann (LBM) method or by Navier-Stokes equations using Computational Fluid Dynamics (CFD) including the Finite Volume Method. For example, Li and Prigiobbe [24] and Zhou et al. [25] used an LBM-DEM approach coupled with an Immersed Boundary Method (IBM) for the migration of fine particles. CFD-DEM is commonly used to simulate particle motion and its feedback on the flow [26].

State-of-the-art Euler-Lagrange methods including CFD-DEM rely on the so-called four-way coupling which states that the fluid transports the particles, the presence of particles impacts the fluid flow, and particle-particle and particle-wall interactions are considered [27]. Two different fluid-particle coupling approaches exist whether particle diameter, $D_p$, is bigger or not than the grid size, $\Delta x$. On the one hand, the unresolved coupling approach is used if the particle size is smaller than the grid resolution. The drag force between fluid and particle is calculated using a sub-grid model. On the other hand, in resolved coupling, the particle size is bigger than the grid resolution and the fluid-particle shear stress forces acting on the particle are an output of the simulation. LBM-DEM is not commonly used for unresolved particles because the locality of LBM calculation is lost and numerical issues may occur when particles are close to each other [28]. Su et al. [29] and Elrahmani et al. [30] use resolved CFD-DEM to simulate sieving and bridging mechanisms. Nan et al. [31] use a coupled Volume Of Fluid method with DEM to investigate clogging mechanisms during concrete seepage under multiphase flow (water+air) with comparison to experiments. For good accuracy of their model, they need cells 3 to 6 times larger than particle sizes to satisfy the relative size requirement.

Limitations exist in both resolved and unresolved approaches (see Table 1). Pirker et al. [32], Marshall and Sala [33] reported that unresolved coupling leads to significant errors and instabilities if $D_p/\Delta x > 1/3$. Despite its faster calculation time, this approach presents discontinuities in the calculation of the particle velocity because the particle oscillates when it crosses two adjacent cells that have different fluid velocities [34]. Importantly, sieving cannot be simulated using unresolved coupling because the grid size is necessarily smaller than the pore-throat and particle size. In resolved coupling, the size ratio has to be at least $D_p/\Delta x > 10$ to ensure an accurate resolution of the particle surface on the Eulerian grid [35,3]. This results in large computational times that increase with the particle number which limits this method for simulating large quantities of particles. Moreover, the complexity of the porous geometry with pores of different sizes often requires computational grids made of cells with different sizes, and the same particle might be resolved in some parts of the domain and unresolved in others.

Some approaches intend to overcome the aforementioned limitations by developing hybrid unresolved-resolved approaches. For example, the semi-resolved CFD-DEM described in Wang et al. [36] corrects the fluid velocity around the particle and the volume fraction in the drag force model using kernel-based approximations. Their model half-theoretical half-empirical leads to inaccurate calculations of the fluid-particle forces for dense suspensions of large particles. Kuruneru et al. [37] proposed a mixed resolved-unresolved CFD-DEM approach that uses an immersed boundary method (IBM), a specific contact handling algorithm to compute particle contact forces, and a Brinkman penalization technique for the momentum sink term of the fluid phase. However, their model works only for regular grids which is limiting for heterogeneous porous media with confined pores.

In this work, we introduce a novel hybrid unresolved-resolved CFD-DEM four-way coupling method to simulate particulate flows in porous media at the pore-scale with a cloud of spherical, rigid particles as the dispersed phase. The proposed model includes a drag force to encounter both particles resolved and unresolved, the associated hydrodynamic torque, an algorithm for identifying particle-fluid cells, and a diffusive layer to smooth the particle-fluid interface. The paper is organized as follows. In section 2, we describe the governing equations of the fluid and particles and the coupling algorithms. In section 3, we show the accuracy, consistency, and efficiency of our approach using cases for which reference solutions exist. Then, in section 4, we use our new CFD-DEM approach to investigate clogging mechanisms due to the injection of particles in porous media. Finally (section 5), we close with a summary and concluding remarks.

## 2. Computational model

In this section, we present our new resolved-unresolved four-way coupling CFD-DEM. First, we introduce the CFD approach to model fluid flow in an Eulerian grid (Section 2.1). Then, we introduce the DEM approach to simulate the particle displacement in a Lagrangian frame (Section 2.2). Finally, we show our strategy to couple CFD and DEM together (Section 2.3).

### 2.1. Fluid motion in the CFD Eulerian grid

The Eulerian computational grid is used to solve fluid flow accounting for the presence of the particles. The latter is described by a local porosity field $\epsilon$, defined as,
<caption>**Table 1**</caption>

Features of the resolved, semi-resolved [36], unresolved and unresolved-resolved CFD-DEM. DNS = Direct Numerical Simulation, FVM = Finite Volume Method.

<table><thead><tr><th></th><th>Resolved<br>CFD-DEM</th><th>Semi-resolved<br>CFD-DEM</th><th>Unresolved<br>CFD-DEM</th><th>Unresolved-resolved<br>CFD-DEM</th></tr></thead><tbody><tr><td>Particle surface</td><td>Resolved</td><td>Not resolved</td><td>Not resolved</td><td>Resolved for D<sub>p</sub>/Δx ≥ 1 and unresolved for D<sub>p</sub>/Δx < 1</td></tr><tr><td>Particle-fluid force</td><td>Particle-resolved DNS</td><td>Drag force model</td><td>Drag force model</td><td>Combined particle-resolved DNS and drag force model</td></tr><tr><td>Background velocity</td><td>Fluid velocities<br>resolved in FVM cells</td><td>Fluid velocities<br>resolved in<br>neighboring FVM<br>cells</td><td>Fluid velocities in the<br>local FVM cell only</td><td>Fluid velocities calculated in FVM cells<br>around resolved particles and in local<br>FVM cells for unresolved particles</td></tr><tr><td>Particle to cell<br>ratio D<sub>p</sub>/Δx</td><td>> 10</td><td>≈ 1</td><td>&lt; 1/3</td><td>No restriction</td></tr></tbody></table>

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>a) Resolved approach</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>b) Unresolved approach</figcaption></figure>

<caption>**Fig. 1.** Mapping of the porosity field, $\epsilon$, on the Eulerian grid for resolved particles (a) and for unresolved particles (b).</caption>

$$ \epsilon = \begin{cases} 1, & \text{if the cell is occupied by fluids only,} \\ ]0, 1[, & \text{if the cell contains a fluid-solid aggregate or a fluid-solid interface,} \\ 0, & \text{if the cell is occupied by solids only.} \end{cases} \qquad (1) $$

The two situations illustrated in Fig. 1 can exist concomitantly whether the particle diameter is larger or smaller than the cell size. On the one hand, in unresolved cases (see Fig. 1.b), particles are smaller than the cell size. A cluster of particles can occupy a computational cell and in this case, $\epsilon \in ]0, 1[$. On the other hand, in resolved cases (see Fig. 1.a), the particle is larger than the cell size, and its shadow covers an ensemble of cells. The particle shadow is obtained by a projection on the Eulerian-covered cells identified through a searching algorithm described in Section 2.3.1. In this situation, intermediate values of the phase indicator $\epsilon \in ]0, 1[$, correspond to cells occupied by the fluid-particle interfaces. If a cell is fully covered by a particle, we impose a minimum porosity value, $\epsilon_{min} = 0.001$, instead of $\epsilon = 0$ to have flow equations defined everywhere in the computational domain regardless of the cell content [38].

The fluid motion is modeled by solving incompressible Volume-Averaged Navier-Stokes (VANS) equations [39,40]. The mass balance equation for the fluid phase reads

$$ \frac{\partial(\epsilon \rho_f)}{\partial t} + \nabla \cdot (\epsilon \rho_f \mathbf{v}^f) = 0, \qquad (2) $$

where $\rho_f$ is the fluid density, and $\mathbf{v}^f$ is the cell-averaged fluid velocity and $t$ is the time.

The fluid momentum balance equation is:

$$ \frac{\partial(\epsilon \rho_f \mathbf{v}^f)}{\partial t} + \nabla \cdot (\epsilon \rho_f \mathbf{v}^f \mathbf{v}^f) = -\epsilon \nabla p + \epsilon \rho_f \mathbf{g} + \epsilon \nabla \cdot (\mu_f (\nabla \mathbf{v}^f + (\nabla \mathbf{v}^f)^T) - \epsilon^2 \mu_f (\mathbf{v}^f - \bar{\mathbf{v}}^p)/K, \qquad (3) $$

where $p$ is the fluid pressure, $\mu_f$ is the fluid dynamic viscosity, $\mathbf{g}$ is the gravitational acceleration, $\bar{\mathbf{v}}^p$ is the averaged particle velocity on the Eulerian grid, and $K$ is the local cell-permeability. The operator $(\cdot)^T$ refers to the transpose of the matrix. The last term of the right-hand side is a drag force corresponding to the flow resistance due to the presence of particles. VANS momentum tends
<figure><img src="image_5.png" /><figcaption>Fig. 2. Illustration of the forces acting on particle i from contacting particle j and the wall w.</figcaption></figure>

asymptotically towards the Navier-Stokes equation in regions that contain fluid only ($\epsilon = 1$ and the drag force vanishes) and toward Darcy's law in cells containing fluid-solid aggregates ($0 < \epsilon < 1$) because inertia and viscous dissipation are negligible in front of the drag force [41].

The transition between these two asymptotic behaviors is obtained using a cell-permeability that varies with the cell porosity in a way that $K^{-1} \to 0$ if $\epsilon = 1$ and to a finite value if $\epsilon < 1$. In the latter case, K is determined differently whether the cell is occupied by resolved or unresolved particles. In resolved cases, particles are seen as low-porosity ($\epsilon = \epsilon_{min}$) low-permeability media and the drag force acts as a penalization term to drop the velocity within the occupied cells near zero value and to approach no-slip condition on the particle-fluid interfaces [42,43]. In this case, K is calculated using a modified Kozeny-Carman formula [44],

$$K = K_0 \frac{e^3}{(1 - \epsilon)^2}, \qquad (4)$$

where $K_0$ is a sufficiently low given permeability. In unresolved cases, K is the permeability field of the ensemble of particles seen as a porous medium and obtained using the drag forces applied on particles (see section 2.2.1).

## 2.2. Particle motion in the DEM Lagrangian frame

The motion of a cloud of particles is solved in a Lagrangian frame using a four-way coupling Discrete Element Method (DEM). In this approach, the particle flow is driven by interactions with the carrier fluid, and particle-particle and particle-wall interactions. The total velocity, $v_i^p$, of particle i reads,

$$v_i^p = U_i^p + \omega_i \times r_i, \qquad (5)$$

where $U_i^p$ and $\omega_i$ are the translational and the angular velocity of particle i, respectively, and $r_i$ is the position vector. The description of the translational and rotational motions in DEM is based on Newton's second law applied to a spherical particle of mass $m_i$ and moment of inertia $I_i$ in contact with j objects ($j = 1, 2, ..., n_i^c$ particles and walls). The balance of forces is illustrated in Fig. 2. We have [45]:

$$m_i \frac{dU_i^p}{dt} = \sum_j^{n_i^c} F_{ij}^c + F_i^f + F_i^g, \qquad (6)$$

and,

$$I_i \frac{d\omega_i}{dt} = \sum_j^{n_i^c} M_{ij}^c + M_i^{hyd}, \qquad (7)$$

where $F_i^f$ and $F_i^g$ are the particle-fluid interactions and gravitational forces acting on particle i at time t, respectively. $F_{ij}^c$, $M_{ij}^c$ and $M_i^{hyd}$ are the contact forces, the contact torques, and the hydrodynamic torque acting on particle i, respectively. $n_i^c$ is the number of objects (particles and walls) in contact with particle i at time t. All these forces and torques are described in the following.

### 2.2.1. Hydrodynamic forces

In the four-way approach, the surrounding fluid in contact with a particle creates particle-fluid interactions on the particle i among which the drag, the pressure gradient, the virtual mass, and the lift forces. We only considered the drag force. The drag force
<caption>**Table 2**</caption>

Different formulations of the fluid-solid momentum exchange coefficient $\beta$ are implemented, where $C_d$ is the drag coefficient, $Re_p$ is the particle Reynolds number, $D_{p,i}$ is the diameter of particle $i$ and $V_{p,i}$ is the volume of particle $i$.

<table><thead><tr><th></th><th>β formulation</th><th>Description</th></tr></thead><tbody><tr><td>Stokes law</td><td>β = 3πμ<sub>f</sub> D<sub>p,i</sub> (1-ε) / (V<sub>p,i</sub><sup>3</sup>)</td><td>For a single particle with Re<sub>p</sub> &lt; 1</td></tr><tr><td>Kozeny-Carman law</td><td>β = 180 (1-ε)<sup>3</sup> / (ε<sup>2</sup>) (μ<sub>f</sub> / D<sub>p,i</sub><sup>2</sup>)</td><td>For dense suspensions (ε &lt; 0.8) and Re<sub>p</sub> &lt; 1</td></tr><tr><td>Ergun law [46]</td><td>β = 150 (1-ε)<sup>2</sup> / (ε) (μ<sub>f</sub> / D<sub>p,i</sub><sup>2</sup>) + 1.75(1-ε) (μ<sub>f</sub> / D<sub>p,i</sub><sup>2</sup>) |v<sup>f</sup> - v<sub>i</sub><sup>p</sup>|</td><td>For dense suspensions (ε &lt; 0.8) and Re<sub>p</sub> &gt; 1</td></tr><tr><td>Wen and Yu law [47]</td><td>β = 3/4 C<sub>d</sub> (ε(1-ε) / D<sub>p,i</sub><sup>2</sup>) ρ<sub>f</sub> |v<sup>f</sup> - v<sub>i</sub><sup>p</sup>| e<sup>-2.65</sup>, with<br>C<sub>d</sub> = <span style="font-size: 2em; vertical-align: middle;">{</span> <table style="display: inline-table; vertical-align: middle; border: none; text-align: center; line-height: 1.2em;"><tbody><tr><td style="border: none; padding: 0 0.2em; border-bottom: 1px solid black;">24 / (Re<sub>p</sub>) (1 + 0.15Re<sub>p</sub><sup>0.687</sup>)</td><td style="border: none; padding: 0 0.2em;">if Re<sub>p</sub> ≤ 1000,</td></tr><tr><td style="border: none; padding: 0 0.2em;">0.44</td><td style="border: none; padding: 0 0.2em;">if Re<sub>p</sub> &gt; 1000.</td></tr></tbody></table></td><td>For dilute suspensions (ε ≥ 0.8)</td></tr></tbody></table>

applied on particles is calculated through the appropriate porosity-permeability relationship if particles are unresolved and with the summation of the stress divergence terms over the cells containing the particle for resolved coupling. We have,

$$F_i^f = \begin{cases} V_{p,i}\beta (v^f - v_i^p) \frac{1}{(1-\epsilon)}, & \text{if unresolved coupling.} \\ \sum_{k}^{n^k} V_{c,k} (1-\epsilon_k) (-\rho_f \nabla p_k + \nabla \cdot \boldsymbol{\tau}_k), & \text{if resolved coupling.} \end{cases} \quad (8)$$

where $n^k$ is the total number of cells covered by the particle $i$, $V_{c,k}$ the volume of cell $k$, $\boldsymbol{\tau}_k = \mu_f (\nabla v^f + (\nabla v^f)^T)$ the fluid shear-rate tensor of cell $k$, $p_k$ and $\epsilon_k$ the fluid pressure and the porosity in cell $k$, $V_{p,i}$ the volume of the particle $i$, and $\beta$ is the fluid-particle momentum exchange coefficient. In the code, we implemented different formulations to calculate the fluid-solid momentum exchange coefficient (see Table 2).

The gravitational force for a spherical particle considering buoyancy effects reads,

$$F_i^g = \frac{4}{3}\pi\rho_i R_i^3 \mathbf{g} - \frac{4}{3}\pi\rho_f R_i^3 \mathbf{g} = m_i \mathbf{g} \left(1 - \frac{\rho_f}{\rho_i}\right), \quad (9)$$

where $\rho_i$, is the particle's density and $R_i$ is the particle's radius.

### 2.2.2. Contact forces

A particle can be in contact with other particles (particle-particle interactions) or with walls (particle-wall interactions). We consider rigid spherical particles, and the contact between two elements is not at a single point but a finite area corresponding to the overlapping of the two objects [27]. The overlapping distance obeys the Hertzian spring-sliding-dashpot model [48,26] in which the spring realizes the elastic deformation, the dashpot realizes the viscous dissipation, and the slider realizes the frictional dissipation. These effects act on particles through the stiffness $k$, the damping coefficient $\eta$, and the friction coefficient $\mu$. These parameters are based on particle properties including radius, mass, Poisson coefficient, and Young modulus (see Table 3). The contact force $F_i^c$ has two components – a normal and a tangential – described as follows:

#### (a) Normal component

The normal component of the inter-particle and particle-wall contact ($F_{nij}^c$) acting on particle $i$ in contact with object $j$ (particle or wall) is given by the sum of the forces related to the spring and dashpot

$$F_{nij}^c = (-k_{nij} |\delta_{nij}|^{3/2} - \eta_{nij} v_{ij}^p \cdot \mathbf{n}_{ij}) \mathbf{n}_{ij}, \quad (10)$$

where $k_{nij}$ and $\eta_{nij}$ are respectively the equivalent normal stiffness and damping coefficients of particles $i$ with object $j$. $|\delta_{nij}|$ is the normal overlapping distance given by,

$$|\delta_{nij}| = \begin{cases} R_i + R_j - |\mathbf{p}_j - \mathbf{p}_i|, & \text{for particle-particle contact,} \\ R_i - |\mathbf{p}_i - \mathbf{p}_w|, & \text{for particle-wall contact,} \end{cases} \quad (11)$$

where $\mathbf{p}_i$ and $\mathbf{p}_j$ are the position vector of particle $i$ and $j$, respectively. The vector $\mathbf{p}_w$ corresponds to the nearest point to $\mathbf{p}_i$ located on the wall. The relative velocity, $v_{ij}^p$, is given by $v_{ij}^p = \mathbf{v}_i^p - \mathbf{v}_j^p$, where for a wall, $v_j^p$ is the slip velocity of any sphere-wall contact point. The unit vector $\mathbf{n}_{ij}$ points either from the center of particle $i$ to that of particle $j$ in particle-particle contact, or points from the wall to the computational domain and is normal to the wall for particle-wall contact.

#### (b) Tangential component

The tangential component of the particle-particle and particle-wall forces ($F_{tij}^c$) acting on particle $i$ depends on the tangential overlap, $\delta_{tij}$, and on the tangential slip velocities, $v_{tij}^p$, according to
<caption>Table 3</caption>

Parameters for the DEM model. For the equivalent quantities, the case of particle-wall interaction is obtained assuming that the wall has infinite radius and mass. v is the Poisson's ratio, E is the Young modulus and α is given by the coefficient of elasticity.

<table><thead><tr><th>Parameter</th><th>Equation</th></tr></thead><tbody><tr><td>Equivalent normal stiffness</td><td>k<sub>ni,j</sub> = <sup>4</sup>&frasl;<sub>3</sub>√<span style="text-decoration: overline;">r<sub>ij</sub></span>E<sub>ij</sub></td></tr><tr><td>Equivalent normal damping coefficient</td><td>n<sub>ni,j</sub> = αδ<sup>0.25</sup><sub>ni,j</sub>√<span style="text-decoration: overline;">m<sub>ij</sub>k<sub>ni,j</sub></span></td></tr><tr><td>Equivalent tangential stiffness</td><td>k<sub>ti,j</sub> = 8G<sub>ij</sub>√<span style="text-decoration: overline;">r<sub>ij</sub></span>δ<sub>ni,j</sub></td></tr><tr><td>Equivalent tangential damping coefficient</td><td>η<sub>ti,j</sub> = η<sub>ni,j</sub></td></tr><tr><td>Equivalent Young modulus</td><td>E<sub>ij</sub> = (<sup>1-v<sub>i</sub><sup>2</sup></sup>&frasl;<sub>E<sub>i</sub></sub> + <sup>1-v<sub>j</sub><sup>2</sup></sup>&frasl;<sub>E<sub>j</sub></sub>)<sup>-1</sup></td></tr><tr><td>Equivalent shear modulus</td><td>G<sub>ij</sub> = (<sup>2(1+v<sub>i</sub>)(2-v<sub>i</sub>)</sup>&frasl;<sub>E<sub>i</sub></sub> + <sup>2(1+v<sub>j</sub>)(2-v<sub>j</sub>)</sup>&frasl;<sub>E<sub>j</sub></sub>)<sup>-1</sup></td></tr><tr><td>Equivalent mass</td><td>m<sub>ij</sub> = (<sup>1</sup>&frasl;<sub>m<sub>i</sub></sub> + <sup>1</sup>&frasl;<sub>m<sub>j</sub></sub>)<sup>-1</sup></td></tr><tr><td>Equivalent radius</td><td>r<sub>ij</sub> = (<sup>1</sup>&frasl;<sub>r<sub>i</sub></sub> + <sup>1</sup>&frasl;<sub>r<sub>j</sub></sub>)<sup>-1</sup></td></tr></tbody></table>

$$F_{tij}^c = -k_{tij} \delta_{tij} - \eta_{tij} \mathbf{v}_{tij}^p, \qquad (12)$$

where $k_{tij}$ and $\eta_{tij}$ are the tangential equivalent stiffness and damping coefficients of particle i with object j. The formula for the tangential overlap $\delta_{tij}$ for particle-particle and particle-wall contacts is found on [49]. The tangential slip velocities are given by $\mathbf{v}_{tij}^p = \mathbf{v}_{tij}^p - (\mathbf{v}_{tij}^p \cdot \mathbf{n}_{tij})\mathbf{n}_{tij} + (R_i\omega_i + R_j\omega_j) \times \mathbf{n}_{tij}$. Notes that for particle-wall contact with immobile walls, $\omega_j = 0$.

If the relation, $|F_{tij}^c| \ge \mu|F_{ni,j}^c|$, is satisfied with $\mu$ the friction coefficient, then particle i is sliding over object j and the tangential force is modeled by Coulomb-type sliding friction,

$$F_{tij}^c = -\mu \frac{\|F_{ni,j}^c\|}{\|\delta_{tij}\|} \qquad (13)$$

To reduce the computational costs, we consider that only the objects j (particles and walls) within a kernel centered on the centroid of particle i and with a radius equal to the particle diameter are candidates for particle-particle and particle-wall interactions. If an object j covers a cell that is overlapping with that radius, the object j is added to the interactions list of particle i.

### 2.2.3. Torques

The torque due to inter-particle and particle-wall contacts acting on a particle i is,

$$M_{ij}^c = R_i \mathbf{n}_{ij} \times \mathbf{F}_{tij}^c. \qquad (14)$$

Rolling is the dominant hydrodynamic mechanism that can cause particle removal from a wall under laminar flow. For an unresolved particle, the hydrodynamic torque is calculated when approaching an object j. For a resolved particle, the hydrodynamic torque is calculated over the cells covered by its shadow. The influence of the hydrodynamic shear on particle i is given by

$$M_i^{hyd} = \begin{cases} 1.4 R_i \mathbf{n}_{ij} \times \mathbf{F}_i^f, & \text{if unresolved coupling,} \\ \sum_{k}^{n_k} \mathbf{s}_k \times V_{c,k} (1 - \epsilon_k) (-\rho_f \nabla p_k + \nabla \cdot \boldsymbol{\tau}_k), & \text{if resolved coupling,} \end{cases} \qquad (15)$$

where $\mathbf{s}_k$ is the position vector relative to the particle i center and pointing to the center of the covered cell k. As the velocity increases with the distance from the obstacle, Torkzaban et al. [50] pointed out that the drag force effectively acts on the particle at a distance equal to $1.4 R_i$.

## 2.3. CFD-DEM coupling strategy and numerical implementation

The unresolved-resolved four-way coupling CFD-DEM is implemented within the open-source finite-volume toolbox OpenFOAM version 9 (https://www.openfoam.org). Our implementation is built on top of the existing denseParticleFoam solver coupled with the OpenFOAM internal DEM package to simulate dense particle packing. This solver has strong limitations for simulating particle transport in porous media at the pore-scale. First, it uses an unresolved approach. Second, numerical instabilities occur when particles cross the Eulerian cell faces. In addition to the DEM model in which we implemented the unresolved-resolved hydrodynamic drag and torque described in Section 2.2, we developed an efficient searching algorithm to identify the Eulerian cells covered by the shadow of the resolved particles (Sec. 2.3.1), we implemented a velocity-pressure solution algorithm free of numerical errors when a particle centroid crosses a cell face (Sec. 2.3.2), and we constrained the time-stepping with appropriate stability criteria (Sec. 2.3.3). The detailed algorithmic procedure of the unresolved-resolved four-way coupling CFD-DEM is illustrated in Fig. 3. Below is a step-by-step guide on how the coupling is made. (i) The coupling strategy is based on the four-way coupling: this includes both particle-fluid
<figure><img src="image_6.png" /><figcaption>Fig. 3. Detailed flowchart of the numerical procedure for solving the unresolved-resolved four-way coupling CFD-DEM. At t = 0, we set appropriate boundary conditions for the fluid and we define the particle-boundary interactions. Within a time-step, the CFD solver: i) updates the particle position and the associated volume fraction and momentum transfer term, ii) solve the pressure-velocity coupling using a PIMPLE algorithm. The particle position is obtained by the DEM package that computes collisions and force balance for every particle. The coupling between CFD and DEM is done by projecting the particle presence (i.e. volume fraction) onto the fixed grid using the search algorithm. The drag force applied on particles is computed differently either the particle is resolved or unresolved.</figcaption></figure>

interactions and particle-particle collisions and is often required in systems with dense particle suspensions or granular flows. (ii) The model uses particle-to-fluid coupling: the DEM solver computes particle positions, velocities, and forces (e.g., from collisions), and updates the fluid with the momentum transfer or volume displacement caused by the particles. This affects the fluid flow fields in the CFD solver. (iii) Coupling time step: typically, the CFD solver takes larger time steps than the DEM solver (since DEM involves particle collisions, the time step for the DEM solver is often much smaller than for the CFD solver to resolve collisions accurately). After every CFD time step, the DEM solver may perform multiple sub-steps. After each complete DEM step, the results are fed back into the CFD solver.

### 2.3.1. Mapping the resolved/unresolved local porosity: searching strategy of covered cells

A key aspect of the CFD-DEM coupling is the projection of the particle presence onto the Eulerian grid. In unresolved cases, the particle is smaller than the grid size and the local porosity is lower than 1. In a resolved case, the particle shadow covers an ensemble of cells. In both cases, the local porosity within a cell is obtained by calculating the sum of each volume of particle contained in the cell using

$$ \epsilon = 1 - \frac{\sum_j \Delta V_{p,j}}{V_c}, \qquad (16) $$

where $\Delta V_{p,j}$ is the volume of each particle j contained within the cell, and $V_c$ is the cell volume. In CFD-DEM, the cell label in which a particle centroid is located is an attribute of the particle. This attribute changes if the particle moves to a neighboring cell. Other attributes include particle mass and diameter. Therefore, in unresolved cases, if the particle is not overlapping two or more cells, then the volume of particles within computational cells is known and the mapping operation is straightforward.

In resolved coupling ($\Delta x < D_p$) and in unresolved coupling with particle overlapping 2 or more cells, however, the cells covered by the particle shadow have to be identified by a searching algorithm knowing the particle position and diameter. The process of identifying all the cells covered by the particles can highly impact the computational cost if a search method such as traversal search (going individually through all the cells of the Eulerian grid) is used, in particular, if the numbers of particles and grid cells are extremely large. Here, we present a peer-to-peer search algorithm whose objective is to identify efficiently the cells neighboring any given cell up to a certain distance.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Principle of the peer-to-peer search algorithm. (a)–(d) Illustrations of the steps for searching nearby grid points in a two-dimensional 4×4 grid. Points correspond to cell centers. Lines correspond to the path to a newly identified cell. Their color changes gradually at each iteration. From an initially identified cell (black empty circle), the algorithm searches for cells adjacent to the newly identified cells, and iterates. (e) The peer-to-peer search is an efficient way to identify cells covered by particles. The red lines describe the path the algorithm took to research covered cells. (f) The peer-to-peer search also works in three-dimensional structured and unstructured grids – here, a 3 × 3 × 3 regular domain. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

An example of the peer-to-peer algorithm applied to a two-dimensional 4 × 4 regular grid is shown in Fig. 4a-d. The black points represent the cell centers. The cell containing the particle centroid is identified as the starting cell (step 0) depicted with the black empty circle. The search algorithm is as follows:

• Step 1: Identification of the cells adjacent to the starting cell (Neumann neighborhood in red points).

• Step 2: Identification of the cells adjacent to the newly identified cells. The algorithm ignores a cell previously identified to avoid duplicates. In the 4 × 4 grid example (Fig. 4a-4d), we see that the 4 cells in the square lattice surrounding the starting cell (Moore neighborhood in red stars) are now identified.

• Step ≥ 3: Repeat the identification of the cells adjacent to the newly identified cells until the grid is mapped.

Note that, as this algorithm is mesh-independent, it can operate on unstructured or non-regular grids, including triangular, tetrahedral, or polyhedral cells without any supplementary steps. To find all the neighboring cells (Neumann and Moore points) of a starting cell using the peer-to-peer search, it takes only 2 steps in 2D and 3 steps in 3D (Fig. 4f) for both structured and unstructured grids.

To search the cells covered by particles, two rules are necessary to unmark newly identified cells: (i) the particle-to-cell-center distance is greater than the particle radius, (ii) the intersection point between the particle surface and the particle-to-cell-center line is outside the cell. Eventually, the algorithm has marked all the cells covered by particles including the particle-fluid interfaces as shown in Fig. 4e.

### 2.3.2. Solution strategy for solving VANS equations in CFD

In this part, we describe our solution strategy for solving the Volume-Averaged-Navier-Stokes equations using the Finite-Volume Method. The system has three unknown variables solved on a collocated grid, namely the void fraction, $ε$, the fluid pressure, $p$, and the fluid velocity, $v^f$. The void fraction results from the projection of the particle shadow onto the grid. The pressure-velocity coupling is solved using a semi-implicit time integration scheme adapted from the Pressure Implicit with Splitting of Operators (PISO) algorithm developed by Issa [51].
A known issue in CFD-DEM coupling is the presence of numerical instability on the velocity profile when the centroid of the particle is crossing a cell boundary [34,33]. A way to stabilize the particle-fluid coupling is to smooth $\epsilon$ after its mapping on the computational grid [32]. We use an isotropic diffusive smoothing controlled by the smoothing length $\lambda_s$. This is achieved using the parabolic filter:

$$ \frac{\partial \epsilon}{\partial t} = \nabla^2 \left( \frac{\lambda_s^2}{\Delta t_{CFD}} \epsilon \right), \qquad (17) $$

where $\Delta t_{CFD}$ is the time step used to solve Eqs. (2) and (3). A smoothing length of $\lambda_s \le \frac{D_p}{10^3}$ is found necessary to improve the stability of the simulations.

The pressure-velocity coupling is solved by forming a pressure equation from the continuity equation, Eq. (2), and a semi-discrete form of the momentum equation, Eq. (3). The latter is obtained using a forward Euler time integration between two successive times $n+1$ and $n$ for each cell. We obtain,

$$ \epsilon \frac{\mathbf{v}_C^{n+1} - \mathbf{v}_C^n}{\Delta t_{CFD}} = -a'_C \mathbf{v}_C^{n+1} + \sum_{NC} a'_{NC} \mathbf{v}_{NC}^{n+1} + \frac{\epsilon^2 \mu_f}{K} \mathbf{v}^p - \epsilon \nabla p', \qquad (18) $$

where the subscript C indicates the cell owner and NC the neighboring cells. The coefficients $a'_C$ and $a'_{NC}$ are the diagonal and off-diagonal coefficients of the space discretization of the momentum equation that includes advection and viscous dissipation effects. The pressure variable $p'$ divided by the fluid density, corresponds to the actual pressure including the hydrostatic pressure as $\nabla p' = \frac{1}{\rho_f}(-\rho_f \mathbf{g} + \nabla p)$.

The semi-discrete momentum equation can be recast into,

$$ a_C \mathbf{v}_C^{n+1} = H(\mathbf{v}) - \epsilon \nabla p', \qquad (19) $$

where $a_C = (\frac{\epsilon}{\Delta t_{CFD}} + a'_C)$ is the diagonal coefficients of the matrix for the velocity, and $H(\mathbf{v}) = \sum_{NC} a'_{NC} \mathbf{v}_{NC}^{n+1} + (\frac{\epsilon}{\Delta t_{CFD}}) \mathbf{v}_C^n + \frac{\epsilon^2 \mu_f}{K} \mathbf{v}^p$, contains the off-diagonal coefficients and the source terms.

Finally, the pressure equation is formed from the combination of the fluid mass balance equation (Eq. (2) divided by the fluid density), the parabolic filter (Eq. (17)), and the semi-discretized momentum (Eq. (19) divided by $a_C$),

$$ \nabla \cdot \left( \frac{\epsilon^2}{a_C} \nabla p' \right) = \nabla \cdot \left( \frac{\epsilon H(\mathbf{v})}{a_C} \right) + \nabla^2 \left( \frac{\lambda_s^2}{\Delta t_{CFD}} \epsilon \right). \qquad (20) $$

PISO is a predictor-corrector scheme for solving pressure-velocity coupling. Within a time step, the procedure is as follows. (i) Solve the discretized momentum equation (Eq. (18)) to compute a predicted velocity field, $\mathbf{v}^*$, for given boundary conditions and the pressure field, $p^n$ calculated at the previous time step. At this stage, the resulting guessed velocity does not satisfy the mass conservation. (ii) Solve the pressure equation (Eq. (20)) and guess the pressure field, $p^{**}$. (iii) Get the corrected velocity field using Eq. (19) and $p^{**}$. (iv) Update the boundary conditions. (v) Repeat the steps 2 to 5 for at least 2 iterations [51]. At the end of these steps, you get the velocity, $\mathbf{v}^{n+1}$, and pressure fields, $p^{n+1}$, for the next time-step.

### 2.3.3. Time-stepping stability criteria

Three numerical stability criteria are necessary for the fluid-particle coupling model proposed in this work: (i) a criterion for the pressure-velocity coupling algorithm, (ii) a constraint due to the calculation of colliding particles, and (iii) a stability criterion related to the fluid-particle interactions.

The PISO-like pressure-velocity algorithm for solving the VANS equations is not unconditionally stable and the time integration, $\Delta t_{CFD}$, is limited by a Courant-Friedrich-Lewy (CFL) condition [52,53],

$$ CFL = \Delta t_{CFD} \max \left( \frac{|\mathbf{v}^f|}{\Delta x} \right) < 1, \qquad (21) $$

where $\Delta x$ is the cell characteristic size.

The time-step requirements to capture the particle-particle and particle-wall collisions is the minimum of the particle relaxation time ($\tau_p = \frac{D_p^2 \rho_p}{18 \mu_f}$) and a portion of the Rayleigh time $T_{Ra}$ [54,55]. The latter corresponds to the time a shear wave takes to propagate through a solid particle [56]. The time-step criterion for DEM reads,

$$ \Delta t_{DEM} = \min(\alpha_t T_{Ra}, \tau_p) \quad \text{with} \quad T_{Ra} = \frac{\pi D_p}{2\Gamma} \sqrt{\frac{\rho_p}{G}}, \qquad (22) $$

where $\alpha_t$ is a positive constant lower than 1, $G = E/4(2-\nu)(1+\nu)$ is the particle shear modulus with $\rho_p$ the particle density, $E$ being the particle Young modulus, $\nu$ the particle Poisson ratio, and $\Gamma$ is a coefficient approximated by $\Gamma = 0.1631\nu + 0.8766$ [56]. As the constraints on the particle dynamics are stronger than on the fluid flow, $\Delta t_{DEM}$ is used as a sub-cycling within a CFD time-step to capture particle collision events accurately as shown in Fig. 5.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Temporal evolution for the fluid flow calculation compared to that of the particle dynamics calculation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Schematic representation of the flow past a stationary single particle, (b)-(d) CFD-DEM approach with different mesh resolution, and (a) CFD approach.</figcaption></figure>

The stability of the unresolved coupling results from the impact of the fluid-solid interaction on the DEM equations of each particle though the drag force is linked to the particle relaxation time. By assuming an explicit scheme on the source term integration, the coupling stability criterion is defined as:

$$ \Delta t_c \le \frac{4 D_p}{3 C_d} \frac{\rho_p}{\rho_f} \frac{1}{|\mathbf{v}^f - \mathbf{v}_i^p|} \qquad (23) $$

In practice, $\Delta t_{CFD}$ is taken as min($\Delta t_{CFD}$, $\Delta t_c$) and is, therefore, satisfying both Eq. (21) and Eq. (23).

### 3. Model verification

In this section, we present test cases to verify the robustness and efficiency of the unresolved-resolved four-way coupling CFD-DEM. First, we verify the implementation of the resolved-unresolved momentum exchange term (Section 3.1). Second, we assess the accuracy of our model in the case of a steel ball sedimenting in a water tube for which experimental data exists (Section 3.2). Then, we analyze the efficiency of our search algorithm for mapping the covered cells (Section 3.3).

#### 3.1. Assessment of the resolved momentum transfer calculation

The purpose of this section is to verify the implementation of the resolved drag force acting on both the particle (in DEM) and the fluid (in CFD). In resolved cases, the mutual drag between the fluid and the particle is not modeled by constitutive laws but it is an output of the simulation. Here, we run two sets of simulations to verify that: (i) in CFD, the fluid velocity profile around a resolved particle is accurate, (ii) in DEM, the drag applied to the particle is correct.

First, we consider a single stationary particle of diameter $D_p = 10$ µm located in the middle of a 50 µm × 50 µm square domain. The inlet is on the left-hand side and the outlet is on the right. A constant velocity $v_0 = 10^{-3}$ m/s ($Re_p = 0.01$) is applied at the inlet and the pressure is set to zero at the outlet. Lateral boundaries are set as walls with no-slip conditions. In CFD-DEM simulations, we use regular grids with different levels of refinement defined as the ratio of the particle diameter $D_p$ to the smallest cell size $\Delta x_{min}$ as illustrated in Fig. 6b-6d. The cells mapping the particle presence are penalized according to Eq. (4). A reference solution, denoted “CFD” hereafter, is obtained by solving the incompressible Navier-Stokes equations (using the so-called simpleFoam OpenFOAM solver) on a refined conformal grid ($\frac{D_p}{\Delta x_{min}} = 40$) in which no-slip condition is applied at the particle surface (see Fig. 6a). The velocity profile along the vertical axis that crosses the particle center is plotted in Fig. 7 for the x and y components. We observe a good agreement of the CFD-DEM with the reference CFD solution, especially for resolutions > 3.

Second, we focus on the calculated drag force acting on the particle surface. The simulation setup is similar to the first case, except that the square box is larger, corresponding to 80$D_p \times 80$$D_p$, to avoid boundary effects, and that the particle-fluid interface is resolved through a local mesh refinement as shown in Fig. 8. The level of refinement at the vicinity of the particle corresponds
<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) and (b), comparison of fluid velocity components $U_x$ and $U_y$ between results obtained from the model and a steady CFD solver of OpenFOAM (simpleFoam) calculated at the vertical line passing by the particle centroid.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Schematic representation of the domain, with (a)-(c) the local meshes refinement.</figcaption></figure>

to $\frac{D_p}{\Delta x_{min}} = 3, 5,$ and $7$, respectively. The inlet velocity, $v_0$, is set such that the Reynolds number varies from 0.001 to 100. In Fig. 9, we compare the resultant drag coefficient, $C_d$, with the empirical model proposed by Schiller and Naumann [57] which is used as a reference solution (the values come from Nguyen et al. [58]). It is calculated using $C_d = \frac{2|F_d|}{\rho_f v_0^2 A_p}$ where $F_d$ is the drag force acting on the particle in DEM calculated using Eq. (8), and $A_p$ is the particle surface projected on the Eulerian grid. We observe a very good agreement between the drag coefficient calculated with our resolved CFD-DEM model and the reference solution.

These simulations highlight the accuracy of the calculation of the drag force in our resolved CFD-DEM both from the CFD and the DEM sides. They also document the optimal mesh refinement in resolved CFD-DEM.

## 3.2. Sedimentation of a steel ball in a water tube

In this test case, we simulate the sedimentation of a steel ball in a water tube and compare the results with the experimental data of Allen [59]. The experiment consists of a 3.18 mm diameter steel (density of 7820 kg/m³) bead falling by gravity in a 11.5 cm long, 3 cm wide, and 28 cm high rectangular water tube. The fall of the particle initially placed at the top center of the domain without initial velocity is simulated using the new unresolved-resolved CFD-DEM coupling on different mesh resolutions ($\frac{D_p}{\Delta x} = \frac{1}{3}, \frac{1}{2}, 1, 2, 3$). The results are also compared with the prediction made by the original unresolved CFD-DEM solver of OpenFOAM (denseParticleFoam) for $\frac{D_p}{\Delta x} = 1$ (which corresponds to a 33 × 9 × 88 regular grid).

We see in Fig. 10 that the unresolved-resolved solver predictions are very close to the experimental value for both the resolved ($\frac{D_p}{\Delta x} \ge 1$) and unresolved ($\frac{D_p}{\Delta x} < 1$) cases. The terminal particle velocity calculated by the standard OpenFOAM unresolved CFD-DEM solver, however, does not match the reference data for $\frac{D_p}{\Delta x} = 1$. These simulations highlight the ability of the hybrid CFD-DEM model to capture accurately the particle trajectory regardless of the mesh resolution.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison between the drag coefficient obtained with the unresolved-resolved CFD-DEM model and the empirical data of Schiller and Naumann [57] for the flow around a stationary spherical particle.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Comparison of the falling velocity obtained with unresolved CFD-DEM coupling on $D_p/\Delta x = 1$, with the resolved-unresolved CFD-DEM coupling on different mesh resolutions.</figcaption></figure>

### 3.3. Searching algorithm efficiency analysis

In this section, we compare the time efficiency of our peer-to-peer search algorithm with the most frequently-used searching strategies, namely, the traverse search, the tree search (quadtree for 2D and octree for 3D) [60], and the linked-list search such as the Hilbert curve search [61,62]. In traditional traverse search, the computation time is devoted to traversing all the cells of the domain, which can become very important with a large number of cells and particles. Tree search continuously divides the current domain into 8 parts until it has no overlapping with the neighborhood area or it only contains one cell. Linked-list search organizes all the cells with the advantage of positioning spatially adjacent cells close to each other in the form of a linked list and then searching for it. The peer-to-peer search belongs to the family of linked-list search. Supposing that there are *N* cells and *N* particles in the domain, the time cost of traverse search, octree search, and Hilbert curve search are $O(N^2)$, $O(N \log N)$, and $O(N)$, respectively [36].
<figure><img src="image_14.png" /><figcaption>Fig. 11. Efficiency of the peer-to-peer search compared with the traverse, octree, and Hilbert curve search results taken from Wang et al. [36]. (a) Efficiency for an increasing number of particles in a case with 10000 cells. (b) Efficiency for an increasing number of cells in a case with 10000 particles.</figcaption></figure>

<figure><img src="image_15.png" /></figure>

The benchmark setup and data come from Wang et al. [36]. In this test case, the cells covered by particles are sought through traverse, octree, and Hilbert searching strategies. Particles of 1 cm diameter are placed randomly in a regular grid made of 1 cm × 1 cm × 1 cm cubic cells. Because the cell and particle center coordinates do not match, a particle can overlap up to 8 cells. Two test cases are considered. First, we investigate the efficiency of the algorithms for an increasing number of particles in a 10 cm × 10 cm × 100 cm grid (i.e. $10^4$ cells). Second, we compare their efficiency for a fixed number of $10^4$ particles within a domain whose size varies from 10 cm × 10 cm × 50 cm (i.e. $5 × 10^3$ cells) to 10 cm × 10 cm × 500 cm (i.e. $5 × 10^4$ cells). In each case, the simulation is run for 100 iterations.

We see in Fig. 11b and Fig. 11a that the peer-to-peer search efficiency is comparable to the Hilbert curve search when the number of cells increases (time cost of $O(N)$) and to the octree search when the particle number increases (time cost of $O(N \log N)$). In both cases, it is more efficient than the traverse search. One of the advantages of the peer-to-peer search is the easy implementation of the algorithm regardless of the structure of the mesh.

## 4. Application: pore-scale modeling of pore-clogging due to particulate flows

In this section, we use our unresolved-resolved four-way coupling CFD-DEM to simulate pore-clogging in porous media at the pore-scale. Section 4.1 investigates pore-clogging by sieving and bridging in a single pore. Section 4.2 discusses retention capacity and permeability reduction in a network of pores under different particle concentrations.

### 4.1. Clogging of a single pore

Sieving and bridging are two of the main pore-clogging mechanisms. In sieving, particles larger than the pore throat block at the pore entrance (exclusion of particles by size). In bridging, particles arrive simultaneously at the pore throat forming an arch. In this part, we simulate such processes in a single pore using our CFD-DEM package. The pore geometry is made of a two-dimensional converging-diverging channel of diameter, $D_c = 2$ mm, and a throat of diameter, $D_t = 500$ µm, gridded with a 150 × 50 conformal mesh. We consider two cases: (i) we inject three particles, one every three seconds, at the middle of the inlet. Particles have different diameters (100 µm, 300 µm, 550 µm) to simulate the sieving of the biggest particle ($D_p > D_t$). (ii) We continuously inject for 50 seconds, 2 particles per second, a polydisperse distribution of particles with sizes ranging from 100 to 300 µm to simulate the arch formation of particles with different sizes.

A constant pressure difference $\Delta p = 2$ µbar is applied between the inlet and the outlet, and the simulations are stopped when the pore is clogged. The fluid is water, the particles are in polystyrene, and the pore walls are in PDMS (Polydimethylsiloxane) (see properties in Table 4).

Snapshots of the two clogging mechanisms are presented in Figs. 12–13. In the first simulation (sieving, Fig. 12), we see that particles smaller than the pore throat ($D_p < D_t$) can pass through the constriction as expected, until a big particle ($D_p > D_t$) is filtered by its size. At that point, the pore permeability decreases and the flow velocity drops to a near-zero value as shown in Fig. 14. Notes that pore-clogging by size exclusion cannot be captured by the standard unresolved CFD-DEM solver of OpenFOAM because
<figure><img src="image_16.png" /><figcaption>Fig. 12. Sieving of a large particle at different times. The background color corresponds to the fluid velocity magnitude.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Bridging formation of 2 particles (214 µm and 287 µm) at different times. The background color corresponds to the fluid velocity magnitude.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Evolution of the permeability damage in the case of sieving and bridging.</figcaption></figure>

the size of particles is larger than the cell size ($D_p > \Delta x$), and the particle trajectories become unphysical. Our resolved-unresolved approach, however, is not constrained by such limitations and can be used in all configurations. In the second scenario (bridging, Fig. 13), the suspension flows in the domain until two particles arrive simultaneously at the pore entrance, forming an arch, and blocking the flow. The time-to-plug by arch formation depends on the probability that such an event happens, the closer the ratio $\frac{D_p}{D_i}$ is to 1, the higher the probability of bridging [20]. Pore-scale modeling of particulate flow with CFD-DEM will bring new insights into the assessment of the time-to-plug according to flow conditions and suspension properties.

## 4.2. Clogging in a pore-network

In this part, we use our CFD-DEM package to investigate particle retention and permeability reduction in a porous medium that represents a PDMS (Polydimethylsiloxane) microfluidic device [63]. The pore geometry and the meshing procedure are found in Soulaine [64]. We obtain a two-dimensional 1050 µm × 310 µm pore-scale domain is shown in Fig. 15. The typical pore size is $d_{50} = 26$ µm, porosity is $\epsilon_i = 0.62$, and permeability is $K_i = 1.1 \times 10^{-11} m^2$. The domain is discretized into an unstructured mesh with 49330 cells using `snappyHexMesh`, the OpenFOAM automatic gridder. The left and right sides are inlet and outlet boundaries, respectively. The top, bottom, and grain surfaces are walls described with no-slip conditions. We consider a cloud of $D_p = 6$ µm diameter monodispersed polystyrene particles. The mean mesh resolution is $D_p/\Delta x \approx 5$, which allows for an accurate calculation of the drag forces (see Section 3.1). Fluid, particle, and PDMS properties as well as simulation parameters are listed in Table 4.
<figure><img src="image_19.png" /><figcaption>Fig. 15. (a) Geometry of the 2D simulation model of the porous media.</figcaption></figure>

<table><caption>**Table 4**

Parameters of the CFD-DEM simulations.</caption><thead><tr><th>CFD and fluid parameters</th><th></th><th>DEM and particle parameters</th><th></th></tr></thead><tbody><tr><td>CFD time-step Δt<sub>CFD</sub></td><td>2 × 10<sup>-5</sup> s</td><td>DEM time-step Δt<sub>DEM</sub></td><td>5 × 10<sup>-9</sup> s</td></tr><tr><td>Fluid density ρ<sub>f</sub></td><td>10<sup>3</sup> kg/m<sup>3</sup></td><td>Particle density ρ<sub>p</sub></td><td>1050 kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity μ<sub>f</sub></td><td>10<sup>-3</sup> Pa.s</td><td>Particle concentration C<sub>0</sub> (mass percentage)</td><td>0.05% ; 0.1% ; 0.15%</td></tr><tr><td>Model depth h</td><td>20 μm</td><td>Poisson ratio ν<sub>p</sub> / ν<sub>w</sub></td><td>0.34/0.5</td></tr><tr><td></td><td></td><td>Young modulus E<sub>p</sub> / E<sub>w</sub></td><td>3 GPa / 2 MPa</td></tr><tr><td></td><td></td><td>Friction coefficient μ<sub>p</sub> / μ<sub>w</sub></td><td>0.4 / 0.84</td></tr></tbody></table>

PV: 2.808400

<figure><img src="image_20.png" /></figure>

PV: 6.820400

<figure><img src="image_21.png" /><figcaption>Fig. 16. Particle deposition and velocity magnitude field for $C_0 = 0.1\%$ during the injection of particles ($PV = 2.8$) and post-injection ($PV = 6.8$). The orange spheres are particles and the color map represents the fluid velocity magnitude. The clogs formed during the particle injection and remain stable when the particle injection is stopped.</figcaption></figure>

Simulations are run in 3 consecutive steps. First, a constant pressure difference, $\Delta p = 0.2$ mbar, is applied between the inlet and the outlet. The fluid (water) flows through the porous formation until the flow field is stable. Pore volume (PV) – a measure of time defined as the ratio of the injected volume of water to the pore-space volume – is set to $PV = 0$. Then, a concentration of particles, $C_0$, is continuously injected from the inlet until $PV = 4$. Finally, the injection of particles is stopped, and the simulation is run until steady-state. Simulations use a Hele-Shaw correction term, $12\mu_f c v^f / h^2$, in Eq. (3) to account for the hydrodynamic effects in the thickness, $h$, of the microfluidic device [65].

Snapshots of the particle migration and retention are shown in Fig. 16 ($C_0 = 0.1\%$) and Fig. 17 ($C_0 = 0.15\%$). We observe that some particles percolate while others remain trapped within the porous medium clogging pores. The retention of particles changes
<figure><img src="image_22.png" /><figcaption>**Fig. 17.** Particle deposition and velocity magnitude field for $C_0 = 0.15\%$ during the injection of particles ($PV = 2.88$) and post-injection ($PV = 6.9$). Arches of particles are formed during the injection, but they are unstable, and most of them are remobilized.</figcaption></figure>

the system porosity and reroutes the local flow lines which affect the system permeability. The particle contact frequency and the probability of simultaneous arrival of particles at a pore entrance to form a clog (arch formation) increase with the particle concentration. For $C_0 = 0.05\%$ (results not shown here), there is no clog inside the porous medium – although, a few isolated particles are trapped in dead-end pores – which is not the case for higher concentrations. For $C_0 = 0.1\%$ (Fig. 16), we see several clogs formed by the bridging of 2 particles at the entrance of small pores and others that grow until they reach their filling limit. The clogs are very stable and remain even after we stopped the injection of particles ($PV > 4$). For $C_0 = 0.15\%$, the possibility of arch formation increases because the particles are more likely to be near each other as they reach a pore entrance due to converging trajectories. We observe (Fig. 17), however, several non-permanent pores clogged. These observations are also seen in the particle cumulative breakthrough curves (Fig. 18a) and the plot of the number of clogs over time. Once the injection stops, we see that almost all particles go through the porous structure for $C_0 = 0.05\%$, and up to 6% of particles remain trapped for $C_0 = 0.1\%$, and 3% for $C_0 = 0.15\%$. For the highest concentration, more clogs are formed but they are less stable and remobilized, eventually (see Fig. 18b).

Pore-clogging leads to permeability reduction. The evolution of the porosity damage severity, $e_r = \langle\epsilon\rangle/\epsilon_i$, with $\langle\epsilon\rangle$ the domain-averaged porosity, and the permeability damage severity, $K_r = K_f/K_i$, is shown in Fig. 19a. The permeability, $K_f$, of the porous system is obtained using Darcy's law, $K_f = \frac{\langle v_x^f \rangle \mu_r L}{\Delta p}$, where $\langle v_x^f \rangle$ is the domain-averaged fluid velocity, $L$ is the length of the domain, and $\Delta p$ is the pressure difference. For identical hydrodynamic conditions, permeability reduction increases with the particle concentration in agreement with experimental observations [66,67,14]. A permeability reduction of 10% is also observed for $C_0 = 0.05\%$ for which there is no pore-clogging. This reduction corresponds to the flow resistance related to the suspension itself. The permeability recovers its initial value after all the particles have been flushed out. For $C_0 = 0.1\%$ and $C_0 = 0.15\%$, we observe a hysteresis between the initial and final permeability/porosity values. It is related to the clogs formed with the retention of particles.

The Probability Density Function (PDF) presented in Fig. 19b gives information about the redistribution of local flow rates during pore-clogging. For $C_0 = 0.05\%$, it is superimposed with the PDF without particles. Indeed, as almost all the particles percolate and the remaining particles are trapped individually in dead-end pores, they do not affect the fluid flow. For $C_0 = 0.1\%$ and $0.15\%$, we observe an increase in the densities of low fluid velocities and a reduction of the highest value of the fluid velocities. This is characteristic of pore-clogging [68]. The multiple permanent clogs formed with $C_0 = 0.1\%$ create additional dead-end pores and subsequent fluid recirculations that increase negative values of the velocity fluid.

These simulations illustrate the potential of the resolved-unresolved four-way coupling CFD-DEM to investigate the complex interplay between particle migration, deposition, and remobilization in porous media.

## 5. Conclusions

Clogging of porous media by particulate flow is a complex process that relies on many different physical phenomena and depends on a wide range of parameters including particle size and concentration, flow conditions, and the geometry of porous microstructure. We proposed an unresolved-resolved four-way coupling CFD-DEM to simulate particle migration, deposition, and retention in porous media at the pore-scale, and therefore to investigate locally the permeability reduction due to the injection of particles.
<figure><img src="image_23.png" /><figcaption>**Fig. 18.** (a) Particle cumulative BTC (Breakthrough Curves) normalized by the total number of injected particles. (b) Evolution of the current number of clogs (Current Nclog) in the system, and the number of clogs formed at different positions (Different Nclog) for different $C_0$. For $C_0 = 0.1\%$, each clog formed at a unique position. The pink zone represents the particle injection period.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 19.** (a) Permeability damage severity, $K_r$, and porosity damage severity, $\epsilon_r$, for different particle concentration, $C_0$. (b) The probability density function of the fluid longitudinal velocity $v_x'$ for different $C_0$ at the end of the simulation compared with the case without particles.</figcaption></figure>

Unlike other CFD-DEM, our approach is independent of the grid resolution and type (i.e. structured and unstructured) and is not limited to particles smaller than the cell size. It relies on a hybrid resolved-unresolved formulation of the drag force calculation and an efficient searching strategy – called peer-to-peer search – to identify the CFD cells covered by DEM particles. Moreover, a diffusive smoothing filter removes unphysical oscillations when particles cross the cell-to-cell interface. The method has been successfully validated using cases for which reference solutions exist. The peer-to-peer search is as efficient as the Hilbert curve search and easier to implement. This searching algorithm as well as the unresolved-resolved CFD-DEM can operate on unstructured or non-regular grids, including triangular, tetrahedral, or polyhedral meshes which are used when refining the grid near complex boundaries.
Our method has great potential to investigate the complex feedback due to the injection, retention, and remobilization of particles in porous media. Our case studies show the ability of the resolved-unresolved four-way coupling CFD-DEM to capture two of the pore-clogging mechanisms, namely clogging by size exclusion (sieving) and by arch formation (bridging, due to the simultaneous arrival of particles at a pore entrance). The physics of the third clogging mechanism by aggregation of particles relies on electrochemical interactions that act on the particle surfaces. The numerical implementation of these forces in the CFD-DEM framework is discussed in Maya et al. [69].

Nomenclature

$\alpha$ Constant related to elasticity coefficient

$\alpha_t$ Constant related to DEM time-step

$\beta$ interphase momentum exchange coefficient

$\delta_n$ Normal overlap

$\delta_t$ Tangential overlap

$\omega$ Particle angular velocity

$\Delta p$ Fluid pressure difference

$\Delta t$ Time-step

$\Delta x$ Cell characteristic length

$\epsilon$ Local cell porosity or void fraction

$\epsilon_r$ Porosity damage severity

$\epsilon_{min}$ Minimum porosity value

$\eta$ Particle damping coefficient

$\eta_n$ Particle normal damping coefficient

$\eta_t$ Particle tangential damping coefficient

$\Gamma$ Coefficient related to DEM time-step

$\lambda_s$ Smoothing length

$\langle\epsilon\rangle$ Domain-averaged porosity

$\langle v_x^f \rangle$ Domain-averaged fluid velocity

$\mu$ Particle friction coefficient

$\mu_f$ Fluid dynamic viscosity

$\nu$ Particle Poisson ratio

$\bar{v}_p$ Averaged particle Velocity

$\rho_f$ Fluid density

$\rho_p$ Particle density

$\tau$ Fluid shear-rate

$\tau_p$ Particle relaxation time

$\mathbf{F}^c$ Particle contact force

$\mathbf{F}^f$ Particle-fluid force

$\mathbf{F}^g$ Particle gravitational force

$\mathbf{F}^c_n$ Particle normal contact force

$\mathbf{F}^c_t$ Particle tangential contact force

$\mathbf{G}$ Particle shear modulus

$\mathbf{g}$ Gravitational acceleration

$\mathbf{M}^c$ Particle contact Torque

$\mathbf{M}^{hyd}$ Particle hydrodynamic Torque

$\mathbf{n}$ unit vector

$\mathbf{p}$ Position vector pointing particle center

$\mathbf{s}$ Position vector pointing cell center

$\mathbf{U}^p$ Particle translational velocity

$\mathbf{v}^*$ Predicted fluid velocity

$\mathbf{v}^f$ Fluid velocity

$\mathbf{v}^p$ Particle velocity

$\mathbf{v}_c^f$ Fluid velocity of cell owner at time n

$\mathbf{v}_{ij}^p$ Particles relative velocity

$\mathbf{v}_t^p$ Particle tangential slip velocity

$a_C$ Diagonal coefficients of matrix for fluid velocity

$a'_C$ Diagonal coefficients of momentum equation sapce discretization

$A_p$ Particle area

$a'_{NC}$ Off-diagonal coefficients of momentum equation space discretization

$C_0$ Affluent particle concentration

$C_d$ Drag coefficient

$C_f$ Effluent particle concentration

$C_r$ Remaining particle concentration

$CFL$ Courant number

$D_c$ Channel diameter

$D_p$ Particle diameter

$D_t$ Throat diameter

$d_{50}$ Mean diameter pore size

$E$ Particle Young modulus

$F_d$ Drag force

$H(v)$ Off-diagonal coefficients and source terms

$h$ Geometry depth

$I$ Particle moment of inertia

$K$ Cell-permeability

$k$ Particle stiffness

$K_0$ Low given permeability

$K_f$ Final domain permeability

$K_i$ Initial domain permeability

$k_n$ Particle normal stiffness

$K_r$ Permeability damage severity

$k_t$ Particle tangential stiffness

$L$ Geometry length

$m$ Particle mass

$n^c$ Number of contacting objects

$n^k$ Number of covered cells

$p$ Fluid pressure

$p'$ Fluid pressure including hydrostatic pressure

$p^{**}$ PISO calculated fluid pressure

$PV$ Pore volume

$R$ Particle radius

$r_{ij}$ Equivalent particle radius

$Re_p$ Particle Reynolds number

$t$ Time

$T_{Ra}$ Rayleigh time

$v_0$ Inlet fluid velocity

$V_c$ Cell volume

$V_p$ Particle volume

$v_x^f$ Longitudinal fluid velocity

CRediT authorship contribution statement

**Laurez Maya Fogouang:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Laurent André:** Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. **Cyprien Soulaine:** Writing – review & editing, Writing – original draft,
Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Cyprien Soulaine reports financial support was provided by European Research Council. Laurez Maya Fogouang reports financial support was provided by Bureau for Geological and Mining Research. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This research is funded by the European Union (ERC, COCONUT, 101043288). LM received a PhD fellowship from BRGM and Région Centre-Val de Loire. The authors benefitted from the Centre de Calcul Scientifique en région Centre-Val de Loire's (CaSciMoDOT) high-performance computational resources. The authors would like to thank Laurent Oxarango, Olivier Liot, and Noura Eddaoui for their interesting discussions and suggestions.

## Data availability

Data will be made available on request.

## References

[1] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle-fluid suspensions, J. Stat. Phys. 104 (2001) 1191-1251, https://doi.org/10.1023/a:1010414013942.

[2] D.J. Rader, A.S. Geller, Transport and deposition of aerosol particles, in: Developments in Surface Contamination and Cleaning, Elsevier, 2008, pp. 23-90.

[3] A. Hager, CFD-DEM on Multiple Scales: An Extensive Investigation of Particle-Fluid Interactions, Ph.D. thesis, Johannes Kepler Universitat Linz, 2014.

[4] L. Oxarango, Heat and mass transfer in multi-scale porous structures: application to diesel particulate filters modelling, 2004.

[5] A.J. Torregrosa, J.R. Serrano, P. Piqueras, Ó. García-Afonso, Experimental and computational approach to the transient behaviour of wall-flow diesel particulate filters, Energy 119 (2017) 887-900, https://doi.org/10.1016/j.energy.2016.11.051.

[6] O.V. Soloveva, S.A. Solovev, R.R. Yafizov, Numerical simulation of the aerosol particle motion in granular filters with solid and porous granules, Processes 9 (2021) 268, https://doi.org/10.3390/pr9020268.

[7] C. Fetter Jr, R. Holzmacher, Groundwater recharge with treated wastewater, J. Water Pollut. Control Fed. (1974) 260-270.

[8] N. Saleh, T. Phenrat, K. Sirk, B. Dufour, J. Ok, T. Sarbu, K. Matyjaszewski, R.D. Tilton, G.V. Lowry, Adsorbed triblock copolymers deliver reactive iron nanoparticles to the oil/water interface, Nano Lett. 5 (2005) 2489-2494, https://doi.org/10.1021/nl0518268.

[9] H. Bouwer, Artificial recharge of groundwater: hydrogeology and engineering, Hydrogeol. J. 10 (2002) 121-142, https://doi.org/10.1007/s10040-001-0182-4.

[10] W. Song, X. Liu, T. Zheng, J. Yang, A review of recharge and clogging in sandstone aquifer, Geothermics 87 (2020) 101857, https://doi.org/10.1016/j.geothermics.2020.101857.

[11] Q. Xie, A. Saeedi, C.D. Piane, L. Esteban, P.V. Brady, Fines migration during CO2 injection: experimental results interpreted using surface forces, Int. J. Greenh. Gas Control 65 (2017) 32-39, https://doi.org/10.1016/j.ijggc.2017.08.011.

[12] N. Eddaoui, M. Panfilov, L. Ganzer, B. Hagemann, Impact of pore clogging by bacteria on underground hydrogen storage, Transp. Porous Media 139 (2021) 89-108, https://doi.org/10.1007/s11242-021-01647-6.

[13] F. Osselin, C. Soulaine, C. Fauguerolles, E.C. Gaucher, B. Scaillet, M. Pichavant, Orange hydrogen is the new green, Nat. Geosci. 15 (2022) 765-769, https://doi.org/10.1038/s41561-022-01043-9.

[14] G.C. Agbangla, É. Climent, P. Bacchin, Experimental investigation of pore clogging by microparticles: evidence for a critical flux density of particle yielding arches and deposits, Sep. Purif. Technol. 101 (2012) 42-48, https://doi.org/10.1016/j.seppur.2012.09.011.

[15] N. Debnath, A. Kumar, T. Thundat, M. Sadrzadeh, Investigating fouling at the pore-scale using a microfluidic membrane mimic filtration system, Sci. Rep. 9 (2019), https://doi.org/10.1038/s41598-019-47096-6.

[16] M. Masoudi, H. Fazeli, R. Miri, H. Hellevang, Pore scale modeling and evaluation of clogging behavior of salt crystal aggregates in CO2-rich phase during carbon storage, Int. J. Greenh. Gas Control 111 (2021) 103475, https://doi.org/10.1016/j.ijggc.2021.103475.

[17] Y. Song, E. Shim, Structure characterization of the clogging process of coarse fibrous filter media during solid particle loading with x-ray micro-computed tomography, Sep. Purif. Technol. 273 (2021) 118980, https://doi.org/10.1016/j.seppur.2021.118980.

[18] R. Ines, K. Eric, F. Marc, P. Nicolas, L. Béatrice, Formation damage induced by clay colloids deposition in Triassic clastic geothermal fields: insights from an experimental approach, Appl. Clay Sci. 234 (2023) 106868, https://doi.org/10.1016/j.clay.2023.106868.

[19] H.M. Wyss, D.L. Blair, J.F. Morris, H.A. Stone, D.A. Weitz, Mechanism for clogging of microchannels, Phys. Rev. E 74 (2006) 061402, https://doi.org/10.1103/PhysRevE.74.061402, https://link.aps.org/doi/10.1103/PhysRevE.74.061402.

[20] E. Dressaire, A. Sauret, Clogging of microfluidic systems, Soft Matter 13 (2017) 37-48, https://doi.org/10.1039/c6sm01879c.

[21] W. Guo, H.-H. Ngo, J. Li, A mini-review on membrane fouling, Bioresour. Technol. 122 (2012) 27-34, https://doi.org/10.1016/j.biortech.2012.04.089.

[22] T. Wang, F. Zhang, J. Furtney, B. Damjanac, A review of methods, applications and limitations for incorporating fluid flow in the discrete element method, J. Rock Mech. Geotech. Eng. 14 (2022) 1005-1024, https://doi.org/10.1016/j.jrmge.2021.10.015.

[23] J. Dennis, Properties of stationary (bubbling) fluidised beds relevant to combustion and gasification systems, 2013, pp. 77-148e, https://doi.org/10.1533/9780857098801.1.77.

[24] Q. Li, V. PrigioBBbe, Numerical simulations of the migration of fine particles through porous media, Transp. Porous Media 122 (2018) 745-759, https://doi.org/10.1007/s11242-018-1024-3.

[25] K. Zhou, J. Hou, Q. Sun, L. Guo, S. Bing, Q. Du, C. Yao, A study on particle suspension flow and permeability impairment in porous media using lbm-dem-imb simulation method, Transp. Porous Media 124 (2018) 681-698, https://doi.org/10.1007/s11242-018-1089-z.

[26] Y. Tsuji, T. Kawaguchi, T. Tanaka, Discrete particle simulation of two-dimensional fluidized bed, Powder Technol. 77 (1993) 79-87, https://doi.org/10.1016/0032-5910(93)85010-7.
[27] C. O'Sullivan, Particulate Discrete Element Modelling: A Geomechanics Perspective, Applied Geotechnics, vol. 4, Spon Press, 2011.

[28] C. Rettinger, U. Rüde, A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann simulations, Comput. Fluids 154 (2017) 74–89,
https://doi.org/10.1016/j.compfluid.2017.05.033.

[29] J. Su, G. Chai, L. Wang, W. Cao, Z. Gu, C. Chen, X.Y. Xu, Pore-scale direct numerical simulation of particle transport in porous media, Chem. Eng. Sci. 199 (2019) 613–627, https://doi.org/10.1016/j.ces.2019.01.033.

[30] A. Elrahmani, R.I. Al-Raoush, H. Abugazia, T. Seers, Pore-scale simulation of fine particles migration in porous media using coupled cfd-dem, Powder Technol. 398 (2022) 117130, https://doi.org/10.1016/j.powtec.2022.117130.

[31] X. Nan, J. Hou, G. Li, Z. Shen, W. Wen, D. Wei, Study of particle clogging pattern in concrete seepage process based on multi-scale vof-dem and experimental comparison, Constr. Build. Mater. 347 (2022) 128496, https://doi.org/10.1016/j.conbuildmat.2022.128496.

[32] S. Pirker, D. Kahrimanovic, C. Goniva, Improving the applicability of discrete phase simulations by smoothening their exchange fields, Appl. Math. Model. 35 (2011) 2479–2488, https://doi.org/10.1016/j.apm.2010.11.066.

[33] J. Marshall, K. Sala, Comparison of methods for computing the concentration field of a particulate flow, Int. J. Multiph. Flow 56 (2013) 4–14, https://doi.org/10.1016/j.jmultiphaseflow.2013.05.009.

[34] Z. Peng, E. Doroodchi, C. Luo, B. Moghtaderi, Influence of void fraction calculation on fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds, AIChE J. 60 (2014) 2000–2018, https://doi.org/10.1002/aic.14421.

[35] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448–476, https://doi.org/10.1016/j.jcp.2005.03.017.

[36] Z. Wang, Y. Teng, M. Liu, A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and Hilbert curve based searching strategy, J. Comput. Phys. 384 (2019) 151–169, https://doi.org/10.1016/j.jcp.2019.01.017.

[37] S.T.W. Kuruneru, E. Marechal, M. Deligant, S. Khelladi, F. Ravelet, S.C. Saha, E. Sauret, Y. Gu, A comparative study of mixed resolved-unresolved CFD-DEM and unresolved CFD-DEM methods for the solution of particle-laden liquid flows, Arch. Comput. Methods Eng. 26 (2018) 1239–1254, https://doi.org/10.1007/s11831-018-9282-3.

[38] C. Soulaine, S. Roman, A. Kovscek, H.A. Tchelepi, Mineral dissolution and wormholing from a pore-scale perspective, J. Fluid Mech. 827 (2017) 457–483, https://doi.org/10.1017/jfm.2017.499.

[39] Z.Y. Zhou, S.B. Kuang, K.W. Chu, A.B. Yu, Discrete particle simulation of particle-fluid flow: model formulations and their applicability, J. Fluid Mech. 661 (2010) 482–510, https://doi.org/10.1017/s002211201000306x.

[40] S. Whitaker, The Forchheimer equation: a theoretical development, Transp. Porous Media 25 (1996) 27–61, https://doi.org/10.1007/bf00141261.

[41] J.-L. Auriault, On the domain of validity of Brinkman's equation, Transp. Porous Media 79 (2009) 215–223, https://doi.org/10.1007/s11242-008-9308-7.

[42] P. Angot, C.-H. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numer. Math. 81 (1999) 497–520,
https://doi.org/10.1007/s002110050401.

[43] C. Soulaine, H.A. Tchelepi, Micro-continuum approach for pore-scale simulation of subsurface processes, Transp. Porous Media 113 (2016) 431–456, https://doi.org/10.1007/s11242-016-0701-3.

[44] J.M. Minto, R.J. Lunn, G.E. Mountassir, Development of a reactive transport model for field-scale simulation of microbially induced carbonate precipitation, Water Resour. Res. 55 (2019) 7229–7245, https://doi.org/10.1029/2019wr025153.

[45] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65, https://doi.org/10.1680/geot.1979.29.1.47.

[46] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94.

[47] C. Wen, Y. Yu, Mechanics of Fluidization, The Chemical Engineering Progress Symposium Series, 1966.

[48] H. Hertz, Uber die berührung fester elastischer korper (on the contact of elastic solids), J. Reine Angew. Math. (1882).

[49] N. Deen, M.V.S. Annaland, M.V. der Hoef, J. Kuipers, Review of discrete particle modeling of fluidized beds, Chem. Eng. Sci. 62 (2007) 28–44, https://doi.org/10.1016/j.ces.2006.08.014.

[50] S. Torkzaban, S.A. Bradford, S.L. Walker, Resolving the coupled effects of hydrodynamics and DLVO forces on colloid attachment in porous media, Langmuir 23 (2007) 9652–9660, https://doi.org/10.1021/la700995e.

[51] R. Issa, Solution of the implicitly discretized fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1986) 40–65, https://doi.org/10.1016/0021-9991(86)90099-9.

[52] R. Courant, K. Friedrichs, H. Lewy, Uber die Partiellen Differenzengleichungen der Mathematischen Physik, Math. Ann. 100 (1928) 32–74, https://doi.org/10.1007/bf01448839.

[53] J.H. Ferziger, M. Perić, R.L. Street, Computational Methods for Fluid Dynamics, Springer International Publishing, 2020.

[54] C. Thornton, C. Randall, Applications of theoretical contact mechanics to solid particle system simulation, Stud. Appl. Mech. 20 (1988) 133–142.

[55] S.J. Burns, P.T. Piirroinen, K.J. Hanley, Critical time step for dem simulations of dynamic systems using a hertzian contact model, Int. J. Numer. Methods Eng. 119 (2019) 432–451, https://doi.org/10.1002/nme.6056, https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6056.

[56] Y. Li, Y. Xu, C. Thornton, A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles, Powder Technol. 160 (2005) 219–228, https://doi.org/10.1016/j.powtec.2005.09.002.

[57] L. Schiller, A. Naumann, A drag coefficient correlation, Z. Ver. Dtsch. Ing. (1935).

[58] G.T. Nguyen, E.L. Chan, T. Tsuji, T. Tanaka, K. Washino, Resolved CFD-DEM coupling simulation using volume penalisation method, Adv. Powder Technol. 32 (2021) 225–236, https://doi.org/10.1016/j.apt.2020.12.004.

[59] H. Allen, The motion of a sphere in a viscous fluid, Lond. Edinb. Dublin Philos. Mag. J. Sci. 50 (1900) 519–534, https://doi.org/10.1080/14786440009463941.

[60] L. Hernquist, N. Katz, TREESPH - a unification of SPH with the hierarchical tree method, Astrophys. J. Suppl. Ser. 70 (1989) 419, https://doi.org/10.1086/191344.

[61] J. Monaghan, Particle methods for hydrodynamics, Comput. Phys. Rep. 3 (1985) 71–124, https://doi.org/10.1016/0167-7977(85)90010-3.

[62] H. Jagadish, Analysis of the Hilbert curve for representing two-dimensional space, Inf. Process. Lett. 62 (1997) 17–22, https://doi.org/10.1016/s0020-0190(97)00014-8.

[63] S. Roman, C. Soulaine, M.A. AlSaud, A. Kovscek, H. Tchelepi, Particle velocimetry analysis of immiscible two-phase flow in micromodels, Adv. Water Resour. 95 (2016) 199–211, https://doi.org/10.1016/j.advwatres.2015.08.015.

[64] C. Soulaine, Introduction to open-source computational fluid dynamics using openfoam® technology - part I, https://doi.org/10.5281/ZENODO.10647983, 2024.

[65] C. Soulaine, J. Maes, S. Roman, Computational microfluidics for geosciences, Front. Water 3 (2021), https://doi.org/10.3389/frwa.2021.643714.

[66] V. Ramachandran, H.S. Fogler, Plugging by hydrodynamic bridging during flow of stable colloidal particles within cylindrical pores, J. Fluid Mech. 385 (1999) 129–156, https://doi.org/10.1017/s0022112098004121.

[67] B. Mustin, B. Stoeber, Deposition of particles from polydisperse suspensions in microfluidic systems, Microfluid. Nanofluid. 9 (2010) 905–913, https://doi.org/10.1007/s10404-010-0613-4.

[68] A. Velásquez-Parra, T. Aquino, M. Willmann, Y. Méheust, T. Le Borgne, J. Jiménez-Martínez, Sharp transition to strongly anomalous transport in unsaturated porous media, Geophys. Res. Lett. 49 (2022), https://doi.org/10.1029/2021gl096280.

[69] L. Maya Fogouang, L. André, P. Leroy, C. Soulaine, Particulate transport in porous media at pore-scale. Part 2: CFD-DEM and colloidal forces, J. Comput. Phys. (2024), https://doi.org/10.1016/j.jcp.2024.113439.