Langmuir 2006, 22, 9879-9893
9879
Particle Deposition onto Charge Heterogeneous Surfaces: Convection-Diffusion-Migration Model Neda Nazemifard,† Jacob H. Masliyah,‡ and Subir Bhattacharjee*,† Department of Mechanical Engineering, UniVersity of Alberta, Edmonton, Alberta T6G 2G8, Canada, and Department of Chemical & Materials Engineering, UniVersity of Alberta, Edmonton, Alberta T6G 2G6, Canada ReceiVed June 13, 2006. In Final Form: September 10, 2006 The role of surface charge heterogeneity of planar collectors on particle deposition and distribution was investigated in the vicinity of a heterogeneous surface for a radial impinging jet flow geometry. The charge heterogeneity was modeled as concentric circular stripes bearing different surface charges. Particle deposition was studied employing the Eulerian approach (convection-diffusion-migration equation). It was observed that when the collector was completely unfavorable, the presence of small amounts of charge heterogeneity in the form of a small fraction of favorably charged stripes enhanced the deposition rate substantially. In contrast, when the collector was completely favorable, the presence of small amounts of charge heterogeneity in the form of a small fraction of unfavorably charged stripes did not affect the particle deposition rate significantly.
Introduction Particle deposition is of great importance in a variety of natural and industrial processes.1-3 Many studies in the literature have theoretically investigated particle deposition over homogeneous surfaces.3-5 However, natural or synthetic solid surfaces often exhibit both physical and chemical heterogeneities.6-13 Even particle deposition on clean and optically flat glass surfaces exhibits local variations owing to the heterogeneity of the glass.14,15 Chemical heterogeneity results in an uneven or heterogeneous surface charge distribution at various length scales.16-18 Numerous large-scale studies on particle deposition and capture onto charge heterogeneous substrates, generally pertaining to transport in * Corresponding author. Telephone: (780)492-6712. Fax: (780)492-2200. E-mail:
[email protected]. † Department of Mechanical Engineering. ‡ Department of Chemical & Materials Engineering. (1) Hirtzel, C. S.; Rajagopalan, R. Chem. Eng. Commun. 1985, 33, 301. (2) Adamczyk, Z. Colloids Surf. 1989, 39, 1. (3) Elimelech, M.; Gregory, J.; Jia, X.; Williams, R. A. Particle Deposition and Aggregation: Measurement, Modelling, and Simulation; Colloid and Surface Engineering Series; Butterworth-Heinemann: Woburn, MA, 1995. (4) van de Ven, T. G. M. Colloidal Hydrodynamics; Academic Press: San Diego, 1989. (5) Yang, J.; Bos, R.; Poortinga, A.; Wit, P. J.; Belder, G. F.; Busscher, H. J. Langmuir 1999, 15, 4671. (6) Vaidyanathan, R.; Tien, C. Chem. Eng. Sci. 1991, 46, 967. (7) Kihira, H.; Ryde, N.; Matijevic, E. J. Chem. Soc.-Faraday Trans. 1992, 88, 2379. (8) Song, L. F.; Johnson, P. R.; Elimelech, M. EnViron. Sci. Technol. 1994, 28, 1164. (9) Johnson, P. R.; Sun, N.; Elimelech, M. EnViron. Sci. Technol. 1996, 30, 3284. (10) Chen, J. Y.; Ko, C. H.; Bhattacharjee, S.; Elimelech, M. Colloids Surf. A 2001, 191, 3. (11) Shellenberger, K.; Logan, B. E. EnViron. Sci. Technol. 2002, 36, 184. (12) Bhattacharjee, S.; Ryan, J. N.; Elimelech, M. J. Contam. Hydrol. 2002, 57, 161. (13) Loveland, J. P.; Bhattacharjee, S.; Ryan, J. N.; Elimelech, M. J. Contam. Hydrol. 2003, 65, 161. (14) Wit, P. J.; Busscher, H. J. J. Colloid Interface Sci. 1998, 208, 351. (15) Luthi, Y.; Ricka, J.; Borkovec, M. J. Colloid Interface Sci. 2005, 206, 314. (16) Richmond, P. J. Chem. Soc.-Faraday Trans. 1 1974, 2, 229. (17) Khachatourian, A. V. M.; Wistrom, A. O. J. Phys. Chem. B 1998, 102, 2483. (18) Vreeker, R.; Kuin, A. J.; Denboer, D. C.; Hoekstra, L. L.; Agterof, W. G. M. J. Colloid Interface Sci. 1992, 154, 138.
granular porous media, have been conducted.3,9,10,12,13 Most of the theoretical studies on this subject resort to a continuum-type approach in describing the charge heterogeneity of the substrate. A commonly adopted technique is to define two types of surface charge locations on a given collector (for instance, positive and negative), assigning the surface area fraction occupied by one type of charge and using a two-site averaging process generally referred to as the patchwise heterogeneity model.8,9 More recently, it has become increasingly clear that while the patchwise heterogeneity models provide reasonably accurate description of macroscopic charge heterogeneity, they are not very accurate when the charge heterogeneous patches have a comparable length scale to the particle size.19,20 This implies that the simple averaging methods, such as patchwise heterogeneity models, may not be adequate to predict the particle deposition in the presence of micro-scale charge heterogeneity. It is therefore of great importance from both a theoretical and a practical perspective to fundamentally explore the effects of micro-scale charge heterogeneity on governing mechanisms of the particle deposition process. Naturally occurring deposition substrates generally contain surface charge heterogeneity that are randomly distributed, of arbitrary geometrical shapes, and have widely varying chemical properties.13,21 Consequently, it becomes extremely difficult to systematically study the influence of charge heterogeneity on the particle deposition behavior onto such substrates. In this context, it might be pertinent to systematically create such heterogeneity by chemically patterning homogeneous model substrates.19 Studying particle deposition onto such a model substrate, where the heterogeneity is artificially created and its nature is known accurately, can lead to considerable insight regarding how the deposition behavior is influenced by the presence of surface charge heterogeneity. If the distribution of the heterogeneous patches is known a priori, elucidation of their influence on particle deposition becomes more tractable.20 In our previous study,20 a theoretical model based on the Lagrangian approach (trajectory analysis) was developed to (19) Elimelech, M.; Chen, J. Y.; Kuznar, Z. A. Langmuir 2003, 19, 6594. (20) Nazemifard, N.; Masliyah, J. H.; Bhattacharjee, S. J. Colloid Interface Sci. 2006, 293, 1. (21) Ryan, J. N.; Elimelech, M. Colloids Surf. A 1996, 107, 1.
10.1021/la061702q CCC: $33.50 © 2006 American Chemical Society Published on Web 10/28/2006
9880 Langmuir, Vol. 22, No. 24, 2006
Nazemifard et al.
investigate the effect of microscopic surface charge heterogeneity on particle deposition inside a radial impinging jet system. According to the trajectory model, there was a region at the leading edge of each favorable stripe that was inaccessible for particle deposition. However, the particle deposition rates predicted by trajectory analysis were significantly larger than those predicted by the patchwise heterogeneity model.9 In this context, one might construe that the trajectory analysis may not be able to provide a holistic picture about the distribution of the particles in the vicinity of the charge heterogeneous patches. In the present study, to achieve a better understanding of the particle deposition behavior in the presence of micro-scale charge heterogeneity, a convection-diffusion-migration model has been developed. This model provides the particle concentration distribution and particle deposition rate in the vicinity of the micro-scale charge heterogeneity in the impinging jet system. The Eulerian model is used to assess the accuracy of the trajectory analysis.20 The Eulerian model allows incorporation of the diffusion dominated deposition process, which was not possible in the trajectory modeling. The predictions of the deposition rate from the present model are compared against those from the Lagrangian model to assess the deposition process onto charge heterogeneous collectors.
Mathematical Model Impinging Jet Flow Geometry. A schematic representation of the radial impinging jet system is illustrated in Figure 1a. A colloidal suspension is introduced to the system through a circular jet of radius (Rjet) with an average velocity of U∞. The fluid impinges vertically on the collector surface and flows radially outward from the impinging jet flow domain. The separation distance between the collector surface and the impinging jet tip is denoted by L. The geometry of the system described here ensures radial symmetry and allows the use of an axisymmetric cylindrical coordinate system. The origin of the cylindrical coordinate system is denoted by O and is located at the stagnation point. Figure 1b schematically depicts the velocity components and the forces acting on a colloidal particle in the flow domain. The particle radius is ap, z denotes the normal distance of the particle center from the collector surface, and s denotes the normal surfaceto-surface distance between the particle and the collector. In the cylindrical coordinate system, the radial and axial (normal to collector surface) particle velocity components are given by Vr and Vz, respectively. Three types of non-hydrodynamic forces are considered, namely, gravity (Fg), van der Waals attraction (Fvdw) between the particle and the collector, and electric doublelayer forces (Fedl), which may be attractive or repulsive, depending on the particle and collector surface potential or charge density. We have considered the deposition to occur against a gravitational field to be consistent with similar studies against which our numerical results are compared. The model can be employed to assess deposition assisted by gravity as well. The dimensions of the system considered are Rjet ) 1 mm, L ) 2 mm, ap ) 1 µm. The fluid velocity field is computed over the domain EFGO, where the radial distance OG is 5Rjet. Convection-Diffusion-Migration Equation. The steadystate transport of particles from a flowing suspension toward a collector surface is governed by the mass conservation equation:4,5,22
∇ B ‚bj ) 0
(1)
Figure 1. Schematic representations of (a) a radial impinging jet system. A colloidal suspension is introduced to the system through a circular jet of radius Rjet, with an average velocity of U∞. The two arrows represent the fluid streamlines. The stagnation point is shown by O, which is the center of the cylindrical coordinate system. The distance between the collector surface and the nozzle tip is represented by L. The box EFGO denotes the computational domain for solving the fluid velocity field. (b) The velocity components and the forces acting on a colloidal particle in the flow domain. The particle radius is ap, z denotes the normal distance of the particle center from the collector surface, and s denotes the normal surface to surface distance between the particle and the collector surface. Particle velocity components are given by Vr and Vz, respectively. Three types of non-hydrodynamic forces are considered, namely, gravity (Fg), van der Waals attraction (Fvdw) between the particle and the collector, and electric double-layer forces (Fedl), which may be attractive or repulsive, depending on the particle and collector surface potential or charge density.
The particle flux (jb) is obtained by considering contributions from convection, diffusion, and migration under the influence of all other non-hydrodynamic forces exerted on the particle. Considering these, the particle flux can be defined as
bj ) b V c - D‚∇ Bc +
c D‚F B kbT
(2)
where b V is the particle velocity, c is the particle concentration, D is the diffusion coefficient tensor, kb is the Boltzmann constant, T is the absolute temperature, and B F is the total force exerted on the particle. To solve the convection-diffusionmigration equation (eq 1), the particle velocity, the diffusion tensor, and the total external forces need to be defined. To determine the particle velocity vector (V b), an accurate evaluation of the fluid velocity field (u b) is necessary. The fluid velocity field inside the impinging jet system was obtained numerically from a solution of the Navier-Stokes and the continuity equations over the entire impinging jet geometry (shaded domain EFGO in Figure 1a) using the finite element method. Detailed description of the numerical scheme and the velocity field calculations are available elsewhere.20,23 A (22) Masliyah, J. H.; Bhattacharjee, S. Electrokinetic and Colloid Transport Phenomena; John Wiley and Sons: New York: 2006. (23) Nazemifard, N. Particle deposition onto charge heterogeneous surfaces in a radial impinging jet system. M.Sc. Thesis, 2006.
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9881
functions and h is the dimensionless surface to surface distance between the particle and the substrate (h ) s/ap) where s is the separation distance between the particle and the substrate as depicted in Figure 1b. To account for the hydrodynamic interactions between the particle and the collector,4,27 the diffusion tensor is used in eq 2. The diffusion tensor can be formulated for the impinging radial flow geometry by incorporating the effects of hydrodynamic interactions as follows
D ) D∞
Figure 2. Velocity field streamlines in the impinging jet flow system, within the domain denoted by EFGO in Figure 1a. These streamlines correspond to Re of 100. Hatched lines represent the collector surface and the jet wall.
streamline depiction of a typical flow field inside the impinging jet geometry is shown in Figure 2. These results are obtained for Re ) 100, which corresponds to an average flow velocity of U∞ ) 0.1 m/s. The flow distribution in the impinging jet geometry can be characterized by the dominant role of normal fluid velocity near the stagnation point (O) with a more pronounced role of tangential (radial) fluid velocity at larger radial distances from the stagnation point. It is worth noting here that analytical expressions for the velocity field inside the impinging jet flow system are available under some limiting conditions, such as, in the vicinity of the stagnation point. These analytical expressions were employed to validate the accuracy of the numerical results for fluid velocity. Near the stagnation point, analytical expressions for the radial and normal components of the fluid velocity (ur and uz) respectively, take the following forms:22,24
ur ) Rsrz
(3a)
uz ) -Rsz2
(3b)
where Rs is a function dependent on the Reynolds number (Re ) U∞Rjet/ν, where ν is the kinematic viscosity) and the geometry of the system, given by
Rs )
R j νRe Rjet3
(4)
The term R j is a function characterizing the intensity of stagnation flow and depends on the Reynolds number and flow geometry.22,25 In the vicinity of the collector surface, the velocity of the particle induced by fluid motion is different from that of the fluid due to the effect of hydrodynamic interactions. Following Spielman and Fitzpatrick,26 the relationship between the fluid and particle velocities is given by
Vr ) ur f3(h)
(5a)
Vz ) uz f1(h)f2(h)
(5b)
where f1(h) to f3(h) are the universal hydrodynamic correction (24) Dabros, T.; van de Ven, T. G. M. Colloid Polym. Sci. 1983, 261, 694. (25) Adamczyk, Z.; Zembala, M.; Siwek, B.; Czarnecki, J. J. Colloid Interface Sci. 1986, 110, 188. (26) Spielman, L. A.; Fitzpatrick, J. A. J. Colloid Interface Sci. 1973, 42, 607.
[
f4(h) 0 0 f1(h)
]
(6)
where D∞ is the Stokes-Einstein diffusion coefficient, while f1 and f4 are the universal hydrodynamic correction functions.28-31 Simplified expressions for the universal hydrodynamic correction functions22 were used in the calculations. The non-hydrodynamic forces acting on the particle are due to gravity (Fg) as well as the van der Waals (Fvdw) and electrostatic double-layer (Fedl) interactions between the particle and the collector. We assume no lateral non-hydrodynamic forces on the particle, Fr ) 0. The total normal force acting on the particle is Fz ) Fg + Fvdw + Fedl. For the particles suspended in a liquid, the net gravitational force can be expressed as
4 b B Fg ) πap3∆Fg 3
(7)
where ∆F is the difference between the particle and fluid densities. From Figure 1b, it is evident that the gravitational force (eq 7) acts along the positive z direction for positive values of ∆F, implying that the particle deposition is opposed by gravity in subsequent calculations, where the particle density is greater than the fluid density. The gravity force can be nondimesionalized with respect to the Brownian force (FBr ) kbT/ap), yielding
F hg )
Fgap 4π ∆Fgap4 ) kbT 3 kbT
(8)
where F h g is called gravity number. An expression for the nonretarded van der Waals force between a sphere and an infinite planar surface based on Hamaker’s approach, and Derjaguin’s approximation is given by32
Fvdw ) -
A 1 6ap h2
(9)
where A is the effective Hamaker constant between the particle and the collector surface in the liquid medium and h is the dimensionless surface to surface separation distance between the particle and the collector surface. The expression for van der Waals interaction can be nondimensionalized with respect to Brownian force as (27) Honig, E. P.; Roeberson, G. J.; Wiersema, P. H. J. Colloid Interface Sci. 1971, 36, 97. (28) Goldman, A. J.; Cox, R. G.; Brenner, H. Chem. Eng. Sci. 1967, 22, 637. (29) Goldman, A. J.; Cox, R. G.; Brenner, H. Chem. Eng. Sci. 1967, 22, 653. (30) Goren, S. L. J. Fluid Mech. 1970, 41, 619. (31) Goren, S. L.; O’Neill, M. E. Chem. Eng. Sci. 1971, 26, 325. (32) Suzuki, A.; Ho, N. F. H.; Higuchi, W. I. J. Colloid Interface Sci. 1969, 29, 552.
9882 Langmuir, Vol. 22, No. 24, 2006
F h vdw )
Nazemifard et al.
Fvdwap 1 ) - Ad 2 kbT h
(10)
where Ad ) A/6kbT is the adhesion number and indicates the strength of van der Waals interaction. The van der Waals force between the particle and the collector is always attractive in our calculations. Assuming constant surface potentials, symmetric (z:z) electrolytes, and linearized Poisson-Boltzmann equation to be applicable, the well-known Hogg, Healy, and Fuerstenau (HHF) expression33 is employed for the electric double-layer (EDL) interaction force. The force is expressed as
( ) [
Fedl ) 4π0(κap)
kbT 2 exp(-κaph) ΨpΨs ze 1 + exp(-κaph)
]
(Ψp - Ψs)2 exp(-2κaph) (11) 2ΨpΨs 1 - exp(-2κaph)
Here, Ψp and Ψs are the scaled surface potentials (Ψ ) zeψ/kbT) of the particle and the substrate, respectively, and κ is the inverse Debye length.34 Equation 11 can be nondimensionalized with respect to the Brownian force as
F h edl )
Fedlap ) kbT Dl(κap)
[
]
exp(-κaph)
exp(-2κaph) (12) -Da 1 + exp(-κaph) 1 - exp(-2κaph)
where Dl is the electrostatic double-layer parameter defined as
Dl )
4π0apkbTΨpΨs z2e2
(13)
and Da is the double-layer asymmetry parameter:
Da )
(Ψp - Ψs)2 2ΨpΨs
(14)
For interacting surfaces with same surface potential, Da assumes a value of zero. The EDL interaction between the collector and the particle can be either attractive or repulsive. Patterned Charge Heterogeneity. For a charge heterogeneous collector, the electric double-layer interaction between the particle and the collector will be a function of spatial location. Such spatial variation of the EDL force will depend on how the surface charge heterogeneity is distributed on the collector surface. In this study, the collector surface charge heterogeneity is modeled as concentric stripes of microscopic dimensions bearing positive and negative charges as shown in Figure 3. Each negative and positive stripe has a width wn and wp, respectively. The total width of a pair of consecutive negative and positive stripes is defined as the pitch or bandwidth, p. For a negatively charged particle, the ratio of positive stripe width (wp) to pitch (p) specifies the fraction of the collector surface area that attracts the particle and, hence, is favorable for deposition. The favorable area fraction of the collector is denoted by λpatterned ) Sf/S where Sf is the oppositely charged area and S is the total collector surface area. The unfavorable area fraction (33) Hogg, R.; Healy, T. W.; Fuerstenau, D. E. Trans. Faraday Soc. 1966, 62, 1638. (34) Hunter, R. J. Zeta Potential in Colloid Science: Principles and Applications; Academic Press: London, 1981.
Figure 3. Schematic representation of the modeled surface charge heterogeneity. The collector consists of concentric circular stripes with alternate negative and positive surface potentials with widths of wn and wp, respectively. The total width of a negative and a positive stripe gives the pitch, p. Nb is the band number, a number assigned to each pair of negative and positive stripe along the radial direction on the collector surface. rNb is the radial position of the leading edge of a pair of negative and positive stripe specified by Nb.
is given by (1 - λpatterned). If sufficiently large number of stripes are considered simultaneously, which occurs as one moves radially away from the stagnation point, the following approximation for λpatterned becomes pertinent:
λpatterned =
wp p
(15)
The consecutive bands are denoted by a counter, Nb, in increasing order of its radial position. Each band spans a consecutive negative and positive stripe (bandwidth ) wn + wp ) p) as shown in Figure 3. In this study, the innermost stripe at the stagnation point is always negatively charged and, hence, is unfavorable with respect to deposition of the negatively charged particle. The charge heterogeneity is incorporated in the model as a variation of Ψs with r in a periodic manner. Over the unfavorable stripe, the surface potential is considered to be Ψs, n. At the boundary between the negative and positive stripes, the surface potential is varied smoothly (using sigmoidal smoothing functions) from its value over the negative stripe, Ψs, n, to its value over the positive stripe, Ψs, p while assuming a value of zero at the boundary of the two stripes.23 The smoothing function was used to prevent discontinuity of the surface potential on the collector. The use of eq 11 in conjunction with the above periodic variation of Ψs constitutes a highly simplified model for incorporating surface charge heterogeneity effects on the EDL interaction force. According to this model, the EDL force is solely dictated by the surface potentials of the particle and region of the collector directly facing it (distance of closest approach). A more rigorous approach for obtaining the EDL force is to numerically solve the Poisson-Boltzmann equation for a
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9883
Table 1. Dimensionless Groups and Parameters Used in the Scaled Convection-Diffusion-Migration Equation dimensionless group
expression
scaled distance scaled velocity scaled concentration scaled particle flux (Sherwood number) scaled normal force Peclet number
jr ) r/ap; jz ) z/ap ujr ) ur/U∞; ujz ) uz/U∞ cj ) c/c∞ jj ) Sh ) jap/(c∞D∞) F h z ) Fzap/(kbT) Pe ) U∞ap/D∞
heterogeneous surface.35 However, based on that study and our numerical calculations for the present geometry, we conclude that the HHF expression provides a reasonable estimate of the normal component of the EDL interaction force when κap > 10. Our use of the HHF expression is primarily dictated by its simplicity and requirement of considerably less computational resources compared to a direct numerical solution of the PoissonBoltzmann equation. The influence of the sigmoidal variation of the surface potential at the edges of the positive and negative stripes was also investigated. It was found to have minimal influence on the calculated EDL force (less than 2% variation) as long as the half-width of the sigmoid was < 1 µm. A halfwidth of 0.1 µm was used in the calculations. Notwithstanding the above confirmations, the HHF approximation might be an oversimplification of the actual EDL interactions between a charge heterogeneous substrate and the particle for smaller values of κap < 10. This methodology also ignores the exact lateral variations of the EDL forces. However, as long as the individual positive or negative stripe width (wp or wn) is larger than the particle size, the particle is sufficiently close to the substrate, and the electric double-layer interactions are sufficiently screened (large κap), use of the above approximation is expected to be reasonably accurate. Particle Transport Equation and Boundary Conditions. The nondimensional form of the particle transport equation (eq 1), written in an axi-symmetric cylindrical coordinate system is
∇ B .B jj )
∂jjz 1∂ (rjjj r) + )0 jr ∂rj ∂zj
∂cj jj r ) ujr f3Pecj - f4 ∂rj jj z ) ujz f1 f2Pecj - f1
∂cj + f1F h zcj ∂zj
(16a)
(16b)
(16c)
where eqs 16b and 16c provide the scaled radial and normal particle fluxes, respectively. The definitions of the pertinent dimensionless groups in the above equations are given in Table 1. F h z is the scaled non-hydrodynamic force along the normal direction. The overall computational geometry is depicted in Figure 1a as EFGO, which was used to numerically calculate the fluid velocity field. However, noting that the particle concentration remains virtually unchanged from its bulk value except in the concentration boundary layer, we use a smaller computation domain to solve the convection-diffusion equation. As suggested in the literature,2 for micron-sized particles, the diffusion boundary layer thickness is of order of 1 to 10 µm. Accordingly, the outer boundary for the particle transport equation is fixed at a distance of z ) 50 µm from the collector surface, which is sufficiently large for the assumption of bulk concentration to be valid at this boundary. The radial limit of the computational domain is (35) Das, P. K.; Bhattacharjee, S. Langmuir 2005, 21, 4755.
Figure 4. Schematic presentation of the computational domain used to solve the convection-diffusion equation. Equation 16 was first solved inside the larger domain ABPO to obtain a coarse estimate of the concentration distribution. After calculating the coarse concentration distribution, a smaller computation domain is defined as CDPO to refine the concentration distribution near the collector surface. The concentration along the line CD was taken from the coarse solution and applied as a Dirichlet boundary condition for the smaller domain CDPO. The line CD in the smaller domain is located at a distance of z ) 2 µm from the collector surface.
chosen to be the jet radius (Rjet). This computational domain is schematically shown in Figure 4 as ABPO. The figure also depicts a smaller (shaded) domain CDPO, which is used to refine the concentration distributions near the collector surface, as will be described later. The boundary conditions for eq 16 over the domain ABPO are given by
cj ) 0
∂Ω ∈OP
b n ‚B jj ) 0
∂Ω ∈OA
cj ) 1
∂Ω ∈AB
b n ‚∇ B cj ) 0
∂Ω ∈BP
(17a) (17b) (17c) (17d)
where b n represents a unit normal vector to the surface (boundary) under consideration. The boundary condition on OP is the perfect sink condition. Although we denote this boundary to be coincident with the collector surface in Figure 4, it is actually placed at a finite cutoff distance, s ) δm, from the collector surface to prevent the divergence of van der Waals force (eq 9) at contact.2,3,22,26,36 We chose this cutoff distance to be 1 nm from the collector surface in all calculations. The boundary condition on OA stems from the axial symmetry, while the condition on BP states that there is no axial concentration gradient, which implies that the flux normal to this boundary is purely convective. On AB, the particle concentration is assumed to attain the bulk concentration value. Particle Deposition Rate. The dimensionless particle deposition rate is represented by the Sherwood number evaluated at the collector surface:
Sh(r) )
apjz(r) ) (jjz)z)ap+δm D∞c∞
(18)
Here, jjz is evaluated at z ) ap + δm. Equation 18 gives the local Sherwood number at any radial position on the collector surface. The average Sherwood number over the entire collector can be expressed as
Sh )
∫
1 Sh(r) dS S
(19)
where S is the collector surface area. Considering the geometry (36) Prieve, D. C.; Ruckenstein, E. AIChE J. 1974, 20, 1178.
9884 Langmuir, Vol. 22, No. 24, 2006
Nazemifard et al.
of the impinging jet system, the above equation can be explicitly written in the cylindrical coordinate system as
Sh )
1 2 π(r2 - r12)
∫r
r2
2πrSh(r) dr )
1
2 2 (r2 - r12)
∫r
r2
rSh(r) dr
1
(20)
Assigning different values for the limits r1 and r2 in eq 20 yields different types of averages. For instance, an average Sherwood number over a collector of radius Rc can be defined as
Sh )
2 Rc2
∫0R rSh(r) dr c
(21)
For heterogeneous collectors, one can define a local average Sherwood number on an individual band (comprising an unfavorable and a favorable stripe) by assigning r1 ) rNb and r2 ) rNb + p, where rNb is the radial position of the leading edge of a pair of unfavorable and favorable stripes with band number Nb and pitch p. Noting that rNb ) (Nb - 1)p, one can define a local average Sherwood number for a single heterogeneous band on the collector surface, Shp(Nb), as
2 Shp(Nb) ) [(rNb + p)2 - rNb2]
∫r
rNb+p
rSh(r) dr )
Nb
p rSh(r) dr ∫(NN -1)p
2 (2Nb - 1)p2
b
(22)
b
This average Sherwood number depends on the band number and, hence, the radial position over the collector. In addition, two other types of dimensionless deposition rates can be defined pertaining to the homogeneous unfavorable and favorable collectors, namely, Shu and Shf, respectively. The different types of Sherwood numbers used in this work are defined in the Nomenclature. These expressions for average Sherwood number are useful in comparing the deposition behavior of homogeneous and heterogeneous collectors.
Numerical Methodology and Model Validation Numerical Solution. Numerical solution of the convectiondiffusion equation (eq 16) subject to the boundary conditions (eq 17) was obtained using finite element analysis. The solution methodology in this study was implemented using the commercially available software COMSOL 3.2 (COMSOL, Inc., USA). The computational domain was discretized using quadrilateral mesh. Lagrangian elements of second order were employed. A highly refined mesh was necessary for the regions having high concentration gradients. Since colloidal interactions lead to shortrange forces, a significant particle concentration gradient exists in the vicinity of the collector surface.26,37 This necessitates the use of extremely small elements with length scales even smaller than Debye length close to the surface.8 For the case of heterogeneous collectors consisting of alternate favorable and unfavorable stripes, the number of elements required to solve the convection-diffusion equation is even higher as compared to a homogeneous collector. For a heterogeneous collector, the EDL force changes over the consecutive stripes. This causes large variations of concentration with r near the edges of each stripe. Consequently, for heterogeneous collectors, in addition to small element size near the collector surface, a (37) Ruckenstein, E.; Prieve, D. C. J. Chem. Soc.-Faraday Trans. 1 1973, 69, 1522.
Table 2. Physical and Chemical Properties of the System Used in the Simulations property
value
particle radius, ap particle density, Fp fluid density, Ff fluid viscosity, µf fluid inlet velocity, U∞ gravitational acceleration, g Reynolds number, Re temperature, T Boltzmann constant, kb electronic charge, e permittivity of vacuum, 0 dielectric constant, scaled particle surface potential, Ψp scaled surface potential of favorable stripes, Ψs,p scaled surface potential of unfavorable stripes, Ψs,n solution ionic strength, I Hamaker constant, A pitch, p valence of ion, z
0.1-2 µm 2300 kg/m3 1000 kg/m3 1.0 × 10-3 N‚s/m2 0.05-0.2 m/s 9.81 m/s2 50-200 298 K 1.38 × 10-23 J/K 1.6 × 10-19 C 8.85 × 10-12 C/Vm 78.54 -1, -2 +1, +2 -1, -2 10-3, 10-3.5 M 10-20 J 4-20 µm 1
high mesh density is also required for the regions where the surface potential changes its sign. This increases the computational loads significantly. To overcome this problem, a multi-geometry procedure was adopted. First, the fluid velocity field was obtained by solving the Navier-Stokes equation over the entire impinging jet geometry (EFGO in Figure 1a). The fluid velocity field was then employed to solve the convective-diffusion equation (eq 16) in the domain ABPO (Figure 4) to obtain a coarse estimate of the concentration distribution. The computational domain ABPO was discretized into 20 000 elements with increasing mesh density near the collector surface. After calculating the coarse concentration distribution, we define another computational domain CDPO (Figure 4). The line CD in this smaller domain was located at a distance of z ) 2 µm from the collector surface. Once again, approximately 40 000 elements were used in this smaller domain. The particle concentration along the line CD was taken from the coarse solution and applied as a Dirichlet boundary condition for the smaller domain. Following this, the convection-diffusion equation was solved in the smaller domain CDPO, keeping all the other boundary conditions unchanged. Shrinking the computational domain from ABPO to CDPO increases the resolution of the geometry near the collector surface. The process of solving the same set of equations inside two different geometries employed in our simulation can be categorized as a class of the general multi-grid methods in numerical simulations. The multi-geometry environment of COMSOL allows us to create these geometries, while extruding or interpolating the results between them. This method minimized memory consumption and accelerated the solution time significantly. Validation of Numerical Results. To investigate the accuracy of the finite element calculations, the numerical results are compared with several analytical and numerical results for homogeneous collectors available in the literature. The operating and physicochemical properties of the modeled system used in this work are listed in Table 2. Available studies in the literature mainly deal with solution of the convection-diffusion equation in the stagnation point region of homogeneous collectors, where the analytical expression for fluid velocity (eq 3) is valid. Incorporating these expressions
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9885
in the convection-diffusion-migration equation, one can simplify eq 16 to
∇ B .B jj )
∂jjz 1∂ (rj jj ) + )0 jr ∂rj r ∂zj
(23a)
∂cj 1 jj r ) f3Perj jz cj - f4 2 ∂rj
(23b)
∂cj 1 h zcj jj z ) f1 f2Pes jz2 cj - f1 + f1F 2 ∂zj
(23c)
where Pes is the stagnation flow Peclet number defined as
Pes )
2Rsap3 D∞
(24)
One of the limiting cases for deposition in the stagnation point region is the diffusion dominated deposition in absence of hydrodynamic and field forces as obtained by Levich.38 For diffusion-controlled deposition (i.e., Pes , 1), the analytical expression for the average Sherwood number near the stagnation point takes the form of38
Sh ) 0.616Pes1/3
(25)
Our numerical code can replicate the assumptions of the Levich analysis by setting F hg ) F h vdw ) F h edl ) 0 and f1 ) f2 ) f3 ) f4 ) 1. The numerical model is first tested by computing Sh in the stagnation point region in absence of energy barriers with all other parameters as stated in Table 2. The average Sherwood number is computed by setting Rc ) 100 µm in eq 21, which is assumed to be the radius of the stagnation point region. The dependence of Sh on stagnation Peclet number (Pes) is depicted in Figure 5a. In this figure, symbols represent the numerical result, while the dashed line represents the predictions of Sh based on eq 25. The two results are in excellent agreement. At larger values of Pes > 1, the Sh obtained using the numerical method deviates from that predicted by Levich equation, since the latter is valid only for small Peclet numbers, while our numerical model captures the interception and gravitational deposition mechanisms at larger Pes.3 In particle deposition models, mechanisms of diffusion, interception, and gravity simultaneously affect the single-collector efficiency,3 although the influence of different mechanisms are manifested differently depending on the particle size. The singlecollector efficiency (η) is related to the collector average Sherwood number by22
η)
Sh Pe
(26)
Several correlations are available to relate the influence of different mechanisms on the single-collector efficiency.39-41 All correlations generally capture the well-known behavior that the variation of the single-collector efficiency with particle radius exhibits a minimum at a critical particle radius of ca. 500 nm. We calculate the particle deposition rate from the numerical model as Sh for different particle radii (ap) ranging from 0.01 to 2 µm. The scaled (38) Levich, V. G. Physicochemical Hydrodynamics; Prentice Hall: Englewood Cliffs, NJ, 1962. (39) Yao, K. M.; Habibian, M. T.; O’Melia, C. R. EnViron. Sci. Technol. 1971, 5, 1105. (40) Rajagopalan, R.; Tien, C. AIChE J. 1976, 22, 523. (41) Tufenkji, N.; Elimelech, M. EnViron. Sci. Technol. 2004, 38, 529.
Figure 5. (a) Variation of overall deposition rate with respect to stagnation point Peclet number for a homogeneous collector inside the stagnation point region. Dashed line denotes the values of Sh calculated by Levich equation (eq 25), whereas symbols represents the Sh obtained by solving convection-diffusion-migration equation (eq 23). (b) Variation of overall deposition efficiency, η, with respect to particle radius for a homogeneous fully favorable collector inside the stagnation point flow domain. Symbols denote the values of η calculated by solving convection-diffusionmigration equation numerically (Ψs ) +1, Ψp ) -1, Re ) 100, I ) 10-3 M). All other parameters are stated in Table 2. Dashed line denotes η obtained using Levich equation (eq 25).
surface potential of the homogeneous collector was Ψs ) +1, while that of the particle was Ψp ) -1. The flow Reynolds number and solution ionic strength were 100 and 10-3 M, respectively. All the other parameters are stated in Table 2. Figure 5b depicts the numerical results (symbols), where the singlecollector efficiency is plotted against the particle radius. Figure 5b shows the same qualitative behavior of η as reported in the literature.3,41-43 As expected, we observe a minimum in the curve corresponding to a particle radius of about 0.5 µm. The predictions of the single-collector efficiency based on the Levich equation (eq 25) is also plotted in the figure (dashed line). It is evident that for small particles ap < 0.5 µm, which correspond to smaller values of Pe < 1, the numerically calculated η is in excellent agreement with the η predicted by the Levich equation. As a third and final test of our numerical solution, we compare the finite element results against those provided by Dabros and van de Ven24 for deposition rates on a homogeneous unfavorable collector in stagnation point flow. Numerical solution of the convection-diffusion equation is considered difficult when repulsive EDL interactions are involved.8 In such situations, the particles accumulate in a narrow region located at a finite distance from the collector surface around the secondary minimum of the DLVO potential22 and form a sharp concentration peak that (42) Rajagopalan, R.; Kim, J. S. J. Colloid Interface Sci. 1981, 83, 428. (43) Nelson, K. E.; Ginn, T. R. Langmuir 2005, 21, 2173.
9886 Langmuir, Vol. 22, No. 24, 2006
Figure 6. Variation of the overall particle deposition rate with respect to flow Reynolds number for a homogeneous unfavorable collector in the presence of an energy barrier in the stagnation point region. Solid line denotes the values of Sh calculated numerically by Dabros and van de Ven for three different values of double-layer parameter, Dl. Dashed line represents the Sh obtained by solving convection-diffusion-migration equation numerically (eq 23).
changes the nature of convection-diffusion equation into stiff equations with possibilities of turning points.8 Figure 6 shows the variation of particle deposition rate, Sh, with respect to flow Reynolds number (Re) for different values of the double-layer parameter (Dl) in the stagnation point region. Solid lines represent the Sh reported by Dabros and van de Ven. Dashed lines represent the Sh obtained by finite element analysis in the present study. For each value of Dl, the finite element calculations were performed for different values of Re. All the parameters are identical to those used in the calculations of Dabros and van de Ven.24 Figure 6 depicts an excellent agreement between the two results, which implies the validity of the present numerical model. The comparisons demonstrate the ability of our finite element model to predict the deposition behavior over homogeneous collectors. We now turn our attention to the model predictions for heterogeneous collectors.
Results and Discussion In this section, the numerical results from the finite element analysis are used to demonstrate the deposition behavior on charge heterogeneous substrates. Radial and axial variations of particle concentration, as well as local and spatially averaged particle fluxes and deposition rates are obtained, which elucidate the influence of chemical heterogeneity on particle deposition. Particle Deposition on Homogeneous Collectors. To assess the influence of collector surface charge heterogeneity on the deposition process, it is first necessary to consider the deposition on homogeneously unfavorable and favorable surfaces. The estimates of the deposition rates from these calculations provide the benchmarks against which the deposition on heterogeneous collectors is compared. Figure 7 depicts the variation of dimensionless particle concentration with axial distance of the particle center from the collector surface (z) at the stagnation point (r ) 0). The particle radius is 1 µm, the flow Reynolds number is 100, and the solution ionic strength is 10-3 M. All other parameters are as stated in Table 2. The solid line in Figure 7 depicts the concentration profile with respect to z for a homogeneous, fully favorable collector where the particle and collector scaled surface potentials are Ψp ) -1 and Ψs ) +1, respectively. The dashed line represents the corresponding particle concentration profile for a homogeneous, fully unfavorable collector where the particle and collector scaled surface potentials are equal (Ψp ) Ψs )
Nazemifard et al.
Figure 7. Variation of the scaled particle concentration with the vertical distance from the collector surface, z, at the stagnation point, r ) 0. Solid line depicts the particle concentration profile over a homogeneous fully favorable collector with Ψs ) +1. Dashed line denotes the particle concentration profile over a homogeneous fully unfavorable collector with Ψs ) -1.
Figure 8. Variation of local particle deposition rate, Sh with the radial distance over the collector surface, r. Solid line depicts the particle deposition rate over a homogeneous fully favorable collector with Ψs ) +1. Dashed line denotes the particle deposition rate over a homogeneous fully unfavorable collector with Ψs ) -1.
-1). Far from the collector surface, the particle concentration approaches the bulk concentration for both favorable and unfavorable collectors. However, at a vertical distance of z ≈ 1.1 µm (s ≈ 100 nm) from the collector surface, the particles start to experience colloidal forces and the concentration profile over favorable and unfavorable collectors show marked differences. In the case of a favorable collector, the competition between the attractive colloidal forces and the repulsive hydrodynamic and gravitational interactions leads to a small and broad peak in the concentration profile. In contrast, the peak is more pronounced for an unfavorable collector, since the electrostatic repulsion additionally prevents the particles from approaching the collector. The particle concentration is zero at the collector surface for both favorable and unfavorable collectors owing to the perfect sink boundary condition. For an unfavorable collector, particles accumulate at a distance of about z ) 1.072 µm (s ) 72 nm) from the collector surface, forming a sharp peak in the concentration profile. The separation distance s ) 72 nm coincides with the location where the total force exerted on the particle in the normal direction (including gravitational, hydrodynamic, and colloidal forces) becomes zero. Figure 8 shows variation of the local Sherwood number, Sh(r), with radial distance from the stagnation point, r, for homogeneously favorable and unfavorable collectors. All conditions are identical to those of Figure 7. On a favorable collector, Sh(r) (solid line) is usually constant and independent of r up to the radial distance of ca. 150 µm from the stagnation point, implying that the stagnation point flow system acts as a uniform collector.
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9887
Figure 9. (a) Variation of the scaled particle concentration with the radial distance from the stagnation point, r, on a heterogeneous collector at a vertical distance of s ) 5 nm from the collector surface. The regions inside the stagnation point (r < 25 µm), the region spanning 540 < r < 580 µm, and the region spanning 940 < r < 980 µm are enlarged and shown in panels b-d, respectively. In each of these panels, the location of the favorable stripes are indicated by the gray horizontal rectangles. The flow Reynolds number and the solution ionic strength are 100 and 10-3 M, respectively. All the other parameters are stated in Table 2.
At larger radial distances, Sh(r) decreases significantly due to the flow distribution in the impinging jet system, which is characterized by the increase in the ratio of the radial to normal velocity components at larger r.44 In other words, the collector becomes nonuniformly accessible in the presence of a tangential flow at large radial distances from the stagnation point. For a homogeneously unfavorable collector, on the other hand, due to the repulsive electrostatic double-layer interactions between the particles and the collector, Sh(r) is zero over the entire collector surface. Particle Deposition on Heterogeneous Collectors. Figures 9 and 10 show the radial dependence of the scaled particle concentration corresponding to two different separation distances (s) from the collector surface. These results are obtained for a collector consisting of alternate negative and positive stripes with scaled surface potentials of Ψs,n ) -1 and Ψs,p ) +1, respectively. The widths of negative and positive stripes are both 10 µm giving λpatterned ) 50%. The particle radius and scaled surface potential are ap ) 1 µm and Ψp ) -1, respectively. The solution ionic strength and Reynolds number are 10-3 M and 100, respectively. Figure 9a depicts the particle concentration close to the collector surface at s ) 5 nm over the entire collector surface (Rc ) 1000 µm). Figure 9, panels b-d, shows enlarged views of the concentration profiles at the stagnation point (r < 25 µm), the region spanning 540 < r < 580 µm, and the region spanning 940 < r < 980 µm, respectively. In each panel, the location of the favorable (positively charged) stripes are indicated by the gray horizontal rectangles. The particle concentration profile varies periodically with r. The concentration is zero over the unfavorable stripes due to the repulsive colloidal interactions. (44) Adamczyk, Z.; Siwek, B.; Warszynski, P.; Musial, E. J. Colloid Interface Sci. 2001, 242, 14.
Figure 10. Variation of the scaled particle concentration with the radial distance from the stagnation point, r, on a heterogeneous collector at a vertical distance of s ) 72 nm from the collector surface. The region spanning 540 < r < 580 µm is shown in the inset. All other parameters are the same as in Figure 9.
In contrast, on each favorable stripe, the concentration is finite, and there is a sharp increase in the particle concentration at the beginning of each favorable stripe. The concentration profile undergoes substantial change with increasing radial distance from the stagnation point. In the stagnation point region (Figure 9b), a sharp concentration front is developed at the leading edge of each favorable stripe. Following this, the concentration decays rapidly and attains a constant (finite) value over the remainder of the favorable stripe. The concentration profile closely emulates the periodic nature of the heterogeneity in this case. At intermediate radial distances (Figure 9c), the concentration profiles do not exhibit the sharp peaks at the leading edge of the favorable stripes any more.
9888 Langmuir, Vol. 22, No. 24, 2006
Figure 11. Schematic representation of the radial particle flux from the unfavorable stripes toward the favorable stripes due to the radial convection and diffusion mechanisms.
Instead, the highest concentration is attained at intermediate locations on the positive stripes. Curiously, a small region near the leading edge of each favorable stripe has virtually the same concentration as over the previous unfavorable stripe. This behavior is in conformity with the presence of an inaccessible region at the leading edge of each favorable stripe as discussed in our earlier study,20 where a trajectory analysis was employed to arrive at the same conclusion. Furthermore, at the trailing edge of each favorable stripe, we observe a second peak in concentration. The trends in the concentration profile observed in this case is replicated at even larger radial positions (Figure 9d), although the trends become smoother with increasing radial distance from the stagnation point. Figure 10 depicts the particle concentration distribution with radial position on the collector corresponding to s ) 72 nm. The region spanning 540 < r < 580 µm is shown in the inset. At this separation distance, the net force acting on the particle normal to the collector surface vanishes. However, the increased radial velocity at this separation smooths out the sharp features of the concentration distribution observed in Figure 9. The concentration increases radially over the unfavorable stripe and decreases over the favorable stripe. The peak value of the concentration decreases with increasing radial distance from the stagnation point. The presence of a concentration peak at the beginning of each favorable stripe can be explained by taking into account the coupled effects of hydrodynamic and colloidal interactions. Figure 11 is a schematic representation of these effects on the radial particle flux at the boundaries between the unfavorable and favorable stripes. At the distance of s ) 72 nm from the collector surface, particles accumulate over the unfavorable stripes since they cannot get closer to the collector due to the presence of the energy barrier (Figure 11a). These particles can, however, be transferred radially toward the next favorable stripe by radial convection and diffusion (Figure 11b). Consequently, the particles that could not deposit on the unfavorable stripes can be swept toward the next accessible favorable stripe, causing an enhanced concentration at its leading edge. In a similar manner, the sharp concentration peaks observed at the trailing edge of the favorable stripes (Figure 9c) can be explained by considering the migration flux component normal to the collector. As the particle approaches the end of a favorable stripe, it experiences less attraction by the collector. As a result, the particle experiences a lower migration flux toward the collector (under some conditions, the direction of the normal migration flux may even reverse). Coupled with radial convection, this causes an accumulation of particles at the
Nazemifard et al.
trailing edge of a favorable stripe at a given separation (s ) 5 nm in this case), since particles closer to the collector (s < 5 nm) are transported toward this plane. In summary, a coupled influence of the two-dimensional velocity field and the charge heterogeneity can result in a highly complex concentration distribution near the collector. For large particle sizes (ap g 1 µm), the convective flux is the dominant mechanism of radial particle transport, while for small particle sizes (ap < 1 µm), both diffusive and convective fluxes play a comparable role in particle transport. The variations of the radial diffusive and convective fluxes with respect to r at the distance of s ) 72 nm are shown in Figure 12. Figure 12, panels a and b, depicts the variations of radial convective and diffusive fluxes for ap ) 1 µm whereas Figure 12, panels c and d, depicts the variation of radial convective and diffusive fluxes for ap ) 0.1 µm. All other parameters are the same as those of Figure 9. Both radial convective and diffusive fluxes increase at the end of each unfavorable stripe. The only difference is that for the 1 µm particle, the radial convective flux is significantly higher than radial diffusive flux while for the 0.1 µm particle, the radial diffusive and convective fluxes are comparable in magnitude. Local Particle Deposition Rate. We now focus on particle deposition rate over heterogeneous collectors with different values of favorable area fraction, λpatterned. Comparing the particle deposition rate over the heterogeneous collectors with those obtained for homogeneous collectors reveals the influence of charge heterogeneity on the particle deposition behavior. Using the particle concentration distributions in eq 18, one can calculate the local particle deposition rate (Sh) over a heterogeneous collector. Figure 13 depicts the radial variation of the local Sherwood number for a heterogeneous collector corresponding to the concentration distribution depicted in Figure 9. It is evident from Figure 13a that the local Sherwood number shows a periodic variation with r similar to what was observed in Figure 9a for particle concentration profile close to the collector surface. The regions of 0 < r < 100 µm and 500 < r < 600 µm in Figure 13a are enlarged and depicted in Figure 13, panels b and c, respectively. It can be noted from Figure 13 that the particle deposition rate is zero over the unfavorable stripes throughout the collector surface. Near the stagnation point, the particle deposition rate increases drastically at the leading edge of each favorable stripe (Figure 13b). In contrast, at larger radial distances from the stagnation point, the deposition rate at the leading edge of each favorable stripe remains zero, while the peak shifts to slightly downstream locations (Figure 13c). A comparison of Figures 9 and 13 clearly shows that the particle deposition rate is zero where the particle concentration is zero, and the particle deposition rate reaches its maximum exactly at the same radial position where the particle concentration reaches its maximum. A significant difference in the deposition rates over homogeneous and heterogeneous collectors becomes immediately evident by comparing the local Sherwood numbers shown in Figures 8 and 13. Comparing the two figures, it is noted that the maximum deposition rate over a homogeneously favorable collector at the stagnation point is about 9 (Figure 8), while the peak value of the deposition rate over the first favorable stripe in Figure 13 is about 110. The enhanced local deposition on the favorable stripe of the heterogeneous collector is caused by the radial transport of particles from the previous unfavorable stripe. It is, however, worth noting that only 50% of the heterogeneous collector encounters particle deposition. Accordingly, although the favorable stripes have a higher local Sh, the Sherwood number averaged over the entire collector surface might be comparable to that for the homogeneously favorable collector.
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9889
Figure 12. Variation of scaled (a) radial convective flux for particle with radius, ap ) 1 µm; (b) radial diffusive flux for particle with radius, ap ) 1 µm; (c) radial convective flux for particle with radius, ap ) 0.1 µm; and (d) radial diffusive flux for particle with radius, ap ) 0.1 µm with the radial distance from the stagnation point, r, at the vertical distance of s ) 72 nm from the collector surface. All the other parameters are the same as in Figure 9. In all these panels, the gray horizontal rectangles represent the favorable stripes while the blank horizontal rectangles represent the unfavorable strips.
Average Particle Deposition Rate. In most deposition applications, one is generally interested in the overall particle deposition rate on a substrate. The overall particle deposition rate is represented by Sh, which can be calculated using eq 19. As described earlier, the average deposition rates can be defined in different ways. First, an average particle deposition rate over the entire collector surface, 0 e r e Rc can be calculated using eq 21. Second, an average Sherwood number over each heterogeneous band (Shp) can be calculated using eq 22. This second averaging procedure provides a “local” average Sherwood number on a heterogeneous substrate, which varies with the band number (Nb) or the radial position of a heterogeneous band (rNb) on the collector. Figure 14 shows a comparison between the Shp obtained for the heterogeneous collector with the corresponding local Sherwood number at different radial positions on a homogeneously favorable collector. The Sherwood numbers are plotted against the band number (Nb), which can be easily converted to radial position by noting that rNb ) (Nb - 1)p, with p ) 20 µm. The results are obtained for four different particle sizes: (a) ap ) 1 µm, (b) ap ) 0.5 µm, (c) ap ) 0.25 µm, and (d) ap ) 0.1 µm. The heterogeneous collector is modeled using parameters as outlined in Figure 9. The scaled surface potential of the homogeneous collector is Ψs ) +1. All other parameters are as stated in Table 2.
It is evident from Figure 14 that the local average Sherwood number (Shp) on a heterogeneous collector exhibits a fairly close resemblance to the radial variation of the local Sherwood number (Sh) on a homogeneously favorable collector. The difference between the two increases as the particle size diminishes. The deposition rates on the heterogeneous surface containing 50% unfavorable sites are quite similar to the corresponding deposition rates on a homogeneously favorable surface. This warrants a closer inspection of a common approach employed for assessing chemical heterogeneity effects on particle deposition rates, namely, the patchwise heterogeneity model.8,9 The predictions of the patchwise heterogeneity model are also shown in Figure 14 (dashed lines), and it is clear that the model does not agree with our numerical predictions of the deposition rate on a heterogeneous collector. In the following, we investigate the reasons for this discrepancy. Patchwise Charge Heterogeneity Model. According to the patchwise heterogeneity model, the overall particle deposition rate (Sh) over a heterogeneous surface consisting of unfavorable and favorable patches can be given in context of the Eulerian approach as8,9
Sh ) λpatterned Shf + (1 - λpatterned) Shu ) λpatterned Shf (27) where Shf and Shu are the average particle deposition rates over the homogeneously favorable and unfavorable surfaces, respec-
9890 Langmuir, Vol. 22, No. 24, 2006
Nazemifard et al.
Figure 13. Variation of local particle deposition rate, Sh, with radial distance from the stagnation point, r, at the collector surface. The region spanning 0 < r < 100 µm and the region spanning 500 < r < 600 µm are enlarged in panels b and c, respectively. All other parameters are the same as in Figure 9.
tively. The final expression is obtained noting that in our theoretical calculations, Shu ) 0. Equation 27 gives a linear relation between the average particle deposition rate on a heterogeneous collector (Sh) and the favorable collector area fraction (λpatterned). This behavior is shown in Figure 14, where the dashed lines represent the patchwise heterogeneity model predictions for λpatterned ) 0.5, with Sh ) 0.5Shf. The patch model assumes that the particle deposition rates on the favorable and unfavorable regions of the collector are unaffected by the heterogeneity. In other words, the favorable and unfavorable sites are assumed to have deposition rates that are equal to the corresponding deposition rates on perfectly homogeneous analogues of these sites with identical properties. However, a comparison between values of local Sh in Figures 8 and 13 shows that at the same radial distance from the stagnation point, the local deposition rate over each favorable stripe on the heterogeneous collector is significantly larger than the corresponding local deposition rate over the homogeneous collector. Stated differently, the average Sherwood number on a favorable stripe of a heterogeneous collector, which may be denoted by Shf,p, can differ substantially from the average Sherwood number on its homogeneous analogue, Shf owing to a coupling between the various transport mechanisms and the charge
heterogeneity. In contrast, the patchwise heterogeneity model assumes that Shf,p ) Shf. Thus, when Shf,p/Shf > 1, the average deposition rate Sh will be greater than the corresponding prediction of the patchwise heterogeneity model. It may also be possible that Shf, p/Shf < 1, in which case Sh for the heterogeneous collector will be less than the corresponding predictions of the patchwise heterogeneity model. The relatively higher values of the local deposition rate over the favorable stripes compensate for the zero deposition rate over the unfavorable stripes during the radial averaging process for calculating the deposition rates on the heterogeneous collector in Figure 14. To summarize, modification of the local deposition rates arising from the coupling of charge heterogeneity with the particle transport mechanisms is not accounted for in the patchwise heterogeneity model. As a consequence, in the presence of microscale charge heterogeneity, the model breaks down. Influence of Favorable Surface Area Fraction. So far, for all of the cases that have been studied, the width of favorable and unfavorable stripes were equal, implying that half of the heterogeneous collector was covered with favorably charged stripes. We now discuss the dependence of the deposition rate on the favorable area fraction (λpatterned) of the collector. For a heterogeneous collector with a constant pitch, different values
ConVection-Diffusion-Migration Model
Langmuir, Vol. 22, No. 24, 2006 9891
Figure 14. Variation of Shp with band number, Nb, for four particle sizes, i.e.: (a) ap ) 1 µm, (b) ap ) 0.5 µm, (c) ap ) 0.25 µm, (d) ap ) 0.1 µm. Solid line with triangular symbols denotes Shp obtained by numerical solution of convection-diffusion-migration equation over a heterogeneous collector. Dashed line denotes Shp predicted by the patchwise heterogeneity model. The heterogeneous collector consists of alternate negative and positive stripes with scaled surface potentials of -1 and +1, respectively. The collector favorable fraction is λpatterned ) 50%. Solid line with circular symbols denotes Shp over a homogeneous fully favorable collector with ψs ) +1. All other parameters are the same as in Figure 9.
of λpatterned can be achieved by changing the ratio of the width of the unfavorable to the favorable stripe (wn/wp). All other parameters used to obtain the following results are identical to those of Figure 9 and Table 2. Figure 15a shows the variation of overall deposition rate (Sh) scaled with respect to the deposition rate over a homogeneously favorable collector (Shf) with the favorable collector area fraction (λpatterned). A completely unfavorable collector is represented by λpatterned ) 0, while a completely favorable collector is represented by λpatterned ) 1. The square symbols denote the Sh/Shf obtained numerically in the present work. The dashed line denotes the Sh/Shf predicted by the patchwise heterogeneity model (eq 27). It is immediately evident that the overall deposition rate on the patterned collector deviates significantly from that predicted by the patchwise heterogeneity model. The numerical simulations in Figure 15a (square symbols) show that Sh is zero for a completely unfavorable collector (λpatterned ) 0). For λpatterned < 0.5, a slight increase in the λpatterned will result in a significant increase in Sh. This has been reported earlier in the literature.45 (45) Adamczyk, Z.; Siwek, B.; Weronski, P.; Jaszczolt, K. Colloids Surf. A 2003, 222, 15.
Furthermore, it can be observed that for λpatterned > 0.5, Sh over the heterogeneous collector becomes quite independent of the favorable area fraction, and closely resembles the deposition behavior on a completely favorable collector. At λpatterned ) 0.5, the average Sherwood number over the heterogeneous collector is within 10% of the corresponding Shf. This behavior was also evident in Figure 14 where the overall deposition rate over a heterogeneous collector with λpatterned ) 0.5 was very close to Shf over a homogeneously favorable collector. Comparing the overall deposition rates obtained by finite element analysis with those predicted by the patchwise heterogeneity model, it is evident that the overall deposition rate (Sh) does not vary linearly with λpatterned as predicted by the patchwise heterogeneity model. The nonlinear relation between overall deposition rate and favorable area fraction of the collector was also observed earlier20 employing trajectory analysis to calculate the single-collector efficiency, η () Sh/Pe). The resulting η versus λpatterned plots are shown in Figure 15b. These results correspond to a flow Reynolds number of 180 (U∞ ) 0.18 m/s), solution ionic strength of 10-3 M, particle radius of 1 µm, and particle surface potential, Ψp ) -2. The scaled surface potentials of the heterogeneous collector are Ψs,n ) -2 and Ψs,p ) +2.
9892 Langmuir, Vol. 22, No. 24, 2006
Nazemifard et al.
the actual favorable area fraction of the collector surface, one may conclude that the overall particle deposition rate should also decrease. However, the limiting trajectory analysis predicted particle deposition rates that were significantly higher than those predicted by the patchwise heterogeneity model. The convectiondiffusion model, on the other hand, provides a comprehensive picture of the particle concentration distribution and deposition rates over heterogeneous collectors. In the present study, we observe that although there is a narrow region at the leading edge of each favorable stripe, which is inaccessible for particle deposition, the particle deposition rates over the favorable stripes are significantly larger than the corresponding deposition rates on the homogeneously favorable collectors. The enhanced deposition rates are caused by the radial particle flux from the unfavorable stripes. The coupling between charge heterogeneity and particle transport mechanisms increases the overall particle deposition rate over the favorable fractions of the heterogeneous collector. This implies that the particle deposition rates obtained earlier20 using the limiting trajectory analysis were correct, and both Lagrangian and Eulerian approaches provide unambiguous predictions of the particle deposition rate over heterogeneous surfaces.
Concluding Remarks
Figure 15. (a) Variation of scaled overall particle deposition rate, Sh/Shf, with favorable area fraction of the collector, λpatterned. The square symbols denote the Sh/Shf obtained numerically in the present work. The dashed line denotes the Sh/Shf predicted by the patchwise heterogeneity model (eq 27). All other parameters are the same as in Figure 9. (b) The dependence of the scaled collector deposition efficiency (η/ηf) on the collector favorable area fraction, λpatterned. The square symbols denote the deposition efficiency for a heterogeneous collector obtained using trajectory model. The dashed line denotes the scaled collector deposition efficiency predicted by the patchwise heterogeneity model. The solution ionic strength is I ) 10-3 M. The scaled particle and collector surface potentials are Ψp ) -2 and Ψs ) (-2, +2), respectively.
The pitch is constant and equal to 10 µm. The numerically obtained particle deposition efficiencies (solid line with square symbols) are calculated by first finding the limiting trajectory corresponding to Rc ) 300 µm for different values of λpatterned. Detailed description of the method for estimating the single-collector efficiency (η) from the limiting trajectories is available elsewhere.20,23 Figure 15b depicts the variation of scaled single-collector efficiency (η/ηf) with respect to λpatterned for a heterogeneous collector with radius of Rc ) 300 µm. The solid line with square symbols represents the particle deposition rate obtained numerically using the trajectory analysis.20 The dashed line represents the particle deposition rate predicted by the patchwise heterogeneity model. Figure 15, panels a and b, shows an identical qualitative dependence of the overall deposition rate on λpatterned. It is worth noting here that although both the convectiondiffusion-migration model and the trajectory analysis predict similar particle deposition rates over charge heterogeneous substrates, the convection-diffusion model provides a better understanding of the particle distributions in the vicinity of the charge heterogeneous substrates. In our earlier study,20 it was observed that a small region at the leading edge of each favorable stripe was inaccessible for particle deposition. Since this decreases
A comprehensive analysis of the particle concentration distributions and corresponding particle deposition rates on a patterned charge heterogeneous collector is presented based on the Eulerian approach. Using this model, it was observed that there is a radial flux at the boundary between unfavorable and favorable stripes that transports the particles accumulated over the unfavorable stripes toward the next accessible favorable region. This causes a sharp increase in the particle concentration, as well as the particle deposition rate on the accessible regions of the favorable stripes. The particle deposition rates on the favorable stripes of a heterogeneous collector are considerably larger than the corresponding deposition rates on their perfectly homogeneous analogues. The large particle deposition rates on the favorable stripes lead to higher overall particle deposition rates on heterogeneous collectors. In contrast, the patchwise heterogeneity model assumes that the particle deposition rate on the favorable patches remains unaffected by charge heterogeneity. This can lead to erroneous predictions of the influence of microscopic charge heterogeneity on overall particle deposition rates. The linear dependence of the overall particle deposition rate on favorable collector area fraction predicted by the patchwise heterogeneity model was not observed in the detailed simulations. The variation of overall particle deposition rate with favorable area fraction obtained using both Eulerian and Lagrangian approaches seem to be in good qualitative conformity. Nomenclature A ) Hamaker constant [J] Ad ) dimensionless adhesion number ap ) particle radius [m] c ) particle concentration [mol m-3] c∞ ) bulk particle concentration [mol m-3] cj ) scaled particle concentration (c/c∞) D∞ ) Stokes-Einstein particle diffusion coefficient [m2 s-1] D ) particle diffusion tensor [m2 s-1] Da ) dimensionless electrostatic double-layer asymmetry parameter Dl ) dimensionless electrostatic double-layer parameter e ) electronic charge [C] B F ) external force acting on a particle [N] Fr ) sum of the forces in r-direction [N] Fz ) sum of the forces in z-direction [N]
ConVection-Diffusion-Migration Model FBr ) Brownian force [N] Fedl ) electrostatic double-layer force [N] Fg ) gravity force [N] Fvdw ) van der Waals force [N] fi ) particle hydrodynamic correction functions (i ) 1 ... 4) b g ) gravitational acceleration [ms-2] h ) dimensionless surface to surface separation distance (s/ap) I ) solution ionic strength [M] bj ) particle flux vector [mol m-2 s-1] jj ) scaled particle flux (jap/(D∞c∞)) kb ) Boltzmann constant [J K-1] L ) vertical distance of impinging jet from collector surface [m] Nb ) band number in radial direction p ) pitch of the charge heterogeneous band () wn + wp) [m] Pe ) Peclet number (U∞ap/D∞) Pes ) stagnation flow Peclet number (2Rsap3/D∞) r ) radial distance [m] jr ) dimensionless radial distance (r/ap) rNb ) radial distance of the leading edge of the band Nb [m] Rc ) collector radius [m] Re ) Reynolds number (U∞Rjet/ν) Rjet ) radius of the jet tip [m] Sh ) local Sherwood number (same as jj) Sh ) average Sherwood number over a collector Shp ) averaged Sherwood number over a band (comprising an unfavorable and a favorable stripe) Shf ) average Sherwood number over a homogeneously favorable collector Shf,p ) average Sherwood number over a favorable stripe in a heterogeneous collector Shu ) average Sherwood number over a homogeneously unfavorable collector
Langmuir, Vol. 22, No. 24, 2006 9893 s ) surface to surface separation distance between the particle and the collector surface [m] S ) total surface area of the collector () πRc2) [m2] Sf ) surface area of the favorable region on the collector [m2] U∞ ) average fluid velocity in jet [ms-1] b u ) fluid velocity vector [ms-1] b V ) particle velocity vector [ms-1] wp ) width of positive stripe [µm] wn ) width of negative stripe [µm] z ) normal distance of particle center from collector surface [m] jz ) dimensionless normal distance (z/ap) Greek Symbols Rs ) impinging jet flow coefficient [1/ms] R j ) dimensionless impinging jet flow coefficient δm ) cutoff separation distance ) dielectric constant 0 ) permittivity of vacuum [C/Vm] η ) single-collector efficiency ηf ) single-collector efficiency of a fully favorable collector λpatterned ) favorable area fraction of the collector κ ) inverse Debye length [m-1] Ψp ) scaled particle surface potential Ψs ) scaled collector surface potential Ψs,p ) scaled surface potential of positive stripes Ψs,n ) scaled surface potential of negative stripes ∂Ω ) boundary of computational domain µf ) dynamic viscosity of the fluid [Ns/m2] ν ) kinematic viscosity of the fluid [m2‚s-1] LA061702Q
