Numerical simulation of electrospray in the cone-jet mode

Numerical simulation of electrospray in the cone-jet mode M. A. Herrada, J. M. L´opez-Herrera, and A. M. Ga˜n´an-Calvo Departamento de Ingenier´ıa Aeroespacial y Mec´anica de Fluidos, Universidad de Sevilla, E-41092 Sevilla, Spain E. J. Vega and J. M. Montanero Departamento de Ingenier´ıa Mec´anica, Energ´etica y de los Materiales, Universidad de Extremadura, E-06006 Badajoz, Spain S. Popinet National Institute of Water and Atmospheric Research, P.O. Box 14-901, Kilbirnie, Wellington, New Zealand (Dated: July 6, 2012) Abstract We present a robust and computationally efficient numerical scheme for simulating steady electrohydrodynamic atomization processes (electrospray). The main simplification assumed in this scheme is that all the free electrical charges are distributed over the interface. Comparison of the results with those calculated with a Volume-of-Fluid method showed that the numerical scheme presented here accurately describes the flow pattern within the entire liquid domain. Experiments were performed to partially validate the numerical predictions. The simulations reproduced accurately the experimental shape of the liquid cone-jet, providing correct values of the emitted electric current even for configurations very close to the cone-jet stability limit. PACS numbers: 47.61.Jd, 47.65.-d, 47.55.-t, 47.15.-x Keywords: Electrohydrodynamic, EHD, Electrospray, VOF 1 I. INTRODUCTION The use of electrostatic fields to spray liquids is a procedure that has been used since the XVIIth century [1]. In the most typical configuration, the liquid is injected at a small flow rate through a metallic needle where an electric potential of the order of kilovolts is imposed with respect to a grounded electrode situated downstream. Although many electrospray modes have been identified [2], the one which has attracted most attention is the cone-jet mode. In this mode, the meniscus adopts a stable quasi-conical shape from whose apex a stationary thin jet is ejected. The jet eventually breaks up into droplets due to the capillary forces, which gives rise to a charged mist. All theories describing the physics behind the cone-jet mode distinguish three regions [4, 5]: (i) the upstream meniscus region attached to the needle, where ohmic conduction is the dominant mechanism of charge transport and the voltage drop is negligible; (ii) the downstream jet region, where the relevant charge transport mechanism is surface convection; and (iii) an intermediate transition region (the cone neck), where the charge transport mechanism evolves from conduction-dominant to convection-dominant. It is generally accepted that both the electric current circulating through the system and the jet’s diameter are essentially governed by this neck region. Two theoretical approaches have been proposed to model this critical region: one considers electrical relaxation phenomena in this region [6], and the other assumes that the electric charges are relaxed over the entire interface, including the cone neck region [4, 7, 8]. They propose different scaling laws for both the jet’s size and issued electric current. Experimental data for the size of the resulting droplets seem to indicate that the jet’s diameter scales with the governing parameters as indicated by the later theory [9]. However, there have as yet been no experimental data on the neck clearly in favour of one or the other of these two approaches. There have been several numerical studies of electrospray in the literature. The first numerical simulation was due to Hartman et al. [10] who used a one-dimensional model to reach steady cone-jet configurations by successive iterations. Their study focused on describing the meniscus and the neck regions. Yan et al. [11] refined that model by solving the axisymmetric Navier-Stokes equations and incorporating grid adaptation into their calculations. An additional constraint imposed in these two studies, even though it is not common operational practice, is that the jet is forced to hit a flat grounded electrode located down- 2 stream. Higuera [12] examined numerically the neck region in order to analyze the electric current circulating through the system. The presence of the needle and the breakup region were obviated by using as up- and down-stream boundary conditions far-field asymptotic behaviors for the conical meniscus [13, 14] and the jet [4], respectively. In all these studies, the flow was assumed to be strictly steady, thus neglecting the growth of perturbations on the free surface and the subsequent jet breakup. Collins et al. [15] proposed a more sophisticated numerical scheme to study the entire electrohydrodynamic (EHD) tip streaming flow, including the formation of the cusp, the emitted jet, and the resulting droplets. Despite the disparity of the length scales present in this phenomenon, it was accurately described by this numerical approach. Reznik et al. [16] studied the temporal evolution of an electrified compound core-shell droplet towards the emission of a coaxial jet. All the aforementioned numerical studies have in common that, in the bulk, the electrical relaxation time was short as compared to any hydrodynamic time present in the process, and thus the liquid bulk was free of electrical charges, which were only allowed to exist at the interface. Consequently, the surface charge density σewas regarded as an unknown governed by the surface balance equation [6] d(2πFusσe) ds= 2πFKEni ,(1) where sis the distance along the interface, Fand usare the interface radial position and velocity, respectively, Kis the liquid electrical conductivity, and Eni is the normal component of the electric field on the inner side of the interface. Besides the aforementioned works, it can be found in the literature numerical analyses of the cone jet mode in which unconventional models for the electrical behaviors have been used. Lastow and Balachandran [17] assumed a dielectric behavior by neglecting the electric current, what seems to be a clear oversimplification since both experimental and theoretical works point to the relevance of the conductivity in the scaling of the emitted droplets. In the same line, Lim et al. [18] used a constant surface charge density over the entire cone-jet structure, adjusting its value by fitting numerical simulations to experiments. Both works showed a strong stretching of the meniscus by electrical means despite the feeble physical ground on which their electric models were built. In the present work, a robust and computationally efficient numerical method is developed to simulate steady EHD atomization processes. As in previous approaches [10–12, 15], the electric charges are confined to the interface. The results are compared with those 3 calculated with Gerris, an open source free software library developed by Popinet [19, 20] and extended by L´opez-Herrera et al. [21] to study EHD phenomena. Experiments are conducted to partially validate the numerical results. The cone-jet transition is examined very close to the stability limit to assess the validity of the hypothesis on the charge behavior. The rest of the paper is organized as follows. The problem is formulated in Sec. II. Sections III and IV give the descriptions of the numerical schemes and the experimental method, respectively. Finally, the results are presented and discussed in Sec. V. II. FORMULATION OF THE PROBLEM The geometrical and electrical configurations considered in this work are represented schematically in Fig. 1. A metallic needle of inner radius Riand outer radius Rois held at a constant voltage V. The needle is brought face to face up close to a planar grounded electrode located at a distance H′. A liquid of density ρ, viscosity µ, electrical conductivity K, and permittivity εis injected through the needle at a constant flow rate Q. The flow is fully developed inside the capillary, so that there is a parabolic Hagen-Poiseuille velocity profile upstream at a distance Lnfrom the needle’s exit. The ambient medium is a perfect dielectric gas of permittivity εoequal to that of the vacuum. The gas-liquid surface tension is γ. Because of the small values of the Bond number and gas/liquid density and viscosity ratios, the effects of gravity and the gas dynamics can be neglected. The equations are made dimensionless with the inner radius Ri, the liquid density ρ, the surface tension γ, and the applied voltage V, which yields the characteristic time and velocity scales tc= (ρR3 i/γ)1/2and vc=Ri/tc, respectively. The dimensionless axisymmetric incompressible Navier-Stokes equations for the liquid velocity v(r, z;t) and pressure p(r, z;t) fields are (ru)r+rwz= 0,(2) ut+uur+wuz=−pr+C[urr + (u/r)r+uzz] + χfr e,(3) wt+uwr+wwz=−pz+C[wrr +wr/r +wzz] + χfz e,(4) where tis the time, r/zis the radial/axial coordinate, u/wis the radial/axial velocity component, χ=εoV2/(Riγ) is the electric Bond number, and C=µ(ργRi)−1/2is the Ohnesorge number. The subscripts t,r, and zhere and henceforth denote the partial 4 r z r' R o z' L nH H' R i R e = V 2 ( r ) 1 ( R e , z ) / z= 1 / z / r =0 r = F ( z ) FIG. 1. (Color online) Sketch of the problem’s formulation. The rectangle (red online) denotes the limits of the computational domain. derivatives with respect to the corresponding variables. The terms fr e/fz estand for the radial/axial volumetric electric force,   fr e fz e =ρe  Er Ez −(Er+Ez)2 2  βr βz ,(5) where ρeis the volumetric charge density, Er/Ezis the radial/axial electric field, and β= ε/εois the relative permittivity. In addition to the Navier-Stokes equations for the liquid flow, the electric potentials ϕi and ϕoin both the inner (liquid) and outer (gas) domains must be calculated. They obey the Poisson and the Laplace equations ϕi zz +ϕi rr +ϕi r/r =ρe/β, (6) ϕo zz +ϕo rr +ϕo r/r = 0,(7) respectively. Taking into account the kinematic compatibility and equilibrium of tangential and normal 5 stresses at the interface r=F(z), one gets the following equations: Fzw−u= 0,(8) C(1 −F2 z)(wr+uz)+2Fz(ur−wz) (1 + F2 z)1/2=χσeEt,(9) p+FFzz −1−F2 z F(1 + F2 z)3/2−2C[ur−Fz(wr+uz) + F2 zwz] 1 + F2 z =χ 2[(Eno)2−β(Eni)2+ (β−1)(Et)2], (10) where σeis the surface charge density, Etis the tangential component of the electric field, and Eni and Eno are the normal components of the electrical field on the inner and outer sides of the interface, respectively. The right-hand sides of Eqs. (9) and (10) are the Maxwell stresses resulting from the accumulation of free electric charges at the interface and the jump of permittivity across this surface. The electrical field and the surface charge density can be calculated from the inner ϕi and outer ϕoelectric potentials, and the free surface shape F: Eni =Fzϕi z−ϕi r √1 + F2 z ,(11) Eno =Fzϕo z−ϕo r √1 + F2 z ,(12) Et=−Fzϕo r+ϕo z √1 + F2 z =−Fzϕi r+ϕi z √1 + F2 z ,(13) σe=Eno −βEni,(14) Note that the continuity of the electric potential across the interface, ϕi=ϕo, has been taken into account in Eq. (13). The interfacial equations are completed by imposing the surface charge conservation at the free surface r=F(z) [Eq. (1)] [22]: Fzu+w 1 + F2 z σez−σe 1 + F2 z[ur+F2 zwz−Fz(uz+wr)]=αEni,(15) where α=β[ρR3 iK2/(γε2)]1/2is the dimensionless electrical conductivity. We shall assume that all characteristic times of hydrodynamic nature are long as compared to the electrical relaxation time (Saville’s leaky dielectric model Saville [23]). Therefore, the volumetric charge density ρeis zero, and the liquid conductivity Kis constant over the entire liquid domain. Also, the liquid electrical permittivity is assumed to be homogeneous, and thus the volumetric dielectrophoretic force proportional to the permittivity 6 gradient is zero. In this case, the electrical effects on the liquid dynamics are exclusively associated with the Maxwell stresses at the interface. As mentioned above, a Hagen-Poiseuille velocity profile is assumed at the entrance of the liquid domain z= 0: u= 0, w = 2ve(1 −r2),(16) where ve=Q/(πR2 ivc) is the dimensionless mean velocity. At the needle wall, the electric potential is fixed and the liquid does not slip, i.e., ϕi=ϕe= 1 and u=w= 0.(17) We impose the standard regularity conditions ϕi r=u=wr= 0 (18) on the symmetry axis, and the outflow conditions uz=wz=Fz= (σe)z= 0 (19) at the right-hand end ze= (H+Ln)/Riof the computational domain (the red rectangle in Fig. 1). For the configuration considered in this work, Ga˜n´an-Calvo et al. [24] found the following analytical solution for the far-field electric potential: ϕ1(r′, z′) = −Kv log(4H′/Ri)log {[r′2+ (1 −z′)2]1/2+ (1 −z′) [r′2+ (1 + z′)2]1/2+ (1 + z′)},(20) where r′and z′are cylindrical coordinates with origin at the intersection between the symmetry axis and the grounded planar electrode (see Fig. 1), and Kvis a dimensionless constant which depends on the ratio H′/Ri. We take advantage of this solution to limit the calculations to the meniscus and the jet, leaving aside the charged spray. Specifically, we assume that the electric potential matches the analytical solution ϕ1at the boundary r=re (re≡Re/Ri). A logarithmic drop of voltage ϕ2= 1 −[1 −ϕ1(re, z′ e)] log r/ log re, z′ e≡(H′+Ln)/Ri,(21) is applied at the boundary z= 0 and 1 < r < re. Finally, the condition ϕz= (ϕ1)z(22) is considered for the right-hand end z=zeof both the liquid and gas computational domains. 7 III. NUMERICAL PROCEDURES We solved the problem formulated in the previous section with two numerical schemes: a tracking interface technique proposed in this work, and an adaptation of the Volume-of- Fluid (VoF) method developed by Popinet [19, 20] and extended by L´opez-Herrera et al. [21] to study EHD phenomena. A. The proposed method Since we search for steady solutions, we drop any time dependence in the Navier-Stokes equations (2)-(4). As mentioned in Sec. II, we assume that the characteristic electric time is much shorter than any hydrodynamic time, and hence all the electrical charges are confined to the free surface (ρe= 0). In this case, the electric conductivity Kcan be regarded as constant. Because the liquid electrical permittivity is also assumed to be homogeneous, the electric forces vanish in the bulk domain, and enter into the problem through the Maxwell stresses in the free surface boundary conditions. A tracking interface technique was implemented to capture the steady liquid-gas interface position, r=F(z). We solved the nonlinear system of equations formed by the bulk equations (2)-(7), the interface equations (8)-(15), and the rest of the boundary conditions (16)-(17). In this way, one obtains the spatial dependence of the meniscus-jet shape F, the velocity and pressures fields (u, w, p) in the liquid domain, the surface charge density σe, and the inner and outer potential fields, ϕiand ϕo. The current intensity Iis the sum of the contributions due to the bulk conduction (Ib) and surface convection (Is), calculated at any axial position zalong the cone-jet as Ib(z) = 2π∫F(z) 0 αEiz(r, z)rdr, Is(z) = 2πσe(z)vs(z),(23) where vsis the dimensionless liquid velocity at the interface. All these quantities are determined as functions of the five dimensionless control parameters C,ve,β,χ, and α. A boundary-fitted coordinate system [25] was used to calculate the solution. Both the liquid and the gas domains were mapped onto the fixed rectangular domain (0 ≤η(l,g)≤ 1,0≤ξ≤ze) by means of the coordinate transformations (ηl=r/F, ξ =z) and (ηg= (r−F)/(re−F), ξ =z) for the liquid and gas domains, respectively. The liquid and gas domains were discretized using nηland nηgChebychev spectral collocation points [26] along 8 the ηland ηgaxes, respectively, while the axial direction ξwas discretized using nξuniformly distributed points. We used fourth-order central finite differences to get the axial derivatives except close to the ξ-boundaries, where second-order upwind finite differences were applied. All the simulations presented here were done with nηl= 11, nηg= 20, and nξ= 160. The equations in the discretized domain yielded a system of (4nηl+nηg)×nξ+ 2nξnon- linear equations which was solved iteratively using a trust-region method [27] implemented in the Matlab subroutine Fsolve. The initial guess for the iterations was obtained as follows. For given values of C,β,χ, and α, we first set a sufficiently high liquid velocity ve. Because a quasi-cylindrical jet is expected in this case, the equations were solved with u= 0, w= 2ve(1 −r2), F= 1, and ϕi=ϕe= 1 as the initial guess. Once the algorithm converged to a solution, the control parameter space was explored using a continuation technique. The algorithm showed itself to be robust and computationally efficient. It converged to the sought solution for all the cases considered in less than 1 hour running on an Intel Xeon(R) [email protected] processor. B. The VoF method We also performed simulations using the open source fluid simulation engine Gerris [28]. Gerris is a VoF method originally conceived to deal with the time-dependent incompressible Euler equations [19]. It has since been improved and extended considerably to suit problems arising in such diverse fields as multiphase [20], geostrophic [29], granular [30], and EHD [21] flows. The VoF method treats immiscible fluids as a single one with spatially varying properties through the interface. Therefore, the balance equation (1) for the surface charge density can not be explicitly considered. To circumvent this difficulty, the VoF method solves the evolution equation (A.8) (see Appendix) for the volumetric charge density ρe with K=const. It simulates a transient regime which eventually reaches the steady conejet state of electrospray (provided that this state is stable for the selected values of the governing parameters). We have verified that the steady states reached in the simulations were consistent with the leaky-dielectric model, i.e., all the charges were confined in a very thin layer next to the interface. In this regard, one can assert that the numerical method converged to the steady solution of the leaky-dielectric model. We set as initial condition a planar interface located inside the capillary tube. As time 9 0.0 0.2 0.4 0.6 0.8 1.0 1.2 4 3 2 1 0 I,F z- ( L n / R i ) FIG. 7. Axial distribution of the total intensity (circles) and the contributions associated with the bulk conduction (triangles) and the surface convection (squares) calculated by the proposed method for the experimental case with Q= 0.5 ml/h. The shape of the meniscus jet (solid line) is also plotted to show where the transition between conduction and convection of charges takes place. then electric forces result from the Maxwell stresses at that surface exclusively. The validity of the discretization scheme has been shown by comparing the results with those obtained with a much more sophisticated VoF method. The free surface position of the meniscusjet was measured very precisely in experiments with 1-octanol. The comparison between the numerical predictions for the liquid contour and the experimental results confirmed the validity of the proposed method. The numerical method also provided correct values of the emitted electric current even for configurations very close to the cone-jet stability limit. Appendix: Equations of the VoF method In this appendix, we derive the equations solved by the VoF method from the Taylor- Melcher leaky-dielectric model [23] . The starting point is the Maxwell electrostatic equations ∇ · (ϵE) = ρe,(A.1) ∇ × E= 0 ,(A.2) the Navier-Stokes equations ∇ · v= 0,(A.3) 16 ρDv Dt =−∇p−1 2E·E∇ϵ+∇ · (ϵE)E+µ∇2v,(A.4) and the equations for the conservation of neutral and charged species Dn(k) Dt =∇ · [−ω(k)ez(k)n(k)E+ω(k)kBT∇2n(k)] + r(k).(A.5) Here, n(k),r(k),z(k), and ω(k)are the concentration, production rate by chemical reaction, valence, and mobility of the kth species, respectively, while kB,e, and Tare the Boltzmann constant, the electron charge, and fluid temperature, respectively. The charge density ρecan be calculated as ρe=∑ k ez(k)n(k).(A.6) A conservation equation for this variable can be readily derived by adding the equations (A.5) for the different species: Dρe Dt =∇ · [−(∑ k ω(k)e2z(k)2n(k))E]+ekBT∑ k z(k)ω(k)∇2n(k)+∑ k ez(k)r(k).(A.7) In Eq. (A.7), migration due to thermal diffusion (the second term of the right-side of the equation) can be neglected as compared to electromigration, provided that ℓeEo/(kBT)≫1 (ℓand Eoare the characteristic length and electric field, respectively). In electrospray, ℓ∼(γε2 o/ρK2)1/3and Eo∼(γ/εoℓ)1/2[4], which yields ℓeEo/(kBT)≃8×103for 1-octanol. Therefore, thermal diffusion is negligible for the configurations considered in the present work. The production rates r(k)depend on the concentrations of the species and the chemical kinetics of the reactions occurring in the system. Note that these reactions can be extremely complex and, in some occasions, unknown. In fact, the rates r(k)depend on each other because the production/destruction of some species implies the destruction/production of others. Therefore, care in the modelling of the chemical reactions must be taken to prevent overdetermination. Nevertheless, three scenarios are simple enough to be easily modeled. The first is the unipolar injection, where a single species generated at an upstream electrode migrates through the fluid with zero production rate. The second one is the evolution of strong electrolytes. In this case, neutral species are fully dissociated and chemical reactions are absent (r(k)= 0). The last case is the evolution of a binary z-z electrolyte analyzed by 17 Saville [23]. Three species are diluted in the fluid: the neutral, the ionic, and the counterionic. In the forward reaction, the neutral species dissociate producing the same amount of ions and counterions, while in the recombination reaction the ionic and counter-ionic species disintegrate at the same rate to produce the neutral one. The source/sink term ∑kez(k)r(k) in Eq. (A.7) vanishes in the three cases described above. Taking into account these considerations, Eq. (A.7) reduces to the expression Dρe Dt =∇ · (−KE), K =∑ k ω(k)e2z(k)2n(k).(A.8) In the leaky-dielectric model [23], all characteristic times of hydrodynamic nature in the bulk are much longer than the electrical relaxation time. Therefore, the volumetric charge density ρeis assumed to be zero, and electrical charges are only allowed to exist at the interface. Therefore, Eq. (A.8) reduces to ∇ · (KE) = 0 .(A.9) It must be noted that Saville [23] keeps the bulk electric forces in Eq. (A.4) although, strictly speaking, they should be removed. Our VoF method solves Eqs. (A.1)-(A.4) and (A.9) with K=const. ACKNOWLEDGMENTS Partial support from the Ministry of Science and Education, Junta de Extremadura, and Junta de Andaluc´ıa (Spain) through Grants Nos. DPI2010-21103, GR10047, and P08-TEP- 04128, respectively, is gratefully acknowledged. [1] W. Gilbert, De Magnete (Book 2, Chapter 2, 1600). [2] A. Jaworek and A. Krupa, J. Aerosol Sci. 30, 873 (1999). [3] J. B. Fenn, M. Mann, C. K. Meng, S. F. Wong, and C. M. Whitehouse, Science 246, 64 (1989). [4] A. M. Ga˜n´an-Calvo, Phys. Rev. Lett. 79, 217 (1997). [5] J. F. de la Mora, Annu. Rev. Fluid Mech. 39, 217 (2007). [6] J. F. de la Mora and I. G. Loscertales, J . Fluid Mech. 260, 155 (1994). 18 [7] A. M. Ga˜n´an-Calvo, J. D´avila, and A. Barrero, J. Aerosol Sci. 28, 249 (1997). [8] A. M. Ga˜n´an-Calvo, J. Aerosol Sci. 30, 863 (1999). [9] A. M. Ga˜n´an-Calvo and J. M. Montanero, Phys. Rev. E 79, 066305 (2009). [10] R. P. A. Hartman, D. J. Brunner, D. M. A. Camelot, J. C. M. Marijnissen, and B. Scarlett, J. Aerosol Sci. 30, 823 (1999). [11] F. Yan, B. Farouk, and F. Ko, J. Aerosol. Sci. 34, 99 (2003). [12] F. J. Higuera, J. Fluid Mech. 484, 303 (2003). [13] G. Taylor, Proc. Roy. Soc. Lond. A, 453 (1969). [14] A. Barrero, A. Ga˜n´an-Calvo, J. D´avila, A. Palacio, and E. G´omez-Gonz´alez, J. Electrostat. 47, 13 (1999). [15] R. T. Collins, J. J. Jones, M. T. Harris, and O. A. Basaran, Nature Phys. 4, 149 (2008). [16] S. N. Reznik, A. L. Yarin, E. Zussman, and L. Bercovici, Phys. Fluids 18, 062101 (2006). [17] O. Lastow and W. Balachandran, J. Electrostatics 64, 850 (2006). [18] L. K. Lim, J. Hua, C.-H. Wang, and K. A. Smith, AICHE J. 57, 57 (2011). [19] S. Popinet, J. Comput. Phys. 190, 572 (2003). [20] S. Popinet, J. Comput. Phys. 228, 5838 (2009). [21] J. M. L´opez-Herrera, S. Popinet, and M. A. Herrada, J. Comput. Phys. 230, 1939 (2011). [22] J. M. L´opez-Herrera, P. Riesco-Chueca, and A. M. Ga˜n´an-Calvo, Phys. Fluids 17, 034106 (2005). [23] D. A. Saville, Annu. Rev. Fluid Mech. 29, 27 (1997). [24] A. M. Ga˜n´an-Calvo, J. C. Lasheras, J. D´avila, and A. Barrero, J. Aerosol Sci. 25, 22 (1994). [25] J. F. Thompson, F. C. Thames, and C. M. Mastin, J. Comput. Phys. 47, 1 (1982). [26] M. R. Khorrami, Int. J. Numer. Methods Fluids 12, 825 (1991). [27] J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research (Springer Verlag, New York, USA, 1999). [28] S. Popinet, “The Gerris flow solver,” http://gfs.sourceforge.net. [29] S. Popinet and G. Rickard, Ocean Modelling 16, 224 (2007). [30] L. Staron, P.-Y. Lagree, and S. Popinet, J. Fluid Mech. 686, 378 (2011). [31] M. G. Cabezas, A. Bateni, J. M. Montanero, and A. W. Neumann, Colloids Surf. A 255, 193 (2005). [32] J. M. Montanero, C. Ferrera, and V. M. Shevtsova, Exp. Fluids 45, 1087 (2008). 19 [33] E. J. Vega, J. M. Montanero, and C. Ferrera, Measurement 44, 1300 (2011). [34] N. Otsu, IEEE Transactions on Systems, Man, and Cybernetics (1979). [35] B. Song and J. Springer, J. Colloid Interface Sci. 184, 77 (1996). [36] I. Hayati, A. I. Bailey, and T. F. Tadros, Nature 319, 41 (1986). [37] A. Barrero, A. M. Ga˜n´an-Calvo, J. D´avila, A. Palacio, and E. G´omez-Gonz´alez, Phys. Rev. E58, 7309 (1998). [38] S. L. Anna, N. Bontoux, and H. A. Stone, Appl. Phys. Lett. 82, 364 (2003). [39] M. A. Herrada, A. M. Ga˜n´an-Calvo, A. Ojeda-Monge, B. Bluth, and P. Riesco-Chueca, Phys. Rev. E 78, 036323 (2008). [40] A. M. Ga˜n´an-Calvo, J. Fluid Mech. 507, 203 (2004). 20 View publication statsView publication stats