Low Reynolds Number Interactions between Colloidal Particles near the Entrance to a Cylindrical Pore

Journal of Colloid and Interface Science 229, 311–322 (2000) doi:10.1006/jcis.2000.6986, available online at http://www.idealibrary.com on

Low Reynolds Number Interactions between Colloidal Particles near the Entrance to a Cylindrical Pore Venkatachalam Ramachandran,∗ Ramachandran Venkatesan,∗ Gr´etar Tryggvason,† and H. Scott Fogler∗,1 ∗ Department of Chemical Engineering and †Department of Mechanical Engineering and Applied Mechanics, University of Michigan, Ann Arbor, Michigan 48109 Received April 2, 1999; accepted May 30, 2000

The interaction between stable colloidal particles arriving at a pore entrance was studied using a numerical method for the case where the particle size is smaller than but of the same order as the pore size. The numerical method was adapted from a front-tracking technique developed for studying incompressible, multifluid flow by S. O. Unverdi and G. Tryggvason (J. Comp. Phys. 100, 25, 1992). The method is based on the finite difference solution of Navier– Stokes equation on a stationary, structured, Cartesian grid and the explicit representation of the particle–liquid interface using an unstructured grid that moves through the stationary grid. The simulations are in two dimensions, considering both deformable and nondeformable particles, and include interparticle colloidal interactions. The interparticle and particle–pore hydrodynamic interactions, which are very difficult to determine using existing analytical and semi-numerical, semi-analytical techniques in microhydrodynamics, are naturally accounted for in our numerical method and need not be explicity determined. Two- and three-particle motion toward a pore has been considered in our simulations. The simulations demonstrate how the competition between hydrodynamic forces and colloidal forces acting on particles dictate their flow behavior near the pore entrance. The predicted dependence of the particle flow behavior on the flow velocity and the ratio of pore size to particle size are qualitatively consistent with the experimental observations of V. Ramachandran and H. S. Fogler (J. Fluid Mech. 385, 129, 1999). °C 2000 Academic Press

INTRODUCTION

Understanding the motion and interaction of colloidal particles at the entrance of pores is crucial for the control and design of diverse processes such as solid–liquid separation by filtration, osmotic phenomena, transport of groundwater colloids, and the permeability reduction due to migrating fines in enhanced oil recovery. The case of a cylindrical pore is relevant because Nuclepore filters having track-etched cylindrical pores are extensively used for aerosol filtration studies (9) and cell separations. In addition, because of their well-defined pore shapes and pore size distributions, Nuclepore filters serve as a model porous medium in experimental studies directed toward understanding retention 1

To whom correspondence should be addressed.

mechanisms of stable colloidal particles within porous media. Two of the authors (V.R. and H.S.F.) have studied the plugging of pores in Nuclepore membranes during flow of dilute, aqueous, stable colloidal suspensions by the mechanisms of hydrodynamic bridging (15) and multilayer deposition (14). The former study, pore blockage by hydrodynamic bridging, is the main motivation for this work and will hereafter be referred to as R-F (1999). Hydrodynamic bridging is the phenomenon of blocking of pores by simultaneously arriving stable particles whose sizes are smaller than the pore size. The term “stable” refers to the particles in the suspension remaining dispersed and not aggregating during Brownian collisions. At a sufficiently high flow velocity yet still in the low Reynolds number regime, hydrodynamic forces at the pore entrance can overcome interparticle and particle–pore surface electrostatic repulsion resulting in particle association and blockage of the pore constriction (see Fig. 1). A set of experiments from R-F (1999) is discussed next to illustrate particle retention by hydrodynamic bridging. Figure 2 shows the effect of velocity on particle retention at a fixed particle concentration during flow of 0.22-µm polystyrene carboxylate latex particles, at a constant volumetric flowrate, through Nuclepore membranes having 1-µm pores (the ratio of actual pore size to actual particle size is 3.7). In this figure, the pressure drop across the membrane during suspension flow, normalized using pressure drop during the flow of deionized water, is plotted as a function of the total volume of suspension injected. In Fig. 2, the increase in the normalized pressure drop during suspension flow indicates that particles are being retained within pores. Particle retention is not expected in these experiments because the small particle size relative to the pore constriction size precludes straining (physical capture of particles larger than the pore constriction) while particle–pore surface electrostatic repulsion (like surface charges) prevents deposition. One also observes in Fig. 2 that as the velocity is increased both the rate and the extent of plugging of the membrane increases for the same number of particles injected. In other words, the efficiency of particle retention by hydrodynamic bridging increases with increasing flow velocity. As discussed in R-F (1999), the observed particle retention in these experiments as well as the observed effect of velocity can be explained on the basis of hydrodynamic bridging. The effect of velocity on retention by

311

0021-9797/00 $35.00

C 2000 by Academic Press Copyright ° All rights of reproduction in any form reserved.

312

RAMACHANDRAN ET AL.

FIG. 1. Particle capture by hydrodynamic bridging.

bridging is unique because particle retention by straining is independent of flowrate (for a fixed particle concentration, curves for the different velocities would coincide if straining was the primary mechanism of particle retention). On the other hand, it was shown by Ramachandran and Fogler (14) that the effect of velocity on deposition is opposite to that observed for retention by bridging. R-F (1999) experimentally also showed that plugging by hydrodynamic bridging will not occur below a certain critical velocity. A critical velocity exists because hydrodynamic forces at the pore entrance must overcome a certain magnitude of the net colloidal repulsion between particles for bridging to occur. The magnitude of this critical velocity is determined by the surface charge properties of the particle and the porous medium, the ratio of the pore size to the particle size, and the flow geometry. They also calculated trajectories of two particles approaching a circular orifice and qualitatively showed why such a critical velocity exists. However, their calculations were approximate in nature because the interparticle hydrodynamic interaction was not included and the hydrodynamic forces acting on the particles were determined based on the hydrodynamic interaction of one particle with a circular orifice. The objective here is to extend that work by including both interparticle hydrodynamic interactions

FIG. 2. Effect of velocity on hydrodynamic bridging (0.22-µm PSC latex spheres).

near the pore entrance and particle pore–wall hydrodynamic interaction when two or more particles are present near the pore. Theoretical calculations of the interaction between particles at the entrance to a pore when the particle size is the same order of magnitude as but smaller than the pore size have not been attempted. This difficult situation is the focus of this paper. We will first briefly discuss relevant studies in the literature prior to describing our method. Relevant to the problem of interaction between particles at the pore entrance is the motion of a particle toward a pore. Because of the strong hydrodynamic interaction between the particle and the pore wall, the hydrodynamic resistance experienced by the particle differs substantially from the Stokes resistance of a particle in an unbounded fluid. The hydrodynamic correction factors for the force and torque are no longer scalars but vary with the position and the direction of the particle motion. The major difficulty in solving such problems is the discontinuous nature of the boundary representing the wall with the orifice (pore). Weinbaum et al. (23) have reviewed three important numerical techniques which have been used for solving complicated boundary value problems in bounded and unbounded Stokes flows. These three techniques are (a) the multipole collocation method, (b) the boundary integral technique, and (c) the multipole moment method. Weinbaum et al. have described the conceptual origin of each technique and discussed some of the boundary-value problems that have been solved using these techniques. Yan et al. (24) were the first to solve the three-dimensional problem of a finite sphere approaching a circular orifice. They used a combined multipole collocation and integral equation method in which the disturbance due to the sphere was represented by a multipole series while the orifice wall was represented using the single-layer potential. Their solution was however not general because of two reasons: (a) in their formulation, the sphere center is located on a plane containing the orifice centerline, and (b) the sphere cannot intersect the plane of the orifice due to limitations of their approach. Furthermore, owing to the excessive computational time required, the hydrodynamic resistance coefficients were only determined for the case of the particle radius equal to one-half the pore radius. Wu et al. (8) have presented an approach that overcomes these limitations. Their technique involves numerically solving a set of boundaryintegral equations for the motion of the sphere using an extension of the Green’s function for a point force moving either parallel or perpendicular to an infinite wall with a hole. Building on the work of Yan et al. (24) Kao et al. (9) developed a hydrodynamic interaction model for the collection efficiency of a Nuclepore filter, while Yan et al. (11) studied the structure of osmosis at the entrance and exit of membrane pores. The implementation of these techniques for the simultaneous motion of two or three particles toward the pore entrance will however be computationally very difficult. Another relevant study is by Vitthal and Sharma (27), who have modeled the deposition of spherical colloidal particles during flow through an assembly of spheres using Stokesian

313

LOW REYNOLDS NUMBER INTERACTIONS

dynamics (for a description of the Stokesian dynamics technique see (3)). In their study, the “pore” was the constriction in the space between the spheres in the assembly, and hence the boundary was not as sharp as in the case of a circular orifice. They calculated the trajectory of one particle at a time taking into account the hydrodynamic interactions and van der Waals attraction between the flowing particle and the sphere assembly including all previously deposited particles. In this manner, they demonstrated that previously deposited particles significantly affect the deposition process by both altering the flow field and acting as additional collectors. They did not consider the case of more than one particle simultaneously approaching the sphere assembly. In the area of computational fluid dynamics, the simulation of the flow of two or more immiscible fluids has received considerable attention. The computational technique used in this work for studying particle interaction at the pore entrance is based on the work of Unverdi and Tryggvason (20). We will refer to their work as U-T (1992) in this paper. U-T (1992) developed a fronttracking method for studying viscous, incompressible multifluid flows. The method is based on a finite difference approximation of the full Navier–Stokes equations and explicit tracking of the interface between the fluids (“front”). Using this method, they successfully studied the unsteady, two- and three-dimensional motion of drops and bubbles rising under a gravity field. A brief description of the technique is included in the next section. The primary advantage of U-T’s method, as can be seen from their simulation of motion of bubbles, is the rather natural way in which interfaces (in the case of more than one drop or bubble) interact during their simultaneous motion. The interaction between the particle and the surrounding fluid is naturally taken into account in the solution of the Navier–Stokes equations for the simultaneous motion of both phases. As a result, this method obviates the explicit calculation of the hydrodynamic forces. We have modified U-T’s method in two dimensions to the case of two and three particles approaching a circular pore. Our goal is to study the interaction between submicrometer, charged colloidal particles as they approach the pore. The interplay between the hydrodynamic forces acting on the particles and colloidal forces such as van der Waals attraction and double-layer repulsion is what dictates the particle flow behavior at the pore entrance. Because U-T’s method is inherently constructed to treat deformable particles, the method cannot be used in its original form for studying solid particles. However, by assigning high viscosities to the particles relative to the surrounding medium, “nearly” solid particles have been simulated. We have also incorporated a particle reconstruction technique to maintain the circular shape which gives physically reasonable results. Finally, the colloidal forces between particles have also been included in the simulations. The main advantage of our method is its ability to simulate the motion of two or more particles near the pore entrance which is very difficult using the semi-analytical, semi-numerical methods mentioned earlier because of the strong interparticle and particle–pore hydrodynamic interactions near the pore entrance. It should be pointed out that

simulations in two dimensions implies that the particles are not spheres but are circular cylinders. As will be shown in this paper, two-dimensional simulations are sufficient to qualitatively demonstrate the hydrodynamic bridging phenomenon reported in R-F (1999). For quantitative comparison with experiments simulations will need to be extended to three dimensions and is discussed later in this paper. The remainder of the paper is organized as follows: the numerical method is described in the next section. A detailed listing of the numerical scheme is included in the Appendix. We then present our simulation results where we first demonstrate the convergence of the numerical solution for the grid used and also verify that the particle–liquid interface has been sufficiently resolved. Next, we discuss the results from the simulations of two and three particles approaching the pore entrance. Finally, we conclude by summarizing the results and discussing the limitations and advantages of the present method. NUMERICAL METHOD

The method used here was developed by Unverdi and Tryggvason (20) to simulate unsteady, incompressible, multifluid flows. They have used this technique to study bubble motion and interaction between bubbles in two and three dimensions. We will briefly discuss the numerical method here; a detailed description of the technique can be found in U-T’s original paper. The schematic in Fig. 3a illustrates the situation

FIG. 3. (a) Computational domain for studying the interaction between particles near the entrance to a pore; (b) staggered, stationary, two-dimensional grid used to discretize the computational domain. u and v are the fluid velocity components in the x and y directions respectively.

314

RAMACHANDRAN ET AL.

being considered in this work. The motion and interaction of submicrometer particles toward a pore is being studied in two dimensions. By an appropriate choice of boundary conditions upstream of the pore, the pore represents a periodic unit in a Nuclepore-type membrane. To determine the motion of the particles, solution of the Navier–Stokes equations governing the unsteady motion of immiscible, viscous fluids is required. For Newtonian fluids, the Navier–Stokes equation in conservative form is given by ∂ (ρu) + ∇ · (ρuu) = −∇ p + ρg + ∇ · µ(∇u + ∇uT ) ∂t + σ κnδ(x − xf ) + Fcol ,

[1]

where u is the velocity field, ρ and µ are the discontinuous density and viscosity fields, respectively, xf denotes the position of the front, k is twice the mean interface curvature, σ is the surface tension, δ(x) is the Dirac-delta function, n is the normal to the interface, g is the acceleration due to gravity, and Fcol is the colloidal force per unit volume acting on the particles. Note that the surface tension force is nonzero only at the fluid interface and hence the representation as a delta function. Similarly, the colloidal force between particles will only be included for those region(s) of the domain occupied by the particles. It is assumed that the convective time scale is much smaller than the time scale for Brownian motion, and hence Brownian motion does not affect the particle flow behavior. This assumption will be acceptable for pore Peclet numbers greater than unity, which is the case in the R-F (1999) experiments. Because the fluids are incompressible, the continuity equation supplementing [1] reduces to ∇ · u = 0.

[2]

Equations [1] and [2] are solved for the rectangular twodimensional domain shown in Fig. 3a with a finite difference method. The spatial derivatives are approximated using secondorder finite differences on a staggered grid. For the advection in time, U-T used either an explicit Euler method or a second-order Adams–Bashforth method. However, the stability criterion that dictates the magnitude of the time step is too restrictive for the micrometer length scales and the large viscosity ratios used in our simulations. We have therefore implemented an Adams– Bashforth/Crank–Nicholson scheme which is implicit in time and second-order accurate (13, p. 164). In this method, the viscous terms which are dominant at low Reynolds numbers are evaluated implicity while the convective terms are evaluated explicity. The details of the implicit scheme are listed in the Appendix. The stationary, structured, staggered grid used in the calculations is shown in Fig. 3b. A no-slip boundary condition was specified for the velocities on the walls bounding the pore (boundaries 4–7). On the boundaries upstream of the pore entrance (boundaries 1–3), the velocity boundary conditions described in (10) were used. Manton determined the flow field ap-

proaching a Nuclepore filter in the low Reynolds number limit in the absence of particles. The boundary conditions used are also listed in the Appendix. The unique feature of U-T’s technique is the explicit representation of the fluid–fluid interface with an interface grid which is superimposed on the stationary Cartesian grid in Fig. 3b. In their study of bubble motion, this additional computational element ensured that the interface was sharp throughout the simulation. To avoid a jump in the magnitudes of ρ and µ across the interface, U-T (1992) “smeared” the interface such that its thickness was of the order of the mesh size. The fluid physical properties change smoothly over this thickness from the value inside the particle to the value outside. The interface is advected by moving each of the individual points in the interface grid. Because the fluids are incompressible, the velocity of a point in the interface grid equals the fluid velocity at that point which can be interpolated from the stationary grid. This technique, called the “immersed boundary technique,” was first developed by Peskin (12) to study the blood flow pattern in the heart. U-T’s treatment of the interface has also been incorporated in computational studies of cell deformation and adhesion (1). The sequence of computational steps involved is as follows: Given u, ρ, and µ at time t, the interface is advected to a new position at time t + 1t. The surface tension and colloidal forces for the new position(s) of the interface(s) are then calculated. Because the scheme is implicit, the elliptic equation for p must be solved simultaneously with the equation for u using ρ and µ at time t and t + 1t. The interfaces are advected to a new position from a knowledge of the new velocity field and the calculations are repeated till the interface has reached the desired final position. Because of the micrometer length scales involved, a de-dimensionalized form of the Navier–Stokes equations was used in the calculations. This scaling ensured that the magnitude of the distances and the velocity components used in the numerical scheme was O(1) for most of the Eulerian grid points. Faster convergence of the solution was obtained when the de-dimensionalized equations were used. As the interface approaches the pore, it will deform in response to the stresses exerted by the surrounding fluid. As portions of the interface stretch or contract during deformation, computational points are added or removed to maintain adequate resolution. While this monitoring of the interface grid was crucial for U-T’s simulations of bubble motion because of the extent of deformation in their work, the high viscosity and submicrometer size of the particles being simulated in this work resulted in such modifications not being frequently needed. Colloidal forces. When the simulations involved two or more particles, the colloidal forces between the particles were included in the simulations. For simulations involving more than two particles, the colloidal force acting on a particle was the sum of its colloidal interaction with all other particles. The magnitude of the colloidal force was calculated based on interactions between solid spherical particles. As mentioned previously, the two-dimensional

LOW REYNOLDS NUMBER INTERACTIONS

nature of the simulations implies that the particles are not spheres but are circular cylinders. Because the goal of this work is to qualitatively demonstrate hydrodynamic bridging reported in R-F (1999), and in anticipation of extending this work to simulations in three dimensions, sphere–sphere colloidal forces were used. The colloidal forces are assigned as body force contributions to all the Eulerian grid points that are within the particles. The colloidal interaction energy between two charged particles dispersed in a polar liquid is calculated in the framework of the DLVO theory and is a combination of van der Waals attraction (for particles of the same material) and electrostatic repulsion (like charge). A short-ranged repulsion attributed to solvation or structural forces (18) is also included here. The expression for the retarded London–van der Waals attraction between two spheres of equal size given by Schenkel and Kitchener (17) has been used here, VA (h) = −

Aa ¢, 12h 1 + 11.12h λ ¡

[3]

where VA is the attractive interaction energy, A is the Hamaker’s constant, h is the surface separation between the particles, λ is the characteristic wavelength of interaction (retardation length) often assumed to be about 100 nm, and a is the particle radius. The range of validity of this expression and its accuracy have been discussed by Gregory (19). The expression for the electrostatic repulsion between two identical colloidal particles used here is based on the linearized Poisson–Boltzmann equation and the Derjaguin approximation for constant potential interaction (16), VR (h) = 2πεaψ02 ln(1 + e−κh ),

[4]

where VR is the electrostatic interaction energy, ε is the permittivity of water, ψ0 is the surface potential of the particles, and κ is the Debye–Huckel parameter (κ −1 is the screening length). For κa ≈ 10 which is the case in our experiments, the interparticle repulsive force calculated using Eq. [4] is within 10% of the exact force for h ≤ 3κ when the potential is small (22). Feke et al. (6) have derived an expression for the short-ranged Born repulsion between two identical spherical particles µ ¶6 ½ 2 σ A 1 R − 14R + 54 60 − 2R 2 + VB (h) = 37800 a R (R − 2)7 R7 ¾ R 2 + 14R + 54 + , [5] (R − 2)7 where VB is the Born interaction energy, σ is the collision di˚ and R is the center-to-center particle separation ameter (4 A), scaled using the particle diameter (R = h/2a + 1). The total interparticle colloidal interaction potential, VT , is the sum of the above three interaction potentials. The magnitude of p–p the interparticle colloidal force, Fcol , is related to VT by the

315

FIG. 4. Comparison of the magnitudes of interparticle and particle–plane wall colloidal forces. In the vicinity of a pore, the particle–pore colloidal force will be less than that calculated based on the interaction between a particle and a solid plane wall.

equation ¯ p–p ¯ ¯F (h)¯ = − d VT (h) = − d [VA (h) + VR (h) + VB (h)], col dh dh [6] where the superscript p–p denotes interparticle interaction. Figure 4 shows a plot of the interparticle colloidal force as a function of the separation distance. The parameters specified in Fig. 4 are the same as those used in the trajectory calculations. It can be seen that the interparticle interaction is strongly repulsive up to very small separations, where the dispersion attraction starts to dominate. For particles to associate near the pore entrance and bridge the pore, the repulsive force barrier seen in Fig. 4 must be overcome. It should be noted that the interparticle potential energy for our system only exhibits a primary minimum because of the strong electrostatic repulsion. We have not considered systems where the interparticle potential is characterized by a primary minimum as well as a secondary minimum. For the case of a particle approaching a pore, exact calculation of the colloidal interaction between the particle and the pore is not straightforward because of the geometry of system (for example, see (2) for the van der Waals interaction and (5) for the electrostatic interaction between a particle and a charged pore in a charged planar surface). In the hydrodynamic bridging experiments of R-F (1999), the magnitude of the particle–pore wall colloidal force is small compared to the interparticle colloidal force because the pore surface charge is smaller than that of the particles. In Fig. 4, we have also plotted the particle–pore wall colloidal force to compare its magnitude with the interparticle colloidal force. The particle–pore wall colloidal force was calculated approximately based on the interaction between a spherical particle and a charged, solid plane wall. It can be seen that the maximum interparticle net repulsion is about 16 times greater than that between a particle and a solid plane wall. As mentioned previously, the hydrodynamic force acting on a particle is equal in magnitude and opposite in direction to the

316

RAMACHANDRAN ET AL.

total colloidal force acting on it, ¡ p–p p–w ¢ Fh = −Fcol = − Fcol + Fcol , where the superscript p–w denotes particle–pore wall colloidal interaction. If the magnitude of the particle–pore wall colloidal force is small compared to the interparticle colloidal force, then the particle velocity as it approaches a pore will be determined p–p primarily by the interparticle colloidal force; i.e., for |Fcol | À p–w |Fcol |, p–p

Fh ≈ −Fcol . We have therefore omitted the particle–pore wall colloidal interaction in our calculations. RESULTS AND DISCUSSION

Single-particle simulations. We first present simulations of a particle approaching the pore entrance. The purpose of the simulations was to validate the numerical technique used in this paper. In all simulations presented in this paper, a particle size of 0.5 µm has been used. The extent of particle deformation in response to flow field at the pore entrance is shown in Fig. 5 for three different cases; in the figures, the initial and final positions of the particles are shown. In these simulations, an aspect ratio (ratio of the pore size to the particle size) of 1.5, a 64 × 64 grid, a particle density to the solvent density ratio (ρp /ρs ) of 2, a Re of 0.05, and a particle resolution of 16 were specified. The Reynolds number is calculated based on the particle diameter and the streamwise velocity far from the pore; the particle resolution, R, is defined as the ratio of the particle diameter to the mesh size. Figure 5a shows the initial and final positions of the particle for the case of µp /µs = 10, where µp /µs is the ratio of the viscosity of the particle to the solvent viscosity. Because of the low relative viscosity of the particle, there is significant deformation of the particle at the pore entrance for the velocity specified. As we are interested in simulating the interactions of nondeformable, “solid” particles at the pore entrance, a high particle relative viscosity of 100 was assigned in our simulations.

FIG. 5. Dependence of the extent of particle deformation on relative viscosity: (a) µp /µs = 10; (b) µp /µs = 100; (c) µp /µs = 100 and particle reconstructed after every time step. In all the above single-particle simulations, Re = 0.05, σ = 1, and particle resolution (R) = 16. In the figures, the initial particle position and the final position when execution was terminated are shown.

It is known that the hydrodynamic drag coefficient for a fluid particle approaches that for a solid particle falling in an unbounded fluid when the ratio of the particle viscosity to the bulk fluid viscosity is greater than 50 (4). However, as can be seen from Fig. 5b, even µp /µs = 100 is insufficient to completely suppress deformation in the flow geometry being studied because of the high shear rates at the pore entrance. To completely eliminate particle deformation, we have implemented a particle reconstruction technique in some of the simulations. In this technique, the center of the particle is identified after each time step and the interface grid points are reassigned about the center based on the specified particle radius. Obviously, whether such a correction is acceptable depends on the extent of deformation in one time step. Simulations in which large particle viscosities and small time steps are specified will be suitable for using this technique in order to eliminate particle deformation. Simulations were performed to study the significance of including the surface tension force for σ in the range of 0–10 dynes/cm. It was seen from the simulations that the extent of deformation of the particle is almost identical for σ in that range. This result can be understood by considering the Eotvos number, Eo = 4ρs ga 2 /σ , characterizing the particle. Under the same flow conditions, the extent of particle deformation decreases as Eo decreases. For a particle radius of 0.5 µm, a solvent density of 1 g/cc, and a surface tension of 1 dyne/cm used in the simulations, the Eotvos number is of the order of 10−5 . Based on these results and the fact that high surface tensions for the submicrometer particle size specified affected the convergence of the solution, we have used σ = 1 dyne/cm in all our simulations. Note that the typical magnitudes of the surface tension for the water–organic polymer solid interface and the water–inorganic solid interface are 10 and 100 dynes/cm, respectively (7). The simulations presented in Fig. 5 use a 64 × 64 grid and are fully converged solutions. This was confirmed by comparing with simulations using grid spacing smaller than that for a 64 × 64 grid (up to 144 × 144 grids were checked). For the above reason, a particle resolution of 16 was used in most of the simulations reported in this paper. As the thickness of the particle–solvent interface is of the order of the mesh size, a resolution R = 16 was found to be sufficient for correct results in our 2- and 3-particle simulations. The need for sufficient resolution is demonstrated in Figs. 6a and 6b for two particles entering a pore. In these figures a comparison is made between two simulations which are identical except for the mesh size or resolution which is R = 8 in Fig. 6a and R = 16 in Fig. 6b. It can be seen that a resolution of R = 8 is insufficient because the smeared interfaces in Fig. 6a are interacting even though the particles are far apart. A resolution of R = 16 was found to be sufficient to obtain both physically reasonable results and practical simulation times. The size of the computation domain used in the simulations depended on the aspect ratio being considered and the initial position of the particles. That the sizes of the computational domain were sufficiently large to obtain correct results was verified by comparing with simulations having larger computational

LOW REYNOLDS NUMBER INTERACTIONS

317

FIG. 7. Initial positions of particles specified in the 2-particle simulations: R = 16. The initial separation distance corresponding to these positions is 50 nm.

FIG. 6. Resolution requirements for simulating particles interaction: (a) R = 8; (b) R = 16. The velocity vectors within the computational domain is shown in each case.

domains. Also, the implicit scheme and the solutions for the time steps used were verified by comparing with suitable simulations using U-T (1992)’s explicit Euler method. A simulation with an 80 × 80 grid-spacing with initial positions specified in Fig. 7 took about 7 h (on HP 9000/782 UNIX workstation with a 200-MHz PA-RISC processor). In the range of grid-spacings used, simulation time scaled as the number of grid points raised to 1.9, i.e., weaker than quadratic dependence. Two- and three-particle simulations. The main focus of this paper while simulating interaction between particles arriving at a pore entrance is to understand the effect of velocity on hydrodynamic bridging shown in R-F (1999). Hydrodynamic bridging refers to the plugging of a pore constriction by simultaneously arriving stable colloidal particles that are smaller in size than the pore constriction. Bridging is the result of competition between the hydrodynamic forces acting on the particles at the pore entrance and the electrostatic repulsion between the particles. If the hydrodynamic forces are strong enough to overcome the net colloidal repulsion near the pore entrance, bridg-

ing will occur. We demonstrate this by simulating two particles approaching a pore under low Reynolds number conditions. Figure 7 shows the initial positions of the particles upstream of the pore entrance. The particles are symmetrically placed about the pore axis and are equidistant from the plane of the pore. Colloidal forces between the particles have been included. The initial separation between the particles is about 50 nm. For particle separation much larger than 50 nm, colloidal forces did not influence the particle flow behavior because of the range of the colloidal forces for the conditions specified. A high zeta potential of 100 mV was used in the simulations to accentuate the influence of colloidal repulsion. Figure 8 shows the separation distance between the particles as a function of the position of their centers as they approach the pore entrance. The aspect ratio in this simulation is 1.5. By using the technique of particle reconstruction, nondeformable particles are being simulated. For a fluid

FIG. 8. Effect of velocity on particle flow behavior: α = 1.5, ψ p = 100 mV, nondeformable particles.

318

RAMACHANDRAN ET AL.

FIG. 9. Effect of aspect ratio on particle flow behavior: Re = 0.1, ψ p = 100 mV, nondeformable particles.

velocity corresponding to Re = 0.02, the particle separation decreases slowly as they approach the pore entrance. The particle separation decreases despite the strong electrostatic repulsion because of the converging nature of the flow into the pore. At a sufficiently high velocity corresponding to Re = 0.1, the particle separation decreases rapidly relative to the Re = 0.02 case and is essentially zero when the particle centers are a distance 1.5a from the pore. Therefore, for Re = 0.1, the hydrodynamic forces are strong enough to overcome the colloidal repulsion and cause particle flocculation. To show the influence of the aspect ratio on particle interaction, simulations were performed for a range of aspect ratios but for the same farfield velocity and particle zeta potential. Figure 9 shows the separation distance between the particles as a function of the particle location for the different aspect ratios. As the aspect ratio increases, at a given axial distance from the pore, the magnitude of the hydrodynamic force acting to push the particles together decreases. As a result, as seen in the simulations with aspect ratios 2.5 and 3, the hydrodynamic forces are insufficient to flocculate the particles as they enter the pore. Simulations showed that the particle zeta potential did not affect the flow behavior significantly. Zeta potentials in the range of −40 to −100 mV were studied. The simulations showed that as the zeta potential decreases, the particles flocculate farther away from the pore. This trend can be understood from the fact that for low Reynolds number particle motion toward a circular orifice, the hydrodynamic force acting on a particle tending to push it toward the pore axis is strongest close to the plane of the pore and decreases with increasing distance away from the pore (24). For a given separation between the particles, it is apparent from Eq. [4] that the magnitude of the interparticle electrostatic repulsion decreases as the zeta potential decreases. As a result, lower magnitudes of the hydrodynamic force existing farther away from the pore are able to overcome the interparticle repulsion and cause flocculation. It is worth noting that the final separation between the particles in the simulations with the dif˚ This value is consistent with ferent zeta potentials was 2.2 A.

FIG. 10. Initial positions of particles specified in the 3-particle simulations: R = 16. The initial separation distance between particles 1–2 and 2–3 is 55 nm.

˚ determined from the location of the primary minimum, 1.7 A, just the interparticle colloidal interaction using Eqs. (3)–(6). In Fig. 10, the initial positions specified in simulations of three particles approaching a pore 3.25 times the particle size are shown. Figure 11 shows the separation distance between the particle on the pore axis and either of the off-axis particles (i.e., referring to Fig. 10, either 1–2 or 2–3) as a function of the position of the particle on the pore axis (particle 2) for two different velocities. The particle flow behavior near the entrance at the two velocities is similar to that shown for the motion of two particles toward a pore. At a farfield velocity corresponding to Re = 0.1, the strong hydrodynamic forces at the pore entrance caused particle flocculation by overcoming the interparticle repulsion. To our knowledge, these simulations are the first to study the interaction between three particles near a pore

FIG. 11. Effect of velocity on particle flow behavior for flow of three particles: α = 3.25, ψ p = 100 mV, deformable particles.

LOW REYNOLDS NUMBER INTERACTIONS

entrance. As mentioned earlier, such simulations are very difficult using the semi-analytic, semi-numerical methods discussed in the Introduction and have not been attempted so far. Our simulations provide interesting insight into the degree of symmetry existing during flow through a pore in a real porous medium. In our simulations of particle interactions at low velocities, it was seen that the particle separation does not decrease to zero because of the strong repulsion. If an aspect ratio of 2 is considered, this result is intriguing because it implies that the particles cannot enter the pore and, therefore, will plug the pore. However, plugging was not observed in the low-velocity experiments of R-F (1999). The absence of plugging in the lowvelocity experiments can be understood by considering the initial conditions specified in our simulations. We have considered a case where both the flow field through the pore and the initial positions of the particles are symmetric about the pore axis. The flow through a pore such as that in the track-etched membrane used in R-F’s (1999) experiments will rarely be symmetric because of two reasons: (a) immediate neighbor pores to any pore are not symmetrically situated, and (2) existence of defects such as noncircular entrance cross-section, overlapping pores and pores not being perpendicular to the plane of the pore entrance. Scanning electron microscope pictures of typical pores in a Nuclepore track-etched membrane are shown in Fig. 12 to support our contention. This membrane is the same as those used in the experiments of R-F (1999). Another important factor not considered in our simulations is the influence of subsequently arriving particles on the motion of the particles near the pore entrance. To demonstrate that asymmetry is essential for nonflocculating particles to flow through the pore without plugging, we have simulated the motion of particles whose initial positions were not symmetric about the pore axis. Figures 13a and 13b

319

FIG. 13. Ability of particles to flow through the pore without plugging: α = 3, ψ p = 100 mV, Re = 0.1; (a) deformable particles, and (b) nondeformable particles.

show simulations with three particles where the on-axis particle (particle 2) is initially placed slightly ahead of the two off-axis particles. An aspect ratio of 3 and a particle zeta potential of 100 mV were specified in the simulations. For both deformable and nondeformable particles, it can be seen that the three particles can flow through the pore without plugging. We performed simulations in which the initial positions of the particles were asymmetric with varying extents of offset of the centers of the particles either in the x direction or in the y direction. An offset in the x direction means that one of the particles is closer to the pore than the other, but both are equidistant from the pore axis; an offset in the y direction means that one of the particles is closer to the pore axis than the other, but both are equidistant from the pore. The objective of these simulations

FIG. 12. Scanning electron microscope pictures of typical pores in a Nuclepore track-etched membrane (nominal pore size 3 µm). For any given pore, the neighboring pores are not symmetrically situated. Also, defects such as overlapping pores and pores not perpendicular to the membrane face are present.

320

RAMACHANDRAN ET AL.

was to identify a critical offset such that for offsets greater than this critical value, particles are able to flow through the pore for velocities less than the critical velocity. In the simulations, when the extent of offset between the particles centers was sufficiently large, the particles did not flocculate and were able to flow through the pore for all Re (similar to results in Fig. 13). For smaller offsets, particles flocculated near the pore for velocities greater than the critical velocity as discussed previously in simulations with symmetric initial positions of particles. We were however unsuccessful in identifying a critical offset for which particles were able to flow through the pore at low velocities but flocculated and plugged the pore at high velocities. The primary reason was the difficulty in resolving the motion of the particles when they are overlapping with the plane of the pore entrance. This difficulty is because of the simulations being in two dimensions. The limitation of two-dimensional simulations can be seen from Fig. 14. In these simulations, the initial positions of the particles are symmetric about the pore axis and the aspect ratio is 2. Figure 14a shows the particles being very near the pore entrance, and Fig. 14b shows the velocity vectors for that configuration of the particles and the pore. One observes fluid circulation above and below particles 1 and 2, respectively. Because

the calculations are in two dimensions, the proximity of the particles to the pore results in the blockage of fluid flow through the pore and hence causes the fluid circulation observed in Fig. 14b. This behavior will not be present in three-dimensional simulations. Also, the accuracy of the numerical simulations can be improved by refining the grid in the regions where the interfaces are interacting. For example, Agresar et al. (1) has studied cell mechanics and cell adhesion using the immersed boundary technique as in this work but on an adaptively refined. Cartesian, unstructured grid. Adaptive refinement allowed the use of sufficiently fine mesh in the neighborhood of the moving interfaces while using a coarse mesh far way from the particles. The inclusion of adaptive refinement may also alleviate to a certain extent the circulation problem when the particles are near the pore. Comparison with experiments. Both two- and three-particle simulations presented so far clearly show that hydrodynamic forces acting on particles in the vicinity of the pore entrance can overcome net interparticle repulsion and cause flocculation. Even though the individual particles are smaller in size than the pore, at sufficiently high flowrates, particles can flocculate at the pore entrance and result in pore blockage. For a given particle size, surface charge, ionic strength, and pH (if H+ is a potential determining ion), it is evident from the simulations that a critical velocity is required to cause particle bridging. This critical velocity, as seen from Fig. 9, is a function of the aspect ratio. These results are in qualitative agreement with the experimental observations of R-F (1999). In the experiments, plugging was observed at Reynolds numbers about an order of magnitude lower than those obtained from the simulations. However, a quantitative comparison between the experiments and the simulations cannot be made for the following reasons: (i) simulations are in two dimensions; (ii) both the pores and the particles used in the experiments have a distribution of sizes (size exclusion however is insignificant); polydisperse particle and pore sizes mean that there exists a distribution of critical velocities for plugging; and (iii) as pores get plugged, the interstitial velocity in the open pores increases which cannot be incorporated in our single-pore simulations. Despite these limitations, the simulations provide insight into the physics of the phenomenon and explain the effect of the various parameters on the flow behavior. CONCLUSIONS

FIG. 14. Limitations of simulations in two dimensions.

The flow behavior of stable colloidal particles approaching the entrance to a pore has been studied using a numerical technique for the case where the particle size is smaller than but of the same order as the pore size. The major advantage of the numerical technique used is its ability to rigorously account for the interparticle and particle–pore hydrodynamic interactions. These interactions are very difficult to determine using the existing analytical and semi-numerical, semi-analytical techniques in microhydrodynamics. The technique was adapted from a front-tracking method developed for studying incompressible,

321

LOW REYNOLDS NUMBER INTERACTIONS

multi-fluid flow by Unverdi and Tryggvason (20). Two- and three-particle motion in two dimensions has been considered in our simulations. The simulations demonstrate that at sufficiently high velocities but in the low Re regime, hydrodynamic forces acting on the particles at the pore entrance can overcome interparticle colloidal repulsion and result in the particles flocculating and plugging the pore. The existence of a critical velocity for plugging under the conditions considered predicted by the simulations has been experimentally demonstrated by Ramachandran and Fogler (15). The predicted dependence of the particle interaction at the pore entrance on the aspect ratio and the particle zeta potential are consistent with expected trends. The simulations also correctly predict the difference in the flow behavior of nondeformable and deformable particles. Another important result is the prediction that for perfect symmetry in both the fluid flow field through the pore and the positions of the particles about the pore axis, pore blockage will occur for any nonzero fluid velocity for aspect ratios equal to or smaller than the number of particles. The absence of plugging observed in the low-velocity experiments of Ramachandran and Fogler (15), with conditions comparable to those specified in the simulations, implies that such symmetry is usually absent in even model experimental systems involving flow through porous media. The simulations provide useful physical insight into the factors affecting the particle motion near the pore entrance even though they are in two dimensions. To obtain quantitative information from the simulations the three-dimensional problem must be solved. The accuracy of the numerical scheme can be improved by increasing the resolution of the portions of the interfaces that are interacting by using adaptive mesh refinement. Furthermore, the particle–pore wall colloidal repulsion and the influence of subsequently arriving particles must be included. APPENDIX

1. Discretization: ˜ n−1 ˜n −A (ρu)n+1 − (ρu)n 3A + 1t 2 ˜Sn+1 + S˜ n + Fn , = −∇h p n+1 + 2

[A4]

where ∇h is the finite difference approximation for ∇. 2. Splitting: ˜n −A ˜ n−1 3A S˜ n ρ n+1 u∗ − (ρu)n + = + Fn , 1t 2 2

[A5]

giving µ ¶ ˜n −A ˜ n−1 S˜ n ρn n 1t 3A n + +F , u = n+1 u + n+1 − ρ ρ 2 2 ∗

[A6]

allowing the explicit evaluation of u∗ . 3. Projection: S˜ n+1 (ρu)n+1 − ρ n+1 u∗ = −∇h p n+1 + . 1t 2

[A7]

Taking the divergence of [A7] and rearranging, we get the Poisson equation for pressure at (n + 1): 1

∇h ·

ρ n+1

∇h p n+1 =

∇ h · u∗ 1 S˜ n+1 + ∇h · n+1 1t ρ 2

[A8]

An equation for un+1 is got from [A7]:

The finite-difference scheme used for solving the Navier– Stokes equations is described here. The scheme is implicit in time and is implemented on a structured, staggered, Eulerian grid shown in Fig. 3. The Navier–Stokes equation to be solved is ∂ ˜ = −∇ p + S˜ + F, (ρu) + A ∂t

in Peyret and Taylor (13, p. 164). In construction, this scheme closely follows the well-known explicit projection method:

[A1]

n+1

u

¶ µ 1t S˜ n+1 n+1 = u + n+1 −∇h p + . ρ 2 ∗

[A9]

This equation is however implicit in un+1 . un+1 and p n+1 are obtained by iteratively solving Eqs. [A8] and [A9] using the successive overrelaxation technique. The boundary conditions used for the velocities are based on those specified in (10):

where ˜ = ∇ · (ρuu) A

[A2]

S˜ = ∇ · µ[∇u + (∇u)T ]

[A3]

and F represents surface tension, gravity, and colloidal forces. Equation [A1] is supplemented by the incompressibility condition ∇ · u = 0. The scheme used is the implicit, second-order accurate Adams–Bashforth/Crank–Nicholson scheme outlined

∂u = 0, v = 0. ∂x ∂u Boundaries 2 and 3: = 0, v = 0. ∂y Boundaries 4–7: u = v = 0. ∂u ∂v Boundary 8: = = 0. ∂x ∂x Boundary 1: u = U,

[A10]

322

RAMACHANDRAN ET AL.

It is ensured that p satisfies the Neumann boundary condition (13) ¯ µ n+1 n+1 ¶¯ ρ (u − u∗ ) 1 ˜ n+1 ¯¯ ∂ p n+1 ¯¯ = − + S ¯ · N, [A11] ∂ N ¯0 1t 2 0 _

_

where N is the unit outward normal to the boundary 0. In the simulations, the magnitudes of distances and velocities at most of the grid points were small because of the micrometer length scale and the low Reynolds number flow conditions involved. We therefore used a de-dimensionalized form of the equation wherein proper scaling lead to distances and velocities being O(1). The de-dimensionalized form of the Navier–Stokes equatins and the incompressibility condition is ∂ ˜0 (ρu0 ) + A ∂t 0

µ ¶ µ ¶ 2 ρp ˜ 0 a = −∇ p + S + F; ∇ 0 · u0 = 0. Re ρs ρs u 2c 0 0

[A12]

The quantities used for scaling are xi0 =

xi ui uc ; u i0 = ; t 0 = t; a uc a

p0 =

p 1 ; ∇ 0 = ∇, rs u 2c a [A13]

where u c is a characteristic velocity which is chosen to be U , the x velocity far from the pore. REFERENCES 1. Agresar, G., Lindermann, J. J., Tryggvason, T., and Powell, K. G., J. Comp. Phys. 143, 346 (1998).

2. Bhattacharjee, S., and Sharma, A., J. Colloid Interface Sci. 171, 288 (1995). 3. Brady, J. F., and Bossis, G., Ann. Rev. Fluid Mech. 20, 111 (1988). 4. Clift, R., Grace, J. R., and Weber, M. E., “Bubbles, Drops, and Particles.” Academic Press, New York, 1978. 5. Bowen, W. R., and Sharif, A. O., J. Colloid Interface Sci. 187, 363 (1997). 6. Feke, D. L., Prabhu, N. D., Amin, J. A., Jr., and Amin, J. A., II, J. Phys. Chem. 88, 5735 (1984). 7. Hiemenz, P. C., and Rajagopalan, R., “Principles of Colloid and Surface Chemistry,” 3rd ed. Dekker, 1997. 8. Wu, W., Weinbaum, S., and Acrivos, A., J. Fluid Mech. 243, 489 (1992). 9. Kao, J., Wang, Y., Pfeffer, R., and Weinbaum, S., J. Collid Interface Sci. 121, 543 (1988). 10. Manton, M. J., Atmos. Environ. 12, 1669 (1978). 11. Yan, Z. Y., Weinbaum, S., and Pfeffer, R., J. Fluid Mech. 162, 415 (1986). 12. Peskin, C., J. Comp. Phys. 25, 220 (1977). 13. Peyret, R., and Taylor, T. D., “Computational Methods for Fluid Flow.” Springer-Verlag, Berlin, 1983. 14. Ramachandran, V., and Fogler, H. S., Langmuir 14, 4435 (1998). 15. Ramachandran, V., and Fogler, H. S., J. Fluid Mech. 385, 129 (1999). 16. Russel, W. B., Saville, D. A., and Schowalter, W. R., “Colloidal Dispersions.” Cambridge Univ. Press, Cambridge, 1989. 17. Schenkel, J. H., and Kitchener, J. A., Trans. Faraday Soc. 56, 161 (1960). 18. Israelachvili, J. N., “Intermolecular and Surface Forces—With Application to Colloidal and Biological Systems. Academic Press, San Diego, 1985. 19. Gregory, J., J. Colloid Interface Sci. 83, 131 (1981). 20. Unverdi, S. O., and Tryggvason, G., J. Comp. Phys. 100, 25 (1992). 21. Vitthal, S., and Sharma, M. M., J. Colloid Interface Sci. 153, 314 (1992). 22. Glendinning, A. B., and Russel, W. B., J. Colloid Interface Sci. 93, 95 (1983). 23. Weinbaum, S., Ganatos, P., and Yan, Z., Annu. Rev. Fluid. Mech. 22, 275 (1990). 24. Yan, Z., Weinbaum, S., Ganatos, P., and Pfeffer, R., J. Fluid Mech. 174, 39 (1987).