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.


Roman letters a
Contact radius a 0 Equilibrium contact radius

Introduction
Microreaction technology has gained popularity and attention beyond the scope of academic labs due to the inherent benefits associated with the decrease in characteristic length scale, such as increased heat and mass transfer, safer handling of hazardous materials, precise control of residence time and reaction conditions (Tabeling 2006). Over the past two decades, this has led to a wide range of commercial equipment and applications (Jensen 2017). However, the use of this technology is still limited by their ability to handle solid material in small-scale devices, which elevates the risk of channel clogging and thus limited life span of these devices. The solid materials can either arise as insoluble starting reagents or products of a reaction, or solid particles can be deliberately added as heterogeneous catalyst for a chemical transformation (Ufer et al. 2011;Liedtke et al. 2015;Schoenitz et al. 2015;Pu and Su 2018). Therefore, the research effort in recent years was dedicated to the experimental investigation of clogging phenomena under various flow conditions in microfluidic devices, and clogging mitigation using external actuation such as ultrasound (Hartman et al. 2010;Kuhn et al. 2011;Hartman 2012;Flowers and Hartman 2012;Wu and Kuhn 2014;Fernandez Rivas and Kuhn 2016). The transport of particles in confined geometries is governed by complex phenomena, and the relevant interactions are sketched in Fig. 1. A microparticle suspended in a flow experiences numerous forces, i.e., inertia, attraction, repulsion, buoyant weight, drag and surface forces (Herzig et al. 1970;McDowell-Boyer et al. 1986;Sharma and Yortsos 1987). Due to the presence of attractive interactions, the microparticles arrange themselves in certain structures/patterns which are called aggregates. These aggregates act as precursor for the formation of even larger aggregates by collision and adhesion, this phenomena is called agglomeration. Depending on the balance between particle-fluid interactions and the inter-particle interactions fragmentation can also take place, i.e., large aggregate structures break down into small aggregates or even into single particles. Both fragmentation and agglomeration often occur simultaneously, but an imbalance between these mechanisms may lead to channel clogging (Henry et al. 2013).
A generic description of the clogging process is given by the interplay between channel constriction and channel bridging. The microparticles in the flow will continuously deposit on the channel walls (particle-surface interaction), and will therefore constrict the free cross-section available for fluid flow. In addition, the particle-particle interaction can lead to the formation of an arch of particles, which then bridges across the width of the constricted channel. And even without wall deposition clogging can occur by aggregation only, as a steady increase in the agglomerate size will eventually block the channel cross-section (Goldsztein and Santamarina 2004;Sharp and Adrian 2005;Mustin and Stoeber 2010).
Early stage aggregates are primarily formed via adhesive particle-particle collisions, and these initial particle aggregates then act as a nucleus for further aggregation, and in subsequent steps several aggregates attach together to generate a single large-sized cluster (Gudipaty et al. 2011). It is observed that the region near the channel walls and in the center of the channel exhibits a higher aggregation rate, along with fragmentation from the walls (Flamm et al. 2011). The deposition of particles on the channel wall can occur as either mono or multilayers, which is influenced by several factors, such as fluid characteristics (e.g., the ionic strength), hydrodynamic conditions (e.g., Reynolds number), substrate (e.g., Zeta potential) and particle properties. In case of high flow rates, multilayer deposition is only observed when the inter-particle repulsion magnitude is low Fogler 1998, 1999;Henry et al. 2012).
Numerical computations are a promising approach to further study the effect of particle-particle interaction and wall deposition in microfluidics. The Discrete Element Method (DEM) was applied to numerically investigate aerosol particle aggregate formation and collision-induced adhesion in a periodic straight channel under laminar flow conditions (Marshall 2007(Marshall , 2009. DEM computes the behavior of each particle corresponding to Newtons second law of motion, and therefore it is straightforward to include additional physics in the model, e.g., a cohesive force model and a contact model (Washino et al. 2013;Gao et al. 2014;Jasevičius et al. 2015). The results obtained by Marshall (2007Marshall ( , 2009 show that the adhesive properties of the particles and the flow rate have an enormous impact on the dynamics of the entire system. At large values of both flow rate and adhesive force, relatively large size aggregates are formed that deposit on the walls. This in turn affects the average velocity profile, which is disturbed due to these large aggregates adhering to the walls (Marshall 2007(Marshall , 2009Shahzad et al. 2016). A force-coupling numerical method was recently implemented to investigate microparticle aggregation in suspension flow in a three-dimensional microchannel (Agbangla et al. 2014). This method automatically considers inter-particle hydrodynamic interactions. However, due to the large computational cost the overall number of particles that can be simulated with this method is significantly lower as compared to DEM. A DLVO force has been applied to accurately model the adhesive, attractive and repulsive inter-particle interactions. Agbangla et al. concluded that inter-particle repulsive forces are the key parameter to initiate aggregate formation in the near wall region, which then will gradually block the channel.
In the current study, we applied a Discrete Element Method (DEM) and CFD-DEM coupling method to numerically simulate aggregation, agglomeration and fragmentation of rigid microparticles suspended in a Newtonian fluid flowing through a microchannel. We employed the JKR (Johnson-Kendall-Roberts) contact model to simulate particle-particle collision and adhesion in the system (Johnson et al. 1971). This contact model has recently been Fig. 1 Interactions governing the behavior of solid particles in microchannels (Wu and Kuhn 2014). a Deposition of particles is initiated by particle-fluid interactions transporting the solid to the microchannel wall where it finally sticks due to a dominating particle-surface interaction. b Increasing the particle-fluid interaction by, e.g., increasing the fluid velocity will lead to resuspension. c The particles will agglomerate in the bulk of the fluid by particle-particle interac-tions; however, agglomerate break-up can again occur when the particle-fluid interactions overcome the inter-particle interactions. d The clogging phenomena itself is governed by all three interactions, and usually occurs via bridging of a constricted microchannel cross-section. Reproduced with permission from Wu and Kuhn (2014). Copyright 2014 Tekno Scienze Srl successfully applied to study aggregation, agglomeration and fragmentation processes in particulate systems (Marshall 2007(Marshall , 2009Shahzad et al. 2016), and to identify the cluster generation in shear-induced coagulation of microparticles flowing in water (Kroupa et al. 2014(Kroupa et al. , 2015. The present work builds on the previous publication by Shahzad et al. (2016) which was the first to implement the JKR contact model to investigate the dynamics of aggregation, agglomeration and fragmentation for micro-gel particles. To achieve a better understanding of these interaction phenomena, particularly in the entrance section of a microchannel, where the initial particle and wall contact occurs, we study the flow of solids and fluid in a non-periodic domain. In this setup, the entire microchannel is initially filled with the fluid, and based on a pre-defined solid volume fraction, particles are injected with the incoming fluid. We systematically address the aggregation and deposition behavior of the particles as a function of the applied Reynolds number and strength of the adhesive force. Furthermore, we compare simulations in two-and three-dimensional domains and we also investigate the influence of one-way (DEM) vs. two-way coupling (CFD-DEM) in microchannels.

Mathematical model
In this work, we investigate a suspension of rigid microparticles flowing through a microchannel, which is represented as a two-(2D) and three-dimensional (3D) cylindrical channel of length L. Figure 2 depicts a schematic overview of both computational domains.
For the 2D channel geometry (Fig. 2a), the total length of the channel is eight times the channel width H, hence a fully developed velocity profile is achieved in the gap. Circular rigid microparticles with a specified diameter d P are suspended in the fluid. A Cartesian coordinate system with the origin in the center of the microchannel is placed in the inlet section, where x represents the fluid flow direction. The 3D cylindrical geometry (Fig. 2b) has a total length of eight times the channel diameter 2R cyl , which also leads to a fully developed velocity profile. Here, we consider a 3D cylindrical coordinate system with the origin at the center of the channel inlet section, where x also represents the fluid flow direction. In the 3D system, spherical microparticles with a specified diameter d P are suspended in the fluid.
In the one-way coupling approach, the flow field is assumed to be unaffected by the particle interactions. However, in the two-way coupling approach, both flow field and particle dynamics are computed simultaneously. Initially, a steady flow field is achieved in the microchannel, upon which particles are randomly injected on the inlet patch by defining the 'Particle Adding Frequency' f add (see Table 1) in such a way that the overall volume fraction remains constant throughout the entire channel. In addition, the velocity of the injected particles is set equal to the local fluid velocity. We impose non-periodic boundary conditions, as soon as a particle crosses the outlet patch it is instantaneously removed from the computation. The developed DEM code used here has been validated in a previous work (Shahzad et al. 2016) by comparing reported results of a particulate aerosol flow (Marshall 2007).

Discrete element method (DEM): one-way coupling approach
In the implemented one-way coupling approach, the microparticles are transported by the fluid while experiencing collisions with each other, however, the flow field of the fluid carrying the suspended microparticles is not modifiable (Marshall 2007). The overall dynamics of the microparticles is simulated by employing the Discrete Element Method (DEM), i.e., we employ a Lagrangian model to calculate the trajectories of the particles by solving the equations of motion for translational and angular velocities, respectively (Cundall and Strack 1979;Landry et al. 2003;You et al. 2004), where, m and I represent the particle mass and momentum of inertia, respectively, v and are the particle translation and angular velocities, respectively. We assume that no additional forces (or torques) arise due to sliding and twisting motions of particles interacting with each other. These assumptions can be justified as the Reynolds and Stokes numbers imposed in the simulations are rather small (Shahzad et al. 2016). Consequently, the only forces experienced by the particles are the fluid drag force F D , and F A the sum of the van der Waals adhesive and the elastic collision forces. The fluid torque M D and the sum of the van der Waals adhesion and elastic collisions torques M A are not considered, Eqs. (1) and (2) are decoupled as the particle translational motion is unaffected by its angular rotation (Shahzad et al. 2016). We also assume a 'freeze' boundary condition, i.e., when a particle gets in contact with the channel wall it will be permanently adhered. In the one-way coupling approach, the fluid-particle interaction close to the wall is minimal, and therefore resuspension of particles is not expected. However, for rather low adhesion forces, very few particles are expected to attach to the walls, but nevertheless, the flowing aggregate dynamics remain unaffected. The force terms used in Eq.
(1) are discussed in Sect. S1 of the Supporting Information, and for further information on the model implementation we refer to Shahzad et al. (2016).

CFD-DEM: two-way coupling approach
The CFD-DEM coupling approach is adopted to study aggregate formation and clogging in the 3D cylindrical channel while incorporating particle-particle, particle-wall, and particle-fluid interactions. In this coupled approach, the continuous liquid phase is represented by solving the volume-averaged Navier-Stokes equations using the opensource toolbox OpenFOAM, while employing DEM to model the particulate phase using the open-source toolbox LIGGGHTS. The open source CFDEM-coupling (Kloss et al. 2012) is used for data exchange between the two solvers at regular intervals; generally several DEM time steps are coupled with a single CFD time step. In the following Sects. 2.2.1 and 2.2.2, the equations and parameters used in the CFD-DEM model are briefly discussed.

DEM for solid phase
As outlined above, DEM is based on solving Newton's second law to compute translational and rotational velocity and position of each particle in time. All the equations in this section are computed using LIGGGHTS implementing Zhou et al. (2010) notations: where m and I are the particle mass and momentum of inertia, respectively, v and are the translational and angular velocities of the particle, respectively, F n and F t are the normal and tangential contact force between particles i and j, respectively, F pf is the particle-fluid interaction force, F g is the gravitational force, and M t and M r are the tangential and rolling friction momentum acting on i and j particles, respectively. In this work, we neglect gravitational and non-contact forces Fluid velocity (Re = 0.75, 1.0, 1.5) U f = 2.91, 3.8889, 5.8333 mm/s Time step size (Re = 0.75, 1.0, 1.5) Δt = 1.16e −6 , 1.56e −6 , 2.36e −6 s (i.e., electrostatic interaction or van der Waals forces), as the hydrodynamic and contact forces are orders of magnitude larger. We used the Johnson-Kendall-Roberts (JKR) contact model (Johnson et al. 1971) as described in Sect. S1.2, which was also validated for hard gel microparticles in our previous work (Shahzad et al. 2016). Table S1 in the Supporting Information summarizes the equations to calculate the forces and parameters in the DEM part of the CFD-DEM coupling.

CFD for liquid phase
Once all the force and torque terms are numerically solved in DEM, i.e., the trajectories and velocities for each particle are known, the CFD computation is performed to solve the motion of an incompressible fluid in the presence of the solid phase. For this, a modified set of the Naiver-Stokes equations based on the local volume average method are used (Gidaspow 1994;Zhou et al. 2010). In the present work, we use 'Model type A' that computes a shared pressure drop between the solid and fluid phase, in contrast to 'Model type B' which considers only the pressure drop for the fluid phase (Blais et al. 2016).
where u f and p are the fluid velocity and pressure terms, respectively, f , and ΔV are the fluid volume fraction in each cell, fluid viscous stress tensor, and the volume of the computational cell, respectively. K pf is a scalar term used to scale the magnitude of the momentum exchange force, F A pf is the volumetric particle-fluid interaction force and F pf,i is the particle-fluid interaction force, f d,i represents drag, f ∇p,i the pressure gradient and f ∇ ,i the viscous stress (or shear stress) on each particle, respectively. Furthermore, we used an unresolved CFD-DEM coupling approach, which is only feasible when the particle size is smaller than the cells of the computational grid.
The CFD-DEM coupling strategy is detailed in Sect. S3 in the Supporting Information. For the discretization of the 3D case we generated a structured mesh consisting of 84,000 cells. The equations are solved using the cfdemSolverPiso solver, which is based on the PISO algorithm, where two PISO loops are performed for each time step and within each PISO loop 2 non-orthogonal correction loops are applied.

Results and discussion
From the geometric details, fluid and particle properties (see Table 1 In the present study, we fixed the geometry and particle size, and systematically varied the fluid velocity U f and the surface energy of the particles to characterize their effect on particle-fluid, particle-particle, and particle-wall interactions. The resulting Reynolds numbers, elasticity and adhesion parameters are listed in Table 2. In the following, the complexity of the simulations is increased incrementally by characterizing particle aggregation and clogging in a 2D microchannel (one-way coupling), 3D microchannel (one-way coupling), and finally a 3D microchannel using two-way coupling.   Table 2. All simulations are performed sufficiently long to achieve steady state, which is confirmed by monitoring the aggregation statistics. The only exception is case C.4 (Fig. 5d), as the channel clogged without ever reaching steady state.

Reynolds numbers
It is generally observed that an increase in the adhesion parameter generates larger aggregates at a fixed Reynolds number. Considering the lowest Reynolds number Re = 0.75 and adhesion parameter = 15 (Fig. 3a), only small size aggregates are present which are mostly formed of 2-3 particles, with the remaining suspended particles leaving the flow domain without interacting with each other or the walls. Gradual increase in the adhesion parameter results in an increased formation of aggregates, especially in the near wall region. This observation is in line with the process of orthokinetic aggregation, in which the aggregation rate scales with the local shear rate, which is highest in the near region for the here considered developed laminar flow (Elimelech et al. 2013). These initially formed aggregates increase in size as they pass through the channel and some wall deposition is observed (Fig. 3c, d). The first monolayer of particles will permanently adhere to the channel walls according to our implemented freeze boundary condition. However, the second layer can attach and detach from this monolayer depending on the balance between adhesive and fluid-induced friction forces.  When increasing the Reynolds number to Re = 1.0 , aggregate formation in the near wall region is already observed at the lowest adhesion parameter (Fig. 4a). Increasing the adhesion parameter at this particular Reynolds number promotes increased aggregate formation and wall deposition. However, the detachment of aggregates induced by particle-fluid interaction is also observed. At the highest value of the adhesion parameter = 15,000 , the formation of dendritic structures attached to the channel wall is observed. While protruding into the flow, these structures stay attached due to the large adhesion force. Free flowing aggregates can collide with these structures, and upon impact either attach and merge with them, or result in breakage and detachment of the dendritic structure.
Increasing the Reynolds number further to the highest considered value of Re = 1.5 , an increase in the deposition and the aggregation dynamics is observed (Fig. 5). A monolayer of particles deposited on the walls is observed for adhesion parameters exceeding = 150 , and furthermore, the particle aggregates formed in the near wall region are larger in size compared to the lower Reynolds number cases. This indicates an increased probability of particle-particle collision when increasing the fluid flow rate (Elimelech et al. 2013). As observed earlier, when increasing the adhesion parameter value the number of aggregates and their sizes increase as well, which in the case of the largest adhesion parameter of = 15,000 leads to increased merging of free flowing and wall attached dendritic structures (Fig. 5d). These structures bridge the entire channel height and drastically reduce the effective cross-sectional area of the microchannel, thus indicating the onset of channel clogging.
To study this particle aggregation dynamics in more detail, the entire flow channel is subdivided into four equal volumes, which are highlighted by different colors in Whenever the center of mass of an aggregate crosses these intersections, its size, defined as the number of particles forming it, is recorded. This information is then used to obtain the spatial distribution of the percentage of aggregates (denoted as %Aggregate ) of a certain size according to where C i,s is the number of aggregates of size s (with s ≥ 1 ) crossing the intersection i, and the denominator represents the total number of aggregates crossing intersection i. Figure 6 depicts the percentage of aggregates %Aggregate consisting of N c particles per aggregate at the three intersections for all considered 2D microchannel cases. Each figure panel contains three color bars, which represent the three intersections at x∕H = 2.0 (green bar), x∕H = 4.0 (blue bar), and x∕H = 6.0 (red bar). The individual figure panel columns represent fixed Reynolds numbers (increasing from left to right), while the rows correspond to fixed values of the adhesion parameter (increasing from top to bottom). These statistics highlight that for all simulated cases the number of non-aggregated single particles is decreasing in the streamwise direction, due to the aggregation processes visually observed in Figs. 3, 4, and 5. For the smallest value of the adhesion parameter ( = 15 ) and Reynolds numbers of Re = 0.75 and Re = 1.0 the aggregation process is still relatively slow, with the majority of aggregates being comprised of 1-6 particles per aggregate. For all other cases, a substantial fraction of the formed clusters consist of seven or more particles ( N c ≥ 7 ), and the percentages of aggregates in this size range increase with increasing Reynolds number and adhesion parameter. It is also worth noting that the formation of these large clusters appears to be primarily taking place in the second half of the microchannel, which can be explained by the merging of smaller free flowing aggregates.
The accumulation of particles inside the microchannel is calculated by comparing the temporal evolution of the total number of particles present in the channel N tot (t) with the steady-state value of the total number of particles present in the channel in the absence of adhesive forces N tot (t → ∞, = 0) according to The temporal evolution of this dimensionless number of accumulated particles is depicted in Fig. 7 for all the simulated cases. Up to a dimensionless simulation time of t 2D = 10 , all curves overlap and approach the steady-state value of N tot = 1 . This transient behavior is caused by the simulation procedure, the microchannel is initially filled with the fluid phase only, and the particles are introduced via the inlet patch according to a specified particle adding frequency.
Hardly any particle accumulation is observed for all flow rates for the lowest values of the adhesion parameter = 15 and = 150 and the two lowest values of the Reynolds number ( Re = 0.75 ; Re = 1.0 ) as N tot ≈ 1 over time for these cases. As depicted in Fig. 6 aggregate formation is also occurring for these five cases, but all formed aggregates are free flowing and exit the microchannel without deposition on the channel walls. Upon a further increase of the flow rate ( Re > 1.0 ) and the adhesion parameter ( > 150 ) particle accumulation is observed, which is increasing over time and with the magnitude of the rate of increase depending on the specific values of Re and . As expected, the largest particle accumulation is found for the largest considered adhesion parameter of = 15,000 . For a Reynolds number of Re = 1.5 , the largest rate of increase in the number of particles accumulated in the microchannel is observed, which indicates the onset of channel clogging as shown in Fig. 5d. The oscillations observed in all the curves are associated with the continuous generation of aggregates and their exit from the channel domain. The amplitude value from those oscillations provide us generic information about aggregate sizes. The averaged particle and streamwise fluid velocity for the 2D microchannel cases are shown in Fig. 8. To obtain the averaged velocities, the microchannel is subdivided into regions with a dimensionless width of y∕H = 0.025 extending along the entire channel length L in flow direction (these regions are indicated in Fig. 8 by vertical lines). The shown particle and fluid velocity profile corresponds to the average velocity of all particles, respectively fluid elements, in each region.
As one-way coupling is used for the 2D microchannel simulations, the streamwise fluid velocity profile equals the profile for fully developed laminar flow. No significant deviation between the fluid and particle velocity is observed in the central region of the microchannel ( −0.25 < y∕H < 0.25 ) for adhesion parameter values ≤ 1500 and all considered flow rates. In the near wall regions ( y∕H < −0.4 and y∕H > 0.4 ), the particle velocity is smaller compared to the fluid velocity, which is caused by aggregate formation and the formation of dendritic structures attached to the microchannel walls. For the largest considered value of the adhesion parameter = 15,000 , a significant decrease of the particle velocity relative to the fluid velocity along the entire microchannel width H is found. This decrease in particle velocity is enhanced with increasing Reynolds number, which is caused by increased wall attachment and merging of free flowing with wall attached aggregates. For Re = 1.5 , the particle velocity approximates zero in the entire channel width which finally leads to clogging. These changes in the particle velocity profiles underline the observations in Figs. 3, 4, and 5. The initial aggregate formation occurs in the near wall region of the microchannel, and developing further into wall deposition and subsequent channel bridging by the growth and merging of aggregates.

One-way and two-way coupling approaches in a cylindrical 3D microchannel
In a next step, the particle aggregation is investigated in a cylindrical 3D microchannel. Figure 9 depicts the particle distribution for a Reynolds number of Re = 0.75 and the four different values of the adhesion parameter = 15, 150, 1500, and 15,000, these simulations are performed using the one-way coupling approach. When comparing with the corresponding 2D microchannel results (Fig. 3), it is observed that the aggregates formed in the 3D case consist of more particles. Furthermore, the monolayer of particles attached to the wall is formed at identical adhesion parameter values ( = 1500 and = 15,000 ), but the deposition process starts earlier in the 3D channel. The averaged particle and streamwise fluid velocity profiles for the 2D and 3D microchannel are compared for Re = 0.75 and = 15,000 in Fig. 10. While the overall averaged particle velocities are similar between the 2D and the 3D case, the relative deviation to the streamwise fluid velocity profile is more pronounced for the 3D case. This behavior is caused by enhanced One-way coupling approach: particle distribution in the 3D channel for Reynolds number Re = 0.75 and adhesion parameters a = 15 , b = 150 , c = 1500 , and d = 15,000 . The particle distribution is plotted at the end of the simulations ( t 3D = 50) aggregation dynamics, and therefore, larger aggregate sizes in the 3D microchannel, which can be quantified by comparing the percentages of aggregates %Aggregate consisting of N c particles per aggregate at the three intersections x∕H = x∕2R cyl = 2.0, 4.0, and 6.0, depicted in Fig. 11. A close resemblance in particle aggregation dynamics is found, with the exception that the aggregates formed in the 3D channel are larger ( N c ≥ 19 ), due to the fact that the particles have an increased cross-section and volume available to arrange themselves into 3D aggregate structures. Correspondingly, the influence of both the Reynolds number and the adhesion parameter values on the percentage of aggregates is more prominent in the 3D case. At the largest considered adhesion parameter = 15,000 , more than 70% large size aggregates are formed at the streamwise position x∕H = x∕2R cyl = 6.0 as shown in Fig. 11g, h. Furthermore, in the 3D cases, very few intermediate size aggregates ( 10 < N c < 18 ) are observed, which additionally highlight the enhanced aggregation dynamics in 3D.
Advancing to the two-way coupling approach, we are simulating the particle-fluid, particle-particle, and particle-wall interactions for Reynolds numbers Re = 0.75 and Re = 1.5 for an adhesion parameter of = 15,000 . The corresponding particle distributions in the 3D channel are depicted in Fig. 12.
Comparing the one-and two-way coupling results for a Reynolds number of Re = 0.75 (Figs. 9d, 12a), it is found that aggregate formation is delayed in the two-way coupling approach. Hence, larger aggregates are formed only in the last section of the microchannel ( x∕2R cyl > 6.0 ), and only limited wall-deposited particles are observed in the channel. This difference can be explained by the fact that the twoway coupling approach is not assuming a freeze boundary condition at the wall, but allows the resuspension of deposited particles if the particle-fluid interaction is able to overcome wall adhesion. This delayed aggregate formation when advancing from the one-way to two-way coupling approach is also visualized in Fig. 13, which depicts the comparison of the percentages of aggregates of a certain size at the three intersections x∕2R cyl = 2.0, 4.0, and 6.0 obtained for the oneand two-way coupling approaches.
Comparing the one-way and two-way coupling approaches, an approximately 10% decrease in large aggregates ( N c ≥ 19 ) in the first section of the microchannel is found for the CFD-DEM case. Small aggregates are observed close to the entry section for the one-way coupled case, whereas in the two-way coupled case, large aggregates are formed while flowing through the microchannel. A further increase in Reynolds numbers to Re = 1.5 also exhibits the same phenomena of delayed aggregate formation. Thus, the observed trends show that high values of the Reynolds number and the adhesion parameter promote wall deposition and aggregate formation, which can finally lead to microchannel clogging.
Furthermore, the two-way coupling approach also allows to study the effect of the particle aggregates on the fluid velocity field. Figure 14 depicts a contour of the normalized instantaneous streamwise velocity V∕U f at the dimensionless time t 3D = 10 for a Reynolds number of Re = 1.5 ( = 15,000).
The presence of the large agglomerates significantly affects the velocity distribution in the channel, reducing the fluid velocity in the center of the microchannel resulting in increased velocities in the near wall region of the channel. This dynamic change of the streamwise velocity profile is shown in more detail in Fig. 15.
This figure depicts the profiles of the instantaneous streamwise fluid velocity extracted at x∕(2R cyl ) = 6.0 at dimensionless times t 3D = 0 , 5, 10, and 15, respectively. At t 3D = 0 , the undisturbed parabolic laminar profile is observed, which is increasingly affected in the near wall region due to the aggregate generation, which results in a higher velocity close to the wall for t 3D = 10 and t 3D = 15 . This increased near wall velocity will also lead to an increased wall shear stress, which will promote the detachment of deposited particles. These changes in the velocity field are also captured in the pressure drop over the computational domain. Figure 16 depicts the normalized pressure drop Δp∕Δp 0 , where Δp 0 represents the pressure drop of the developed flow without particles, for both considered Reynolds numbers. It is observed that the pressure drop will initially ( t 3D < 10 ) increase non-linearly with time due to the presence of particles, and subsequent agglomeration and wall deposition. For longer simulation times, the pressure drop starts to fluctuate, which is explained by the formation of larger agglomerates in the bulk of the fluid flowing through the channel, leading to large disturbances of the velocity field. The finding of fluctuating pressure drops over the entire domain is in line with experimental observations of solid forming reactions (Hartman et al. 2010;Kuhn et al. 2011).

Conclusions
In this study, we investigated the behavior of suspensions containing rigid microparticles flowing through rectangular (2D) and cylindrical (3D) microchannels using a Discrete Element Method (DEM) and a CFD-DEM coupling method. Simulations were performed at different values of the Reynolds number and the adhesion parameter, with an overall volume fraction of particles of 0.1 (2D) and 0.05 (3D). The particle-particle, particle-wall, and particle-fluid interactions were characterized in terms of statistics of aggregate formation and wall deposition.
The results reveal that particle aggregation is initiated in the near wall region of the microchannel. Furthermore, wall attachments are already observed at considerably low values of the Reynolds numbers and the adhesion parameter. As expected, increasing the value of the adhesion parameter results in larger agglomerate sizes and in their earlier formation. A further increase in the Reynolds numbers results in increased particle-aggregate, aggregate-aggregate and aggregate-wall interactions. For the condition of largest considered Reynolds numbers and adhesion parameters, large-sized aggregates are formed together with enhanced wall deposition, which finally leads to microchannel blockage. Comparing the 2D and 3D one-way coupled simulation cases, larger agglomerate sizes are observed for the 3D case at the same conditions.
Advancing to the two-way coupled CFD-DEM simulations, the influence of the particles on the fluid velocity field is characterized. Compared with the one-way coupled 3D case, a delay in wall deposition is observed, however, the overall size of the formed agglomerates is similar. The presence of large size agglomerates bridging the channel leads to the fluid bypassing them, which in turn results in increased near wall velocities and shear stress, which can resuspend attached particles from the wall. And also in the CFD-DEM simulations, we observe increased aggregation dynamics with increasing Reynolds number.
In conclusion, as the three-dimensional two-way coupling approach is able to resolve the dynamics in wall shear stress and associated delayed particle deposition, it represents the more realistic approach to capture particle-fluid dynamics in microfluidic applications and to predict clogging. As such, these CFD-DEM simulations provide the opportunity to guide experimental studies, as they allow to evaluate the clogging risk of a particular microfluidic device, and its ability to handle different particles with different adhesion parameters. Furthermore, the two-way coupling approach also enables the prediction of the optimum fluid dynamic conditions for the production of particles and agglomerates, with applications in, e.g., nanomaterial synthesis.