Journal of Colloid and Interface Science 417 (2014) 72–79
Contents lists available at ScienceDirect
Journal of Colloid and Interface Science www.elsevier.com/locate/jcis
Direct numerical simulation of AC dielectrophoretic particle–particle interactive motions Ye Ai a,⇑, Zhenping Zeng b, Shizhi Qian c,⇑ a
Pillar of Engineering Product Development, Singapore University of Technology and Design, Singapore 138682, Singapore School of Electronic and Optical Engineering, Nanjing University of Science & Technology, Nanjing 210094, China c Institute of Micro/Nanotechnology, Old Dominion University, Norfolk, VA 23529, USA b
a r t i c l e
i n f o
Article history: Received 6 September 2013 Accepted 11 November 2013 Available online 19 November 2013 Keywords: Arbitrary Lagrangian–Eulerian (ALE) Dielectrophoresis Microfluidics Particle chaining Particle assembly
a b s t r a c t Under an AC electric field, individual particles in close proximity induce spatially non-uniform electric field around each other, accordingly resulting in mutual dielectrophoretic (DEP) forces on these particles. The resulting attractive DEP particle–particle interaction could assemble individual colloidal particles or biological cells into regular patterns, which has become a promising bottom-up fabrication technique for bio-composite materials and microscopic functional structures. In this study, we developed a transient multiphysics model under the thin electric double layer (EDL) assumption, in which the fluid flow field, AC electric field and motion of finite-size particles are simultaneously solved using an Arbitrary Lagrangian–Eulerian (ALE) numerical approach. Numerical simulations show that negative DEP particle–particle interaction always tends to attract particles and form a chain parallel to the applied electric field. Particles usually accelerate at the first stage of the attractive motion due to an increase in the DEP interactive force, however, decelerate until stationary at the second stage due to a faster increase in the repulsive hydrodynamic force. Identical particles move at the same speed during the interactive motion. In contrast, smaller particles move faster than bigger particles during the attractive motion. The developed model explains the basic mechanism of AC DEP-based particle assembly technique and provides a versatile tool to design microfluidic devices for AC DEP-based particle or cell manipulation. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Dielectrophoresis refers to the motion of polarized particles under a spatially non-uniform electric field. Due to the simplicity in miniaturized integration of electric fields into microfluidics, dielectrophoresis has become one of the most popular techniques for particle manipulation in microfluidics [1]. One of the most interesting applications is directed assembly of individual colloidal particles or biological cells into microscopic functional structures. For example, electrically functional microwires can be fabricated from the assembly of metallic nanoparticles suspended in an aqueous medium [2]. Two-dimensional (2D) crystals have also been constructed from individual particles using a similar dielectrophoretic (DEP) assembly technique [3–5]. In addition, this technique has been exploited to assemble hematopoietic stem cells into a multilayered structure, mimicking the hematon for blood cell production [6]. Assembly of live biological cells and functionalized microparticles could open a new door for the development of novel tools for tissue engineering and biosensing [7,8].
⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Ai), [email protected] (S. Qian). 0021-9797/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcis.2013.11.034
Attractive chaining phenomenon of individual particles is the basis of DEP assembly technique. When an electric field is imposed in a particle suspension, the presence of particles locally distorts the electric field. If individual particles are close to each other, the non-uniform electric field around each particle becomes asymmetric with respect to its center, and accordingly induces a mutual DEP force acting on each other. If the mutual DEP force arising from the particle–particle interaction is attractive, particles tend to approach each other and eventually form a chain. Hwang et al. [9] performed an experimental study on the DEP particle–particle interaction, in which a pair of spherical particles, initially presenting an angle with respect to the applied electric field, were attracted to form a chain parallel to the electric field. Additionally, the alignment and assembly of two rod-shaped particles, as a result of the DEP particle–particle interaction, were also experimentally observed [10]. Besides growing experimental studies on the DEP assembly technique, further efforts have also been made on numerical simulations to gain an insightful understanding of the DEP particle– particle interaction. Aubry and her co-workers [11–14] developed a Lagrange multiplier-based numerical model to study the particle motion resulting from the DEP particle–particle interaction. The DEP force was calculated using a simplified point dipole (PD)
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
method in their model, which is sufficiently accurate when the gap between two particles is larger than the particle size [15]. Later, Kang and Li [16] investigated the DEP particle–particle interaction by balancing the DEP force and Stokes drag force. However, the approximation solution for the DEP particle–particle interactive force adopted in their study is also valid exclusively when the gap between particles is larger than the particle size. Maxwell stress tensor (MST) method has been demonstrated as the most rigorous approach for DEP force calculation [17–19]. Ai and Qian [20] adopted the MST method to account for the distortion of the electric field due to the presence of particles, and investigated particle–particle interactive motions under DC electric fields based on an arbitrary Lagrangian–Eulerian (ALE) method. It was found that the negative DEP particle–particle interaction always tends to chain particles parallel to the applied DC electric field. Recently, Ai and Qian’s model was further implemented by other numerical methods, such as smoothed profile method [21], boundary element method [22], and immersed boundary method [23], which predicted similar attractive motions of initially separated particles. In particular, House et al. [22] found that non-spherical particles could self-align with the applied electric field during the attractive chaining. It was also demonstrated that positive DEP particle–particle interaction tends to chain particles perpendicular to the applied electric field [23]. Note that all these studies considering finite-size particles with the MST method are under DC electric fields. However, AC electric fields are typically preferred in the application of DEP-based particle assembly, because electrophoresis of inherently charged particles can be eliminated in AC electric fields. In addition, it is feasible to switch the DEP motion from positive to negative, and vice versa, by tuning the frequency of AC electric fields. AC electric fields can also inhabit the electrolysis of water, a very common disadvantage for the DEP technique under DC electric fields. In this study, we developed a transient numerical model to investigate DEP particle–particle interactive motions under AC electric fields. The MST method is adopted to calculate the timeaveraged AC DEP force acting on each particle. The fluid flow field, AC electric field and particles’ motions are simultaneously solved using the ALE moving mesh algorithm. Transient positive and negative DEP motions of a single particle undergoing spatially nonuniform AC electric fields are first simulated to demonstrate the capability of the developed numerical model. AC DEP interactive motions of a pair of particles of identical and dissimilar sizes are subsequently studied to reveal the mechanism of AC DEP-based particle assembly technique.
2. Mathematical model We consider two circular particles suspended in an incompressible Newtonian fluid confined in a square with a side length of L, as illustrated in Fig. 1. The center of the square fluid coincides with the midpoint of the connecting line of the two particles, which is selected as the origin of the 2D Cartesian coordinate system (x, y). The radii of the particles initially located in the third quadrant and the first quadrant are b and c, respectively. The two particles are initially separated with a center-to-center distance of R, presenting an angle of h with respect to the x-axis. An AC electric field is applied along the x-axis to induce the DEP particle–particle interaction, and meanwhile excluding electrophoretic and electroosmotic effects in this study. It has been theoretically demonstrated that the effect of Brownian motion is usually negligible in short-range DEP particle–particle interactions [20]. In addition, our numerical method requires a finite separation distance between the two particles, which is larger than a typical thickness of electric double layer (EDL). For example, the EDL thickness of
73
Fig. 1. Two particles are suspended in a fluid medium under an externally applied ~ f . The origin of the Cartesian coordinate systems (x, y) is located at AC electric field E the midpoint of the connecting line of the two particles and also the center of the square fluid domain ABCD. The distance between the two particles and the angle between the connecting line of the two particles and the external electric field are, respectively, R and h. The fluid domain and particle domain are denoted by Xf and Xp, respectively. The surfaces of the particles in the first quadrant and the third quadrant are denoted by C and K, respectively.
a charged surface in a 0.1 M KCl solution at 25 °C is approximately 1 nm. Therefore, the EDL interactive force and van der Waals force are ignored in the present study. The net charge density in the entire computational domain is zero because of the thin EDL assumption. Therefore, the distribution of the quasi-static electric potential is governed the Gauss’s law, given by
~ f Þ ¼ 0 in Xf ; r ð~ef r/
ð1Þ
and
~ p Þ ¼ 0 in Xp ; r ð~ep r/
ð2Þ
where ~ef ¼ ef jrf =x and ~ep ¼ ep jrp =x are, respectively, the ~ f and / ~ p are, respeccomplex permittivity of the fluid and particle, / tively, the complex potential in the fluid and particle. In the above, ef and rf are the permittivity and conductivity of the fluid, respectively. Similarly, ep and rp are the permittivity and conductivity of the particle, respectively. x is the angular frequency of the AC elecpffiffiffiffiffiffiffi tric field. j ¼ 1 is the imaginary unit. The superscript ‘‘’’ represents complex variables. Electric potential applied to generate the AC electric field is given by
~ f ¼ /0 / 2
on AB;
ð3Þ
and
~ f ¼ /0 / 2
on CD:
ð4Þ
Electric insulation is applied on all the other boundaries
~ f ¼ 0 on AD and BC; n ~ef r/
ð5Þ
where n is the unit normal vector on the corresponding boundary. The electric potential and the normal component of the electric displacement are both continuous at the interface between the fluid and particle, given by
~f ¼ / ~p /
on K and C;
~ f ¼ n ~ep r/ ~p n ~ef r/
ð6Þ on K and C:
ð7Þ
74
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
The Reynolds number of the fluid flow in the present study is very small, therefore the fluid inertia can be neglected. As a result, the fluid flow is governed by the continuity equation and the Stokes equations, given by
r u ¼ 0 in Xf ;
ð8Þ
and
qf
@u ¼ rp þ gr2 u in Xf ; @t
ð9Þ
where u is the fluid velocity vector, qf is the fluid density, g is the dynamic viscosity, and p is the pressure. The fluid velocity on the particle surface is related to its rigidbody motion, expressed as
ui ¼ Upi þ xpi ðxsi xpi Þ on K and C;
ð10Þ
where Upi and xpi are, respectively, the translational velocity and rotational velocity of the ith particle, xsi and xpi are, respectively, the position vector of the surface and center of the ith particle. Segments AB and CD are assumed no slip boundary condition. Normal flow with zero pressure is imposed on segments AD and BC. The hydrodynamic force, FHi, and the time-averaged AC DEP force, FDEPi, acting on the ith particle are obtained, respectively, by integrating the hydrodynamic stress tensor, TH, and the timeaveraged MST [24], TM, over the surface of the ith particle, given by
FHi ¼
Z
TH ndSi ¼
Z h
T
i
pI þ gðru þ ðruÞ Þ ndSi ;
Z
TM ndSi ¼
Z
i ef h e e 0 e 0 e e f j2 I ndSi : ðEf Ef þ Ef Ef Þ jE 4
ð12Þ
e f is the electric field strength on the fluid side, obtained by Here E e 0 is the complex conjugate of E e f . The translation and ~ f ¼ r/ ~f . E E f rotation of the ith particle are determined by
dUpi ¼ FHi þ FDEPi ; dt
ð13Þ
and
dxpi ¼ dt
Z
ðxsi xpi Þ ½ðTH þ TM Þ ndSi ;
ð14Þ
where mpi and Ipi are, respectively, the mass and moment of inertia of the ith particle. Accordingly, the position and orientation of the particle can be obtained by solving the following equations,
dxpi ¼ Upi ; dt
ð15Þ
and
dhpi ¼ xpi : dt
ð16Þ
All the governing equations can be normalized by letting ~¼/ / ~ and t ¼ at =U 1 , where x ¼ ax ; u ¼ U 1 u ; p ¼ gU 1 p =a; / 1 the asterisk denotes dimensionless variables. The following quantities, a, U 1 ¼ ðef /21 Þ=ðgaÞ, and /1 = R0T/F, are respectively, the characteristic length, velocity, and electric potential, where R0 is the universal gas constant, T is the absolute temperature of the fluid medium, and F is the Faraday constant. Therefore, the dimensionless governing equations are expressed as
~ ¼ 0 in Xf ; r ~ef r / f
ð17Þ
~ Þ ¼ 0 in Xp ; r ð~ep r / p
ð18Þ
r u ¼ 0 in Xf ;
ð19Þ
in Xf ;
ð20Þ
where ~ef ¼ 1 jrf =ðxef Þ; ~ep ¼ 1 jrp =ðxep Þ, and Re = qfaU1/g. The dimensionless boundary conditions become
~ ¼ /0 / f 2/1
on AB;
~ ¼ /0 / f 2/1
ð21Þ
on CD;
ð22Þ
~ ¼ 0 on AD and BC; n ~ef r / f
ð23Þ
~ ¼ / ~ / f p
ð24Þ
on K and C;
~ ¼ n r / f
~ep ~ n r / p ~ef
on K and C;
ui ¼ Upi þ xpi ðxsi xpi Þ on K and C:
ð25Þ
ð26Þ
The equations governing the particle motion become
FHi ¼
Z
FDEPi ¼
FDEPi ¼
Ipi
@u ¼ r p þ r2 u @t
ð11Þ
and
mpi
Re
mpi
Ipi
TH ndSi ¼ Z
dUpi dt
dxpi dt
Z h
TM ndSi ¼
Z
i T p I þ ðr u þ ðr u Þ Þ ndSi ; i 1 h e e 0 e 0 e e j2 I ndS ; Ef Ef þ Ef Ef jE i f 4
¼ FHi þ FDEPi ;
¼
Z
ðxsi xpi Þ ðTH þ TM Þ n dSi ;
ð27Þ
ð28Þ
ð29Þ
ð30Þ
where the force, stress tensor, mass, moment of inertia, and rotational velocity are normalized as FHi ¼ gU 1 FHi ; TH ¼ gU 1 TH =a; mpi ¼ gampi =U 1 ; Ipi ¼ ga3 Ipi =U 1 , and xpi ¼ U 1 xpi =a. The presence of multiple particles induces a non-uniform electric field around them, and the resulting AC DEP effect gives rise to the particle motion. Meanwhile, the particle motion alters the surrounding flow field and electric field, which in turn affects the hydrodynamic and DEP forces and their corresponding torques acting on the particles. Therefore, the strongly coupled particle–fluid– electric field problem is numerically solved using the ALE method, which solves the fluid and electric fields in an Eulerian framework, and simultaneously tracks the particle motion in a Lagrangian manner. Detailed theoretical derivation of the ALE technique can be found in Ref. [25]. Numerical implementation of the mathematical model based on the ALE method is performed by a commercial finite-element package COMSOL (version 3.5a, www.comsol.com) operated with MATLAB (version 2009a, www.mathworks.com) in a high-performance cluster. This approach has been utilized to simulate pressure-driven particle motion in microchannels [25– 27], DC electrokinetic particle transport in microchannels [28– 32] and through nanopores [33,34]. Remarkable agreements with theoretical solutions, as well as various experimental results, have demonstrated the versatility and accuracy of this numerical approach [35]. In this study, quadratic triangular elements are generated in both fluid and particle domains. A denser mesh is created around the particles to accurately capture the nearby flow field and non-uniform electric field for precise calculation of the force and torque exerting on each particle. The total element number is typically around 15,000 to obtain converged and mesh-independent results.
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
75
3. Results and discussion The fluid medium and particles used in the simulations are water and polystyrene beads, respectively. The fluid density and viscosity are qf = 1.0 103 kg/m3 and g = 1.0 103 kg/(m s), respectively. To neglect the gravity effect, we set the density of the particle to be the same as that of the fluid. The conductivity and permittivity of the fluid are, respectively, rf = 2.0 102 S/m and ef = 80e0, where e0 is the absolute permittivity of vacuum. The conductivity and permittivity of the particle are, respectively, rp = 4.0 104 S/m and ep = 2.6e0 without special indication. The applied peak-to-peak potential and frequency are, respectively, /0 = 20 V and x/(2p) = 2 kHz. The characteristic length is the radius of the bigger particle, a = 5 lm, and the side length of the square fluid domain is L = 20a. 3.1. DEP motion of a single particle under non-uniform AC electric fields We first investigate the DEP motion of a single particle under non-uniform AC electric fields. In this section, the electric potential, /0/2, is selectively applied on the central one tenth of segment AB. The remaining portion of segment AB is specified as electric insulation. An electric potential, /0/2, is imposed along the entire segment CD, which accordingly creates a spatially non-uniform electric field across the entire fluid domain. The developed model is validated by comparing the time-averaged DEP forces acting on a single particle located in the center of the square with a side length of 20 lm, using the PD method and the MST method. The time-averaged DEP force of a circular particle calculated by the PD method is given by [36]
e rms j2 ; FPD ¼ 0:75pef R2p Reð~f CM Þrj E
ð31Þ
where Rp is the particle radius, ~f CM ¼ ð~ep ~ef Þ=ð~ep þ 2~ef Þ is the Claue rms j is the rms magnitude of the electric sius–Mossotti factor, and j E field within the fluid medium in the absence of particles. Reð~f CM Þ represents the real part of the complex Clausius–Mossotti factor, and a negative (positive) value represents a negative (positive) DEP motion, in which particles are pushed to the region with a lower (higher) electric field strength. The default properties of fluid and particle at an AC frequency of 2 kHz yield Reð~f CM Þ ¼ 0:485, indicating a negative DEP. Actually, polystyrene beads are always within the negative DEP region in a typical frequency range. In order to consider a positive DEP case, the conductivity of the particle is specified as rp = 0.2 S/m to get a positive value of Reð~f CM Þ ¼ 0:75. Fig. 2 shows the time-averaged DEP force as a function of the particle radius obtained by the two methods. Obviously, the DEP force increases accordingly as the particle radius increases. Smaller particles have a less distortion on the electric field than bigger particles. Therefore, the DEP forces calculated by the two methods are in agreement when the particle size is relatively small. Due to the distortion of the electric field by the particle, the difference between the two methods increases as the particle radius increases. The relative differences for positive and negative DEP force when the particle radius is 3 lm are 21.90% and 19.75%, respectively. In addition, the PD method overestimates the positive DEP force, and in contrast, underestimates the negative DEP force. When the particle size becomes comparable to the characteristic length of the channel or device, it is preferred to use the MST method for the calculation of the DEP force instead of the PD method. As the electric field is stronger near the electrode on segment AB, the negative (positive) DEP force tends to push the particle toward segment CD (AB). For the particle experiencing negative DEP motion and initially located at (x*, y*) = (8, 0), once it is released from rest, the negative DEP force accelerates the particle from sta-
Fig. 2. Time-averaged AC DEP force acting on a circular particle located at the origin as a function of the particle radius. Symbols and lines are the AC DEP force obtained by the PD method and MST method, respectively. Solid line and circles represent the positive DEP force; while dashed line and squares represent the negative DEP force.
tionary state to its steady velocity very rapidly (Fig. 3a). The negative DEP force pushes the particle away from the higher electric field region, where the non-uniformity of the electric field is also higher. Therefore, the x-component DEP force decreases as the particle moves further away from segment AB (Fig. 3b). Accordingly, the x-component particle velocity reaches a maximum after the initial rapid acceleration process and then gradually decreases (Fig. 3a). Note that the initial rapid acceleration process is purposely excluded from other velocity profiles. The electric field is symmetric with respect to the line of y = 0, therefore, the y-component DEP force is zero, and the particle only translates along the line of y = 0 (Fig. 3c and Movie 1 in the Supplementary material). For the particle with rp = 0.2 S/m and initially located at (x, y ) = (0, 0), it experiences positive DEP motion. Similarly, the particle moves along the line of y = 0 because of a zero y-component DEP force. The x-component DEP force increases as the particle is attracted to the higher electric field region (Fig. 3b). As a result, the x-component particle velocity gradually increases until the particle is about one particle radius away from segment AB. Subsequently, the x-component particle velocity starts to decrease quickly due to a rapid increase in the repulsive hydrodynamic force arising from the particle–wall interaction (Fig. 3a). Once the particle is eventually attached to the electrode, the particle is expected to become stationary. When the particle is initially placed at (x, y) = (0, 4), the positive DEP force also laterally attracts it toward the line of y = 0 during the movement toward segment AB (Fig. 3c and Movie 2 in the Supplementary material). It is worth mentioning that varying the electrical property of the fluid medium could also tune the Clausius–Mossotti factor and in turn the DEP force, and such effect is similar to the effect of the electrical property of the particle on its DEP motion. Conclusively, positive DEP motion tends to attract all suspended particles toward the shorter electrode. On the contrary, negative DEP motion tends to push particles away from the shorter electrode. If dissimilar particles or biological cells could be manipulated to have opposite DEP motions by controlling the properties of fluid medium and/or AC electric field, it is able to achieve highly selective particle separation or enrichment by AC DEP. The developed numerical model provides a versatile tool to design AC DEP-based particle or cell manipulation in microfluidics with complex geometries.
76
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
(a)
(b)
(c) Fig. 3. x-Component velocity (a), x-component time-averaged AC DEP force (b) and trajectory of a single particle subjected to a positive (solid line) and negative (dashed line) DEP motion. The particle experiencing negative DEP motion is initially located at (x, y) = (8, 0). The particle experiencing positive DEP motion is initially located at (x, y) = (0, 0) in (a and b), and (x, y) = (0, 4) in (c).
3.2. DEP particle–particle interaction: Parallel orientation In the following, we will investigate the relative motion of two particles arising from the negative DEP particle–particle interaction under an AC electric field. A uniform electric potential is imposed along the entire segment AB. Therefore, the non-uniform electric field is solely induced by the presence of particles. Under this condition, the DEP force calculated by the PD method is zero due to the neglect of the distortion of the electric field by the particle. Therefore, the PD method could not be used to calculate the DEP force in this study. Initially, two particles are placed at (x, y) = (±2.5, 0). When the two particles have an identical size (b = c = 1), they move toward each other at an identical speed straightly along their connecting line. Fig. 4a shows that the electric field around each individual particle is not symmetric with respect to the particle’s center, leading to a non-zero DEP force acting on each of them. The absence of either particle would break down the asymmetric distribution of the electric field around each particle. Hence, the DEP force purely arises from the electrical particle–particle interaction. As the electric field in the gap between the two particles is weakened, the negative DEP force results in an attractive motion, which is the
Fig. 4. Distribution of the electric field (a), fluid velocity, (b) and pressure, (c) around a pair of particles (b = c = 1) located at (x, y) = (±1.47, 0). Arrows in (b) represent the direction of the flow field. The negative DEP force shown in (a), FDEP, tends to attract the two particles to form a particle chain; while the pressure force denoted in (c), FP, repels the attractive motion.
basic mechanism of DEP-based particle chaining and assembly. The attractive motion of the two particles squeezes the fluids between them, and thus generates a symmetric circular flow field with respect to the origin (Fig. 4b). Additionally, the pressure between the two particles is enhanced, which retards the attractive motion (Fig. 4c). Fig. 5a shows the x-component velocity of the two particles along their corresponding trajectory. As mentioned previously, two identical particles (b = c = 1) move at an identical speed toward each other. At the first stage of the attractive motion, the
77
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
2
(a)
Starting point
y*
1
0
-1
Starting point -2 -2
-1
0
1
2
x* 3
(b) 2
y*
1
Fig. 5. x-Component velocity (a) and x-component time-averaged AC DEP force (b) of a pair of particles initially located at (x, y) = (±2.5, 0). Solid line, dashed line and dash-dotted line represent b = c = 1, b = 5c/4 = 1, and b = c = 4/5, respectively. Lines and symbols in (b) represent the magnitude of the DEP forces acting on the particles initially located at (x, y) = (2.5, 0) and (x, y) = (2.5, 0), respectively.
0
-1
-2
-3
particle velocity increases as they move closer to each other due to the increasing DEP force. After a certain traveling distance, the particle velocity starts to decrease due to the repulsive hydrodynamic pressure force that increases faster than the attractive DEP force thereafter. Therefore, the particle velocity decreases monotonously at the second stage of the attractive motion until forming a particle chain (Movie 3 in the Supplementary material). When the two particles are both reduced in size (b = c = 4/5), they move at an identical speed toward each other with a similar acceleration–deceleration process. However, the maximum particle velocity is smaller than that for the case of b = c = 1, and this is attributed to a lower DEP force on smaller particles. When the two particles are in different sizes (b = 5c/4 = 1), they still experience the acceleration–deceleration process during the attractive motion. However, it is found that the bigger particle moves more slowly than the smaller particle, indicating an asymmetric attractive motion of two particles of different sizes (Movie 4 in the Supplementary material). As a result, the final positions of the two particles are not symmetric with respect to the origin any more. Fig. 5b shows the x-component DEP force acting on the particle along its traveling time. It is intuitive to understand that the DEP forces acting on two particles of the same size (b = c = 4/5 and b = c = 1) are identical. As aforementioned, the mutual DEP interactive force arises from the breakdown of symmetric distribution of non-uniform electric field around each particle due to the presence of nearby particle. The bigger particle has a greater impact on the asymmetric distribution of electric field around the smaller particle. In other words, the electric field around the smaller particle is more asymmetrically non-uniform than that around
-3
-2
-1
0
1
2
3
x* Fig. 6. Trajectory of a pair of particles with an initial orientation of h = 45° (a) and h = 85° (b). Solid line, dashed line and dash-dotted line in (a) represent b = c = 1 (R = 4), b = 5c/4 = 1 (R = 4), and b = 5c/3 = 1 (R = 4), respectively. Solid line, dashed line and dash-dotted line in (b) represent b = c = 1 (R = 4), b = c = 1 (R = 3), and b = 5c/4 = 1 (R = 3), respectively.
the bigger particle. Therefore, two particles of different sizes (b = 5c/4 = 1) also experience an identical DEP force in the DEP particle–particle interaction. However, the bigger particle experiences a larger repulsive hydrodynamic force, resulting in a lower particle velocity. 3.3. DEP particle–particle interaction: Arbitrary orientation Generally, particles are randomly suspended in the fluid medium with any arbitrary orientation. Therefore, we further consider the negative DEP particle–particle interaction with a non-zero initial orientation with respect to the x-axis. Fig. 6a shows the trajectories of a pair of particles initially located at R = 4 and h = 45° with different particle sizes. When the two particles have an identical size (b = c = 1), the two particles start to rotate clockwise with respect to each other (Movie 5 in the Supplementary material). This process gradually reduces the orientation with respect to the x-axis until h = 0°. Eventually, the two particles end up with an attractive motion along the line of y = 0, just like they were initially located with a parallel orientation. It is found that the
78
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
trajectories of the two particles are actually antisymmetric with respect to the y-axis. When the size of the particle initially located in the first quadrant is reduced (b = 5c/4 = 1), the two particles also end up with an attractive motion, however, eventually along the line of y = 0.114. When the particle size is further reduced (b = 5c/3 = 1), the final position of the particle chain is along the line of y = 0.235, further shifted from the line of y = 0 (Movie 6 in the Supplementary material). Therefore, the difference in particle size breaks down the antisymmetry of particle trajectory with respect to the y-axis. It is also confirmed that a smaller particle moves faster than a bigger particle in the DEP particle–particle interaction. The initial orientation of the two particles with respect to the x-axis is further increased to h = 85°. As a result, the electric field in the gap between the two particles is actually enhanced. When b = c = 1 and R = 4, the DEP particle–particle interactive force first pushes the two particles away from each other, indicated by an increasing center-to-center distance between the two particles. However, the x-component DEP force still drives the two particles to rotate clockwise with respect to each other, reducing the orientation of the two particles with respect to the x-axis. As the two particles rotate further, the DEP particle–particle interaction reverts back to an attractive force. Therefore, the two particles start to approach each other, and eventually end up with the attractive motion along the line of y = 0 (Movie 7 in the Supplementary material). A very similar experimental study on the particle motion under an AC electric field was conducted recently [9], in which two spherical particles initially presented a non-zero angle with respect to the applied electric field. When an AC electric field was activated, the repulsive-to-attractive interaction between the two particles was observed, and eventually led to the formation of a chain parallel to the electric field. Our numerical results are in qualitative agreement with the experimental observations reported in the literature [9]. When the initial center-to-center distance is reduced to R = 3, the two particles make a smaller turn prior to the attractive motion along the line of y = 0. Antisymmetric motion of the two particles with respect to the y-axis is observed in both cases. When the particle initially located in the first quadrant is reduced in size (b = 5c/4 = 1), as expected, the antisymmetric particle motion is broken down (Movie 8 in the Supplementary material). When the connecting line of the two particles is perfectly perpendicular to the applied electric field, the electric field between the two particles is enhanced, and the resulting DEP particle– particle interactive force always pushes the two particles away from each other. However, a very small perturbation from the unavoidable Brownian motion can break down this extremely unstable orientation. Therefore, the DEP particle–particle interaction always tends to attract particles experiencing negative DEP motion to form a chain that is parallel to the applied electric field, independent of the initial particle orientation.
4. Concluding remarks Despite the fact that AC dielectrophoresis has been widely used for particle assembly in microfluidics, direct numerical simulations of the DEP motion of finite-size particles are only considered under DC electric fields [20–23]. In this work, a transient numerical model based on the ALE algorithm has been developed to investigate the DEP particle motion under AC electric fields in DEP-based particle assembly process [2–6]. The model accounts for the strongly coupled particle–fluid–electric field interactions and the AC DEP force is calculated by the MST method, considered as the most rigorous approach for DEP force calculation [17–19]. Numerical simulation of the motion of a single particle confirmed that negative
(positive) DEP motion tends to attract suspended particles toward the region with a lower (higher) electric field. Therefore, particles or cells with opposite DEP motions could be efficiently separated using the AC DEP effect. When two particles are close to each other, the asymmetric distortion of the surrounding electric field due to the presence of multiple particles gives rise to the mutual DEP particle–particle interaction. It was found that the negative AC DEP particle–particle interaction force always tends to chain particles parallel to the applied electric field, independent of the initial particle orientation with respect to the applied electric field, which is in good agreement with existing experimental results [9]. During the attractive motion, the two particles first accelerate, and subsequently decelerate due to a faster increase in the repulsive hydrodynamic pressure force as they further approach each other. When the two particles have the same size, they follow an identical trajectory pattern at the same speed. When the two particles have dissimilar size, the AC DEP forces acting on the two particles are still identical. However, the smaller particle moves faster than the bigger particle, and they follow different trajectory patterns. The final assembly site is slightly closer to the initial location of the bigger particle. Conclusively, the developed numerical model is of great value to understand the AC DEP-based particle assembly technique, and can be widely used to design microfluidic devices for AC DEP-based particle or cell manipulation. Acknowledgments This work was supported, in part, by SUTD-MIT International Design Center IDG11300101 (YA), the State Scholarship Fund of China under Grant No. 2011684502 (ZZ), and NSF DMS-1319078 (SQ).
Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jcis.2013.11.034. References [1] R. Pethig, Biomicrofluidics 4 (2010) 022811. [2] K.D. Hermanson, S.O. Lumsdon, J.P. Williams, E.W. Kaler, O.D. Velev, Science 294 (2001) 1082. [3] O.D. Velev, K.H. Bhatt, Soft Matter 2 (2006) 738. [4] S. Gangwal, O.J. Cayre, O.D. Velev, Langmuir 24 (2008) 13312. [5] S. Gangwal, A. Pawar, I. Kretzschmar, O.D. Velev, Soft Matter 6 (2010) 1413. [6] G.H. Markx, L. Carney, M. Littlefair, A. Sebastian, A.M. Buckle, Biomed. Microdevices 11 (2009) 143. [7] S. Gupta, R.G. Alargova, P.K. Kilpatrick, O.D. Velev, Soft Matter 4 (2008) 726. [8] OD. Velev, S. Gangwal, DN. Petsev, Ann. Rep. Sec. ‘‘C’’ (Phys. Chem.) 105 (2009) 213. [9] H. Hwang, J.J. Kim, J.K. Park, J. Phys. Chem. B 112 (2008) 9903. [10] M. Janjua, S. Nudurupati, I. Fischer, P. Singh, N. Aubry, Mech. Res. Commun. 36 (2009) 55. [11] A.T.J. Kadaksham, P. Singh, N. Aubry, Electrophoresis 25 (2004) 3625. [12] J. Kadaksham, P. Singh, N. Aubry, J. Fluids Eng. – Trans ASME 126 (2004) 170. [13] J. Kadaksham, P. Singh, N. Aubry, Mech. Res. Commun. 33 (2006) 108. [14] N. Aubry, P. Singh, Electrophoresis 27 (2006) 703. [15] N. Aubry, P. Singh, Europhys. Lett. 74 (2006) 623. [16] K.H. Kang, D.Q. Li, Langmuir 22 (2006) 1602. [17] A. Al-Jarro, J. Paul, D.W.P. Thomas, J. Crowe, N. Sawyer, F.R.A. Rose, K.M. Shakesheff, J. Phys. D – Appl. Phys. 40 (2007) 71. [18] C. Rosales, K.M. Lim, Electrophoresis 26 (2005) 2057. [19] X.J. Wang, X.B. Wang, P.R.C. Gascoyne, J. Electrostat. 39 (1997) 277. [20] Y. Ai, S. Qian, J. Colloid Interface Sci. 346 (2010) 448. [21] S. Kang, R. Maniyeri, J. Mech. Sci. Technol. 26 (2012) 3503. [22] D.L. House, H. Luo, S. Chang, J. Colloid Interface Sci. 374 (2012) 141. [23] M.R. Hossan, R. Dillon, A.K. Roy, P. Dutta, J. Colloid Interface Sci. 394 (2013) 619. [24] X. Wang, X.-B. Wang, P.R.C. Gascoyne, J. Electrostat. 39 (1997) 277. [25] H.H. Hu, N.A. Patankar, M.Y. Zhu, J. Comput. Phys. 169 (2001) 427. [26] Y. Ai, S.W. Joo, Y. Jiang, X. Xuan, S. Qian, Biomicrofluidics 3 (2009) 022404.
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79 [27] N. Al Quddus, W.A. Moussa, S. Bhattacharjee, J. Colloid Interface Sci. 317 (2008) 620. [28] Y. Ai, A. Beskok, D.T. Gauthier, S.W. Joo, S. Qian, Biomicrofluidics 3 (2009) 044110. [29] Y. Ai, S.W. Joo, Y. Jiang, X. Xuan, S. Qian, Electrophoresis 30 (2009) 2499. [30] Y. Ai, S. Park, J.J. Zhu, X.C. Xuan, A. Beskok, S.Z. Qian, Langmuir 26 (2010) 2937. [31] Y. Ai, S. Qian, S. Liu, S.W. Joo, Biomicrofluidics 4 (2010) 013201. [32] Y. Ai, B. Mauroy, A. Sharma, S.Z. Qian, Electrophoresis 32 (2011) 2282.
79
[33] Y. Ai, S. Qian, Electrophoresis 32 (2011) 996. [34] Y. Ai, S. Qian, Phys. Chem. Chem. Phys. 13 (2011) 4060. [35] S. Qian, Y. Ai, Electrokinetic Particle Transport in Micro/Nanofluidics: Direct Numerical Simulation Analysis, CRC Press Taylor & Francis Group, Boca Raton, FL, USA, 2012. [36] H. Morgan, N.G. Green, AC Electrokinetic: Colloids and Nanoparticles, Research Studies Press, Philadelphia, PA, USA, 2002.
Contents lists available at ScienceDirect
Journal of Colloid and Interface Science www.elsevier.com/locate/jcis
Direct numerical simulation of AC dielectrophoretic particle–particle interactive motions Ye Ai a,⇑, Zhenping Zeng b, Shizhi Qian c,⇑ a
Pillar of Engineering Product Development, Singapore University of Technology and Design, Singapore 138682, Singapore School of Electronic and Optical Engineering, Nanjing University of Science & Technology, Nanjing 210094, China c Institute of Micro/Nanotechnology, Old Dominion University, Norfolk, VA 23529, USA b
a r t i c l e
i n f o
Article history: Received 6 September 2013 Accepted 11 November 2013 Available online 19 November 2013 Keywords: Arbitrary Lagrangian–Eulerian (ALE) Dielectrophoresis Microfluidics Particle chaining Particle assembly
a b s t r a c t Under an AC electric field, individual particles in close proximity induce spatially non-uniform electric field around each other, accordingly resulting in mutual dielectrophoretic (DEP) forces on these particles. The resulting attractive DEP particle–particle interaction could assemble individual colloidal particles or biological cells into regular patterns, which has become a promising bottom-up fabrication technique for bio-composite materials and microscopic functional structures. In this study, we developed a transient multiphysics model under the thin electric double layer (EDL) assumption, in which the fluid flow field, AC electric field and motion of finite-size particles are simultaneously solved using an Arbitrary Lagrangian–Eulerian (ALE) numerical approach. Numerical simulations show that negative DEP particle–particle interaction always tends to attract particles and form a chain parallel to the applied electric field. Particles usually accelerate at the first stage of the attractive motion due to an increase in the DEP interactive force, however, decelerate until stationary at the second stage due to a faster increase in the repulsive hydrodynamic force. Identical particles move at the same speed during the interactive motion. In contrast, smaller particles move faster than bigger particles during the attractive motion. The developed model explains the basic mechanism of AC DEP-based particle assembly technique and provides a versatile tool to design microfluidic devices for AC DEP-based particle or cell manipulation. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Dielectrophoresis refers to the motion of polarized particles under a spatially non-uniform electric field. Due to the simplicity in miniaturized integration of electric fields into microfluidics, dielectrophoresis has become one of the most popular techniques for particle manipulation in microfluidics [1]. One of the most interesting applications is directed assembly of individual colloidal particles or biological cells into microscopic functional structures. For example, electrically functional microwires can be fabricated from the assembly of metallic nanoparticles suspended in an aqueous medium [2]. Two-dimensional (2D) crystals have also been constructed from individual particles using a similar dielectrophoretic (DEP) assembly technique [3–5]. In addition, this technique has been exploited to assemble hematopoietic stem cells into a multilayered structure, mimicking the hematon for blood cell production [6]. Assembly of live biological cells and functionalized microparticles could open a new door for the development of novel tools for tissue engineering and biosensing [7,8].
⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Ai), [email protected] (S. Qian). 0021-9797/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcis.2013.11.034
Attractive chaining phenomenon of individual particles is the basis of DEP assembly technique. When an electric field is imposed in a particle suspension, the presence of particles locally distorts the electric field. If individual particles are close to each other, the non-uniform electric field around each particle becomes asymmetric with respect to its center, and accordingly induces a mutual DEP force acting on each other. If the mutual DEP force arising from the particle–particle interaction is attractive, particles tend to approach each other and eventually form a chain. Hwang et al. [9] performed an experimental study on the DEP particle–particle interaction, in which a pair of spherical particles, initially presenting an angle with respect to the applied electric field, were attracted to form a chain parallel to the electric field. Additionally, the alignment and assembly of two rod-shaped particles, as a result of the DEP particle–particle interaction, were also experimentally observed [10]. Besides growing experimental studies on the DEP assembly technique, further efforts have also been made on numerical simulations to gain an insightful understanding of the DEP particle– particle interaction. Aubry and her co-workers [11–14] developed a Lagrange multiplier-based numerical model to study the particle motion resulting from the DEP particle–particle interaction. The DEP force was calculated using a simplified point dipole (PD)
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
method in their model, which is sufficiently accurate when the gap between two particles is larger than the particle size [15]. Later, Kang and Li [16] investigated the DEP particle–particle interaction by balancing the DEP force and Stokes drag force. However, the approximation solution for the DEP particle–particle interactive force adopted in their study is also valid exclusively when the gap between particles is larger than the particle size. Maxwell stress tensor (MST) method has been demonstrated as the most rigorous approach for DEP force calculation [17–19]. Ai and Qian [20] adopted the MST method to account for the distortion of the electric field due to the presence of particles, and investigated particle–particle interactive motions under DC electric fields based on an arbitrary Lagrangian–Eulerian (ALE) method. It was found that the negative DEP particle–particle interaction always tends to chain particles parallel to the applied DC electric field. Recently, Ai and Qian’s model was further implemented by other numerical methods, such as smoothed profile method [21], boundary element method [22], and immersed boundary method [23], which predicted similar attractive motions of initially separated particles. In particular, House et al. [22] found that non-spherical particles could self-align with the applied electric field during the attractive chaining. It was also demonstrated that positive DEP particle–particle interaction tends to chain particles perpendicular to the applied electric field [23]. Note that all these studies considering finite-size particles with the MST method are under DC electric fields. However, AC electric fields are typically preferred in the application of DEP-based particle assembly, because electrophoresis of inherently charged particles can be eliminated in AC electric fields. In addition, it is feasible to switch the DEP motion from positive to negative, and vice versa, by tuning the frequency of AC electric fields. AC electric fields can also inhabit the electrolysis of water, a very common disadvantage for the DEP technique under DC electric fields. In this study, we developed a transient numerical model to investigate DEP particle–particle interactive motions under AC electric fields. The MST method is adopted to calculate the timeaveraged AC DEP force acting on each particle. The fluid flow field, AC electric field and particles’ motions are simultaneously solved using the ALE moving mesh algorithm. Transient positive and negative DEP motions of a single particle undergoing spatially nonuniform AC electric fields are first simulated to demonstrate the capability of the developed numerical model. AC DEP interactive motions of a pair of particles of identical and dissimilar sizes are subsequently studied to reveal the mechanism of AC DEP-based particle assembly technique.
2. Mathematical model We consider two circular particles suspended in an incompressible Newtonian fluid confined in a square with a side length of L, as illustrated in Fig. 1. The center of the square fluid coincides with the midpoint of the connecting line of the two particles, which is selected as the origin of the 2D Cartesian coordinate system (x, y). The radii of the particles initially located in the third quadrant and the first quadrant are b and c, respectively. The two particles are initially separated with a center-to-center distance of R, presenting an angle of h with respect to the x-axis. An AC electric field is applied along the x-axis to induce the DEP particle–particle interaction, and meanwhile excluding electrophoretic and electroosmotic effects in this study. It has been theoretically demonstrated that the effect of Brownian motion is usually negligible in short-range DEP particle–particle interactions [20]. In addition, our numerical method requires a finite separation distance between the two particles, which is larger than a typical thickness of electric double layer (EDL). For example, the EDL thickness of
73
Fig. 1. Two particles are suspended in a fluid medium under an externally applied ~ f . The origin of the Cartesian coordinate systems (x, y) is located at AC electric field E the midpoint of the connecting line of the two particles and also the center of the square fluid domain ABCD. The distance between the two particles and the angle between the connecting line of the two particles and the external electric field are, respectively, R and h. The fluid domain and particle domain are denoted by Xf and Xp, respectively. The surfaces of the particles in the first quadrant and the third quadrant are denoted by C and K, respectively.
a charged surface in a 0.1 M KCl solution at 25 °C is approximately 1 nm. Therefore, the EDL interactive force and van der Waals force are ignored in the present study. The net charge density in the entire computational domain is zero because of the thin EDL assumption. Therefore, the distribution of the quasi-static electric potential is governed the Gauss’s law, given by
~ f Þ ¼ 0 in Xf ; r ð~ef r/
ð1Þ
and
~ p Þ ¼ 0 in Xp ; r ð~ep r/
ð2Þ
where ~ef ¼ ef jrf =x and ~ep ¼ ep jrp =x are, respectively, the ~ f and / ~ p are, respeccomplex permittivity of the fluid and particle, / tively, the complex potential in the fluid and particle. In the above, ef and rf are the permittivity and conductivity of the fluid, respectively. Similarly, ep and rp are the permittivity and conductivity of the particle, respectively. x is the angular frequency of the AC elecpffiffiffiffiffiffiffi tric field. j ¼ 1 is the imaginary unit. The superscript ‘‘’’ represents complex variables. Electric potential applied to generate the AC electric field is given by
~ f ¼ /0 / 2
on AB;
ð3Þ
and
~ f ¼ /0 / 2
on CD:
ð4Þ
Electric insulation is applied on all the other boundaries
~ f ¼ 0 on AD and BC; n ~ef r/
ð5Þ
where n is the unit normal vector on the corresponding boundary. The electric potential and the normal component of the electric displacement are both continuous at the interface between the fluid and particle, given by
~f ¼ / ~p /
on K and C;
~ f ¼ n ~ep r/ ~p n ~ef r/
ð6Þ on K and C:
ð7Þ
74
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
The Reynolds number of the fluid flow in the present study is very small, therefore the fluid inertia can be neglected. As a result, the fluid flow is governed by the continuity equation and the Stokes equations, given by
r u ¼ 0 in Xf ;
ð8Þ
and
qf
@u ¼ rp þ gr2 u in Xf ; @t
ð9Þ
where u is the fluid velocity vector, qf is the fluid density, g is the dynamic viscosity, and p is the pressure. The fluid velocity on the particle surface is related to its rigidbody motion, expressed as
ui ¼ Upi þ xpi ðxsi xpi Þ on K and C;
ð10Þ
where Upi and xpi are, respectively, the translational velocity and rotational velocity of the ith particle, xsi and xpi are, respectively, the position vector of the surface and center of the ith particle. Segments AB and CD are assumed no slip boundary condition. Normal flow with zero pressure is imposed on segments AD and BC. The hydrodynamic force, FHi, and the time-averaged AC DEP force, FDEPi, acting on the ith particle are obtained, respectively, by integrating the hydrodynamic stress tensor, TH, and the timeaveraged MST [24], TM, over the surface of the ith particle, given by
FHi ¼
Z
TH ndSi ¼
Z h
T
i
pI þ gðru þ ðruÞ Þ ndSi ;
Z
TM ndSi ¼
Z
i ef h e e 0 e 0 e e f j2 I ndSi : ðEf Ef þ Ef Ef Þ jE 4
ð12Þ
e f is the electric field strength on the fluid side, obtained by Here E e 0 is the complex conjugate of E e f . The translation and ~ f ¼ r/ ~f . E E f rotation of the ith particle are determined by
dUpi ¼ FHi þ FDEPi ; dt
ð13Þ
and
dxpi ¼ dt
Z
ðxsi xpi Þ ½ðTH þ TM Þ ndSi ;
ð14Þ
where mpi and Ipi are, respectively, the mass and moment of inertia of the ith particle. Accordingly, the position and orientation of the particle can be obtained by solving the following equations,
dxpi ¼ Upi ; dt
ð15Þ
and
dhpi ¼ xpi : dt
ð16Þ
All the governing equations can be normalized by letting ~¼/ / ~ and t ¼ at =U 1 , where x ¼ ax ; u ¼ U 1 u ; p ¼ gU 1 p =a; / 1 the asterisk denotes dimensionless variables. The following quantities, a, U 1 ¼ ðef /21 Þ=ðgaÞ, and /1 = R0T/F, are respectively, the characteristic length, velocity, and electric potential, where R0 is the universal gas constant, T is the absolute temperature of the fluid medium, and F is the Faraday constant. Therefore, the dimensionless governing equations are expressed as
~ ¼ 0 in Xf ; r ~ef r / f
ð17Þ
~ Þ ¼ 0 in Xp ; r ð~ep r / p
ð18Þ
r u ¼ 0 in Xf ;
ð19Þ
in Xf ;
ð20Þ
where ~ef ¼ 1 jrf =ðxef Þ; ~ep ¼ 1 jrp =ðxep Þ, and Re = qfaU1/g. The dimensionless boundary conditions become
~ ¼ /0 / f 2/1
on AB;
~ ¼ /0 / f 2/1
ð21Þ
on CD;
ð22Þ
~ ¼ 0 on AD and BC; n ~ef r / f
ð23Þ
~ ¼ / ~ / f p
ð24Þ
on K and C;
~ ¼ n r / f
~ep ~ n r / p ~ef
on K and C;
ui ¼ Upi þ xpi ðxsi xpi Þ on K and C:
ð25Þ
ð26Þ
The equations governing the particle motion become
FHi ¼
Z
FDEPi ¼
FDEPi ¼
Ipi
@u ¼ r p þ r2 u @t
ð11Þ
and
mpi
Re
mpi
Ipi
TH ndSi ¼ Z
dUpi dt
dxpi dt
Z h
TM ndSi ¼
Z
i T p I þ ðr u þ ðr u Þ Þ ndSi ; i 1 h e e 0 e 0 e e j2 I ndS ; Ef Ef þ Ef Ef jE i f 4
¼ FHi þ FDEPi ;
¼
Z
ðxsi xpi Þ ðTH þ TM Þ n dSi ;
ð27Þ
ð28Þ
ð29Þ
ð30Þ
where the force, stress tensor, mass, moment of inertia, and rotational velocity are normalized as FHi ¼ gU 1 FHi ; TH ¼ gU 1 TH =a; mpi ¼ gampi =U 1 ; Ipi ¼ ga3 Ipi =U 1 , and xpi ¼ U 1 xpi =a. The presence of multiple particles induces a non-uniform electric field around them, and the resulting AC DEP effect gives rise to the particle motion. Meanwhile, the particle motion alters the surrounding flow field and electric field, which in turn affects the hydrodynamic and DEP forces and their corresponding torques acting on the particles. Therefore, the strongly coupled particle–fluid– electric field problem is numerically solved using the ALE method, which solves the fluid and electric fields in an Eulerian framework, and simultaneously tracks the particle motion in a Lagrangian manner. Detailed theoretical derivation of the ALE technique can be found in Ref. [25]. Numerical implementation of the mathematical model based on the ALE method is performed by a commercial finite-element package COMSOL (version 3.5a, www.comsol.com) operated with MATLAB (version 2009a, www.mathworks.com) in a high-performance cluster. This approach has been utilized to simulate pressure-driven particle motion in microchannels [25– 27], DC electrokinetic particle transport in microchannels [28– 32] and through nanopores [33,34]. Remarkable agreements with theoretical solutions, as well as various experimental results, have demonstrated the versatility and accuracy of this numerical approach [35]. In this study, quadratic triangular elements are generated in both fluid and particle domains. A denser mesh is created around the particles to accurately capture the nearby flow field and non-uniform electric field for precise calculation of the force and torque exerting on each particle. The total element number is typically around 15,000 to obtain converged and mesh-independent results.
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
75
3. Results and discussion The fluid medium and particles used in the simulations are water and polystyrene beads, respectively. The fluid density and viscosity are qf = 1.0 103 kg/m3 and g = 1.0 103 kg/(m s), respectively. To neglect the gravity effect, we set the density of the particle to be the same as that of the fluid. The conductivity and permittivity of the fluid are, respectively, rf = 2.0 102 S/m and ef = 80e0, where e0 is the absolute permittivity of vacuum. The conductivity and permittivity of the particle are, respectively, rp = 4.0 104 S/m and ep = 2.6e0 without special indication. The applied peak-to-peak potential and frequency are, respectively, /0 = 20 V and x/(2p) = 2 kHz. The characteristic length is the radius of the bigger particle, a = 5 lm, and the side length of the square fluid domain is L = 20a. 3.1. DEP motion of a single particle under non-uniform AC electric fields We first investigate the DEP motion of a single particle under non-uniform AC electric fields. In this section, the electric potential, /0/2, is selectively applied on the central one tenth of segment AB. The remaining portion of segment AB is specified as electric insulation. An electric potential, /0/2, is imposed along the entire segment CD, which accordingly creates a spatially non-uniform electric field across the entire fluid domain. The developed model is validated by comparing the time-averaged DEP forces acting on a single particle located in the center of the square with a side length of 20 lm, using the PD method and the MST method. The time-averaged DEP force of a circular particle calculated by the PD method is given by [36]
e rms j2 ; FPD ¼ 0:75pef R2p Reð~f CM Þrj E
ð31Þ
where Rp is the particle radius, ~f CM ¼ ð~ep ~ef Þ=ð~ep þ 2~ef Þ is the Claue rms j is the rms magnitude of the electric sius–Mossotti factor, and j E field within the fluid medium in the absence of particles. Reð~f CM Þ represents the real part of the complex Clausius–Mossotti factor, and a negative (positive) value represents a negative (positive) DEP motion, in which particles are pushed to the region with a lower (higher) electric field strength. The default properties of fluid and particle at an AC frequency of 2 kHz yield Reð~f CM Þ ¼ 0:485, indicating a negative DEP. Actually, polystyrene beads are always within the negative DEP region in a typical frequency range. In order to consider a positive DEP case, the conductivity of the particle is specified as rp = 0.2 S/m to get a positive value of Reð~f CM Þ ¼ 0:75. Fig. 2 shows the time-averaged DEP force as a function of the particle radius obtained by the two methods. Obviously, the DEP force increases accordingly as the particle radius increases. Smaller particles have a less distortion on the electric field than bigger particles. Therefore, the DEP forces calculated by the two methods are in agreement when the particle size is relatively small. Due to the distortion of the electric field by the particle, the difference between the two methods increases as the particle radius increases. The relative differences for positive and negative DEP force when the particle radius is 3 lm are 21.90% and 19.75%, respectively. In addition, the PD method overestimates the positive DEP force, and in contrast, underestimates the negative DEP force. When the particle size becomes comparable to the characteristic length of the channel or device, it is preferred to use the MST method for the calculation of the DEP force instead of the PD method. As the electric field is stronger near the electrode on segment AB, the negative (positive) DEP force tends to push the particle toward segment CD (AB). For the particle experiencing negative DEP motion and initially located at (x*, y*) = (8, 0), once it is released from rest, the negative DEP force accelerates the particle from sta-
Fig. 2. Time-averaged AC DEP force acting on a circular particle located at the origin as a function of the particle radius. Symbols and lines are the AC DEP force obtained by the PD method and MST method, respectively. Solid line and circles represent the positive DEP force; while dashed line and squares represent the negative DEP force.
tionary state to its steady velocity very rapidly (Fig. 3a). The negative DEP force pushes the particle away from the higher electric field region, where the non-uniformity of the electric field is also higher. Therefore, the x-component DEP force decreases as the particle moves further away from segment AB (Fig. 3b). Accordingly, the x-component particle velocity reaches a maximum after the initial rapid acceleration process and then gradually decreases (Fig. 3a). Note that the initial rapid acceleration process is purposely excluded from other velocity profiles. The electric field is symmetric with respect to the line of y = 0, therefore, the y-component DEP force is zero, and the particle only translates along the line of y = 0 (Fig. 3c and Movie 1 in the Supplementary material). For the particle with rp = 0.2 S/m and initially located at (x, y ) = (0, 0), it experiences positive DEP motion. Similarly, the particle moves along the line of y = 0 because of a zero y-component DEP force. The x-component DEP force increases as the particle is attracted to the higher electric field region (Fig. 3b). As a result, the x-component particle velocity gradually increases until the particle is about one particle radius away from segment AB. Subsequently, the x-component particle velocity starts to decrease quickly due to a rapid increase in the repulsive hydrodynamic force arising from the particle–wall interaction (Fig. 3a). Once the particle is eventually attached to the electrode, the particle is expected to become stationary. When the particle is initially placed at (x, y) = (0, 4), the positive DEP force also laterally attracts it toward the line of y = 0 during the movement toward segment AB (Fig. 3c and Movie 2 in the Supplementary material). It is worth mentioning that varying the electrical property of the fluid medium could also tune the Clausius–Mossotti factor and in turn the DEP force, and such effect is similar to the effect of the electrical property of the particle on its DEP motion. Conclusively, positive DEP motion tends to attract all suspended particles toward the shorter electrode. On the contrary, negative DEP motion tends to push particles away from the shorter electrode. If dissimilar particles or biological cells could be manipulated to have opposite DEP motions by controlling the properties of fluid medium and/or AC electric field, it is able to achieve highly selective particle separation or enrichment by AC DEP. The developed numerical model provides a versatile tool to design AC DEP-based particle or cell manipulation in microfluidics with complex geometries.
76
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
(a)
(b)
(c) Fig. 3. x-Component velocity (a), x-component time-averaged AC DEP force (b) and trajectory of a single particle subjected to a positive (solid line) and negative (dashed line) DEP motion. The particle experiencing negative DEP motion is initially located at (x, y) = (8, 0). The particle experiencing positive DEP motion is initially located at (x, y) = (0, 0) in (a and b), and (x, y) = (0, 4) in (c).
3.2. DEP particle–particle interaction: Parallel orientation In the following, we will investigate the relative motion of two particles arising from the negative DEP particle–particle interaction under an AC electric field. A uniform electric potential is imposed along the entire segment AB. Therefore, the non-uniform electric field is solely induced by the presence of particles. Under this condition, the DEP force calculated by the PD method is zero due to the neglect of the distortion of the electric field by the particle. Therefore, the PD method could not be used to calculate the DEP force in this study. Initially, two particles are placed at (x, y) = (±2.5, 0). When the two particles have an identical size (b = c = 1), they move toward each other at an identical speed straightly along their connecting line. Fig. 4a shows that the electric field around each individual particle is not symmetric with respect to the particle’s center, leading to a non-zero DEP force acting on each of them. The absence of either particle would break down the asymmetric distribution of the electric field around each particle. Hence, the DEP force purely arises from the electrical particle–particle interaction. As the electric field in the gap between the two particles is weakened, the negative DEP force results in an attractive motion, which is the
Fig. 4. Distribution of the electric field (a), fluid velocity, (b) and pressure, (c) around a pair of particles (b = c = 1) located at (x, y) = (±1.47, 0). Arrows in (b) represent the direction of the flow field. The negative DEP force shown in (a), FDEP, tends to attract the two particles to form a particle chain; while the pressure force denoted in (c), FP, repels the attractive motion.
basic mechanism of DEP-based particle chaining and assembly. The attractive motion of the two particles squeezes the fluids between them, and thus generates a symmetric circular flow field with respect to the origin (Fig. 4b). Additionally, the pressure between the two particles is enhanced, which retards the attractive motion (Fig. 4c). Fig. 5a shows the x-component velocity of the two particles along their corresponding trajectory. As mentioned previously, two identical particles (b = c = 1) move at an identical speed toward each other. At the first stage of the attractive motion, the
77
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
2
(a)
Starting point
y*
1
0
-1
Starting point -2 -2
-1
0
1
2
x* 3
(b) 2
y*
1
Fig. 5. x-Component velocity (a) and x-component time-averaged AC DEP force (b) of a pair of particles initially located at (x, y) = (±2.5, 0). Solid line, dashed line and dash-dotted line represent b = c = 1, b = 5c/4 = 1, and b = c = 4/5, respectively. Lines and symbols in (b) represent the magnitude of the DEP forces acting on the particles initially located at (x, y) = (2.5, 0) and (x, y) = (2.5, 0), respectively.
0
-1
-2
-3
particle velocity increases as they move closer to each other due to the increasing DEP force. After a certain traveling distance, the particle velocity starts to decrease due to the repulsive hydrodynamic pressure force that increases faster than the attractive DEP force thereafter. Therefore, the particle velocity decreases monotonously at the second stage of the attractive motion until forming a particle chain (Movie 3 in the Supplementary material). When the two particles are both reduced in size (b = c = 4/5), they move at an identical speed toward each other with a similar acceleration–deceleration process. However, the maximum particle velocity is smaller than that for the case of b = c = 1, and this is attributed to a lower DEP force on smaller particles. When the two particles are in different sizes (b = 5c/4 = 1), they still experience the acceleration–deceleration process during the attractive motion. However, it is found that the bigger particle moves more slowly than the smaller particle, indicating an asymmetric attractive motion of two particles of different sizes (Movie 4 in the Supplementary material). As a result, the final positions of the two particles are not symmetric with respect to the origin any more. Fig. 5b shows the x-component DEP force acting on the particle along its traveling time. It is intuitive to understand that the DEP forces acting on two particles of the same size (b = c = 4/5 and b = c = 1) are identical. As aforementioned, the mutual DEP interactive force arises from the breakdown of symmetric distribution of non-uniform electric field around each particle due to the presence of nearby particle. The bigger particle has a greater impact on the asymmetric distribution of electric field around the smaller particle. In other words, the electric field around the smaller particle is more asymmetrically non-uniform than that around
-3
-2
-1
0
1
2
3
x* Fig. 6. Trajectory of a pair of particles with an initial orientation of h = 45° (a) and h = 85° (b). Solid line, dashed line and dash-dotted line in (a) represent b = c = 1 (R = 4), b = 5c/4 = 1 (R = 4), and b = 5c/3 = 1 (R = 4), respectively. Solid line, dashed line and dash-dotted line in (b) represent b = c = 1 (R = 4), b = c = 1 (R = 3), and b = 5c/4 = 1 (R = 3), respectively.
the bigger particle. Therefore, two particles of different sizes (b = 5c/4 = 1) also experience an identical DEP force in the DEP particle–particle interaction. However, the bigger particle experiences a larger repulsive hydrodynamic force, resulting in a lower particle velocity. 3.3. DEP particle–particle interaction: Arbitrary orientation Generally, particles are randomly suspended in the fluid medium with any arbitrary orientation. Therefore, we further consider the negative DEP particle–particle interaction with a non-zero initial orientation with respect to the x-axis. Fig. 6a shows the trajectories of a pair of particles initially located at R = 4 and h = 45° with different particle sizes. When the two particles have an identical size (b = c = 1), the two particles start to rotate clockwise with respect to each other (Movie 5 in the Supplementary material). This process gradually reduces the orientation with respect to the x-axis until h = 0°. Eventually, the two particles end up with an attractive motion along the line of y = 0, just like they were initially located with a parallel orientation. It is found that the
78
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79
trajectories of the two particles are actually antisymmetric with respect to the y-axis. When the size of the particle initially located in the first quadrant is reduced (b = 5c/4 = 1), the two particles also end up with an attractive motion, however, eventually along the line of y = 0.114. When the particle size is further reduced (b = 5c/3 = 1), the final position of the particle chain is along the line of y = 0.235, further shifted from the line of y = 0 (Movie 6 in the Supplementary material). Therefore, the difference in particle size breaks down the antisymmetry of particle trajectory with respect to the y-axis. It is also confirmed that a smaller particle moves faster than a bigger particle in the DEP particle–particle interaction. The initial orientation of the two particles with respect to the x-axis is further increased to h = 85°. As a result, the electric field in the gap between the two particles is actually enhanced. When b = c = 1 and R = 4, the DEP particle–particle interactive force first pushes the two particles away from each other, indicated by an increasing center-to-center distance between the two particles. However, the x-component DEP force still drives the two particles to rotate clockwise with respect to each other, reducing the orientation of the two particles with respect to the x-axis. As the two particles rotate further, the DEP particle–particle interaction reverts back to an attractive force. Therefore, the two particles start to approach each other, and eventually end up with the attractive motion along the line of y = 0 (Movie 7 in the Supplementary material). A very similar experimental study on the particle motion under an AC electric field was conducted recently [9], in which two spherical particles initially presented a non-zero angle with respect to the applied electric field. When an AC electric field was activated, the repulsive-to-attractive interaction between the two particles was observed, and eventually led to the formation of a chain parallel to the electric field. Our numerical results are in qualitative agreement with the experimental observations reported in the literature [9]. When the initial center-to-center distance is reduced to R = 3, the two particles make a smaller turn prior to the attractive motion along the line of y = 0. Antisymmetric motion of the two particles with respect to the y-axis is observed in both cases. When the particle initially located in the first quadrant is reduced in size (b = 5c/4 = 1), as expected, the antisymmetric particle motion is broken down (Movie 8 in the Supplementary material). When the connecting line of the two particles is perfectly perpendicular to the applied electric field, the electric field between the two particles is enhanced, and the resulting DEP particle– particle interactive force always pushes the two particles away from each other. However, a very small perturbation from the unavoidable Brownian motion can break down this extremely unstable orientation. Therefore, the DEP particle–particle interaction always tends to attract particles experiencing negative DEP motion to form a chain that is parallel to the applied electric field, independent of the initial particle orientation.
4. Concluding remarks Despite the fact that AC dielectrophoresis has been widely used for particle assembly in microfluidics, direct numerical simulations of the DEP motion of finite-size particles are only considered under DC electric fields [20–23]. In this work, a transient numerical model based on the ALE algorithm has been developed to investigate the DEP particle motion under AC electric fields in DEP-based particle assembly process [2–6]. The model accounts for the strongly coupled particle–fluid–electric field interactions and the AC DEP force is calculated by the MST method, considered as the most rigorous approach for DEP force calculation [17–19]. Numerical simulation of the motion of a single particle confirmed that negative
(positive) DEP motion tends to attract suspended particles toward the region with a lower (higher) electric field. Therefore, particles or cells with opposite DEP motions could be efficiently separated using the AC DEP effect. When two particles are close to each other, the asymmetric distortion of the surrounding electric field due to the presence of multiple particles gives rise to the mutual DEP particle–particle interaction. It was found that the negative AC DEP particle–particle interaction force always tends to chain particles parallel to the applied electric field, independent of the initial particle orientation with respect to the applied electric field, which is in good agreement with existing experimental results [9]. During the attractive motion, the two particles first accelerate, and subsequently decelerate due to a faster increase in the repulsive hydrodynamic pressure force as they further approach each other. When the two particles have the same size, they follow an identical trajectory pattern at the same speed. When the two particles have dissimilar size, the AC DEP forces acting on the two particles are still identical. However, the smaller particle moves faster than the bigger particle, and they follow different trajectory patterns. The final assembly site is slightly closer to the initial location of the bigger particle. Conclusively, the developed numerical model is of great value to understand the AC DEP-based particle assembly technique, and can be widely used to design microfluidic devices for AC DEP-based particle or cell manipulation. Acknowledgments This work was supported, in part, by SUTD-MIT International Design Center IDG11300101 (YA), the State Scholarship Fund of China under Grant No. 2011684502 (ZZ), and NSF DMS-1319078 (SQ).
Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jcis.2013.11.034. References [1] R. Pethig, Biomicrofluidics 4 (2010) 022811. [2] K.D. Hermanson, S.O. Lumsdon, J.P. Williams, E.W. Kaler, O.D. Velev, Science 294 (2001) 1082. [3] O.D. Velev, K.H. Bhatt, Soft Matter 2 (2006) 738. [4] S. Gangwal, O.J. Cayre, O.D. Velev, Langmuir 24 (2008) 13312. [5] S. Gangwal, A. Pawar, I. Kretzschmar, O.D. Velev, Soft Matter 6 (2010) 1413. [6] G.H. Markx, L. Carney, M. Littlefair, A. Sebastian, A.M. Buckle, Biomed. Microdevices 11 (2009) 143. [7] S. Gupta, R.G. Alargova, P.K. Kilpatrick, O.D. Velev, Soft Matter 4 (2008) 726. [8] OD. Velev, S. Gangwal, DN. Petsev, Ann. Rep. Sec. ‘‘C’’ (Phys. Chem.) 105 (2009) 213. [9] H. Hwang, J.J. Kim, J.K. Park, J. Phys. Chem. B 112 (2008) 9903. [10] M. Janjua, S. Nudurupati, I. Fischer, P. Singh, N. Aubry, Mech. Res. Commun. 36 (2009) 55. [11] A.T.J. Kadaksham, P. Singh, N. Aubry, Electrophoresis 25 (2004) 3625. [12] J. Kadaksham, P. Singh, N. Aubry, J. Fluids Eng. – Trans ASME 126 (2004) 170. [13] J. Kadaksham, P. Singh, N. Aubry, Mech. Res. Commun. 33 (2006) 108. [14] N. Aubry, P. Singh, Electrophoresis 27 (2006) 703. [15] N. Aubry, P. Singh, Europhys. Lett. 74 (2006) 623. [16] K.H. Kang, D.Q. Li, Langmuir 22 (2006) 1602. [17] A. Al-Jarro, J. Paul, D.W.P. Thomas, J. Crowe, N. Sawyer, F.R.A. Rose, K.M. Shakesheff, J. Phys. D – Appl. Phys. 40 (2007) 71. [18] C. Rosales, K.M. Lim, Electrophoresis 26 (2005) 2057. [19] X.J. Wang, X.B. Wang, P.R.C. Gascoyne, J. Electrostat. 39 (1997) 277. [20] Y. Ai, S. Qian, J. Colloid Interface Sci. 346 (2010) 448. [21] S. Kang, R. Maniyeri, J. Mech. Sci. Technol. 26 (2012) 3503. [22] D.L. House, H. Luo, S. Chang, J. Colloid Interface Sci. 374 (2012) 141. [23] M.R. Hossan, R. Dillon, A.K. Roy, P. Dutta, J. Colloid Interface Sci. 394 (2013) 619. [24] X. Wang, X.-B. Wang, P.R.C. Gascoyne, J. Electrostat. 39 (1997) 277. [25] H.H. Hu, N.A. Patankar, M.Y. Zhu, J. Comput. Phys. 169 (2001) 427. [26] Y. Ai, S.W. Joo, Y. Jiang, X. Xuan, S. Qian, Biomicrofluidics 3 (2009) 022404.
Y. Ai et al. / Journal of Colloid and Interface Science 417 (2014) 72–79 [27] N. Al Quddus, W.A. Moussa, S. Bhattacharjee, J. Colloid Interface Sci. 317 (2008) 620. [28] Y. Ai, A. Beskok, D.T. Gauthier, S.W. Joo, S. Qian, Biomicrofluidics 3 (2009) 044110. [29] Y. Ai, S.W. Joo, Y. Jiang, X. Xuan, S. Qian, Electrophoresis 30 (2009) 2499. [30] Y. Ai, S. Park, J.J. Zhu, X.C. Xuan, A. Beskok, S.Z. Qian, Langmuir 26 (2010) 2937. [31] Y. Ai, S. Qian, S. Liu, S.W. Joo, Biomicrofluidics 4 (2010) 013201. [32] Y. Ai, B. Mauroy, A. Sharma, S.Z. Qian, Electrophoresis 32 (2011) 2282.
79
[33] Y. Ai, S. Qian, Electrophoresis 32 (2011) 996. [34] Y. Ai, S. Qian, Phys. Chem. Chem. Phys. 13 (2011) 4060. [35] S. Qian, Y. Ai, Electrokinetic Particle Transport in Micro/Nanofluidics: Direct Numerical Simulation Analysis, CRC Press Taylor & Francis Group, Boca Raton, FL, USA, 2012. [36] H. Morgan, N.G. Green, AC Electrokinetic: Colloids and Nanoparticles, Research Studies Press, Philadelphia, PA, USA, 2002.
