Aggregation and clogging phenomena of rigid microparticles in microfluidics

We developed a numerical tool to investigate the phenomena of aggregation and clogging of rigid microparticles suspended in a Newtonian fluid transported through a straight microchannel. In a first step, we implement a time-dependent one-way coupling Discrete Element Method (DEM) technique to simulate the movement and effect of adhesion on rigid microparticles in two- and three-dimensional computational domains. The Johnson–Kendall–Roberts (JKR) theory of adhesion is applied to investigate the contact mechanics of particle–particle and particle–wall interactions. Using the one-way coupled solver, the agglomeration, aggregation and deposition behavior of the microparticles is studied by varying the Reynolds number and the particle adhesion. In a second step, we apply a two-way coupling CFD–DEM approach, which solves the equation of motion for each particle, and transfers the force field corresponding to particle–fluid interactions to the CFD toolbox OpenFOAM. Results for the one-way (DEM) and two-way (CFD–DEM) coupling techniques are compared in terms of aggregate size, aggregate percentages, spatial and temporal evaluation of aggregates in 2D and 3D. We conclude that two-way coupling is the more realistic approach, which can accurately capture the particle–fluid dynamics in microfluidic applications. Electronic supplementary material The online version of this article (10.1007/s10404-018-2124-7) contains supplementary material, which is available to authorized users.


S1 Discrete Element Method (DEM): One-way coupling approach
In the following we briefly discuss the force terms used in Eq. (1).

S1.1 Fluid forces
The fluid force experience by each particle in laminar flow in a Newtonian fluid is given by the drag force F D : where, u f and µ is the fluid velocity and dynamic viscosity. We consider fully developed laminar flow, and the velocity profile for 2D channel flow can be calculated as and for a 3D cylindrical channel as

S1.2 Collision and adhesion forces
We neglect sliding and twisting motions of interacting particles, and therefore only the normal component of adhesion and collision forces F A in Eq. (1) are considered, where F n represents the normal force magnitude and n is the unit vector connecting the centers of two particles: where x i and x j are the position vectors to two generic particles i and j. The normal force magnitude F n is the sum of adhesion and elastic deformation of the particles and energy losses during the impact of the particles. The latter contribution can be neglected for very small particles [4].
During the computations we also make use of the effective radius R and the effective elastic modulus E of two generic particles i and j upon impact, which are defined as where r i and r j are the radii of two particles with σ i and σ j their Poisson ratios and E i and E j their elastic moduli. We consider the Johnson-Kendall-Roberts (JKR) contact model [2] to capture the aggregate generation in the flowing suspension. We assume that particle-particle collisions are limited to their contact region, and therefore the contact upon particle collision is computed only by the instantaneous forces occurring at that time but not by their previous history. The normal elastic force magnitude F n and particle normal overlap δ N can be calculated as a function of the contact region radius a [1,2]: where F C is the critical force which can be calculated using the Van der Waals surface potential energy of the particle γ The contact region radius a is linked to the particle normal overlap δ N via and by the following relation between δ N and critical overlap δ C in which the term δ C provides the critical overlap distance between two particles The equilibrium contact radius a 0 used in above Eqs. (S8), (S11) and (S12) is the contact region suggested in the JKR model [2] and is expressed as: To speed-up the simulations pre-computation was performed for Eqs. (S8) and (S11). The results are stored in a look-up table, i.e., both δ N /δ C and F N /F C are solved for variuos a/a 0 values. During the time-stepping computation the generated δ N /δ C value is then used (with appropriate interpolation) to extract the corresponding value of F N /F C from the table. We also adopt a neighbor-list algorithm approach, in such a way that particles can experience only collision/adhesion forces at a certain distance from each other. For this, the entire channel is divided into small blocks of size d + δ C , and the collision/adhesion force computation is only performed if the particles are present in three consecutive blocks. Table S1: Equations to calculate the forces and parameters in the DEM part of the CFD-DEM coupling. The indices i and j denote two particles.

S3 CFD-DEM coupling strategy
To apply the CFD-DEM coupling, both the DEM solver (LIGGGHTS) and an OpenFOAM CFD solver are run in parallel. The data exchange between LIGGGHTS and OpenFOAM is done after a predefined number of time-steps (here: 10 DEM time-steps). The coupling approach is based on several computation steps [3]: 1. The particle positions and velocities are calculated by the DEM solver (LIGGGHTS) 2. The particle positions and velocities are transferred to the CFD solver (OpenFOAM) 3. The corresponding cell according to the provided position of each particle in the CFD mesh is found 4. The particle volume fraction and mean particle velocity is calculated for each cell 5. The fluid forces acting on each particle are calculated according to the particle volume fraction in that cell 6. The particle-fluid momentum exchange terms are computed by averaging all the particle based forces in each CFD cell 7. The fluid velocity is computed by using the local volume fraction and momentum exchange terms 8. The data are returned to the DEM solver for the next time step after calculating the fluid forces on each particle 9. The entire routine is repeated starting from Eq. (3)