Immersed electrokinetic finite element method

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2007; 71:379–405 Published online 13 December 2006 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.1941 Immersed electrokinetic finite element method Yaling Liu, Wing Kam Liu∗, † , Ted Belytschko, Neelesh Patankar, Albert C. To, Adrian Kopacz and Jae-Hyun Chung‡ Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3111, U.S.A. SUMMARY A new method is proposed for modelling the electrokinetic-induced mechanical motion of particles in a fluid domain under an applied electric field. In this method, independent solid meshes move in a fixed background field mesh that models the fluid and the electric field. This simple strategy removes the need for expensive mesh updates. Furthermore, the reproducing kernel particle functions enable efficient coupling of various immersed deformable solids with the surrounding viscous fluid in the presence of an applied electric field. The electric force on a particle is calculated by the Maxwell stress tensor method. For the first time, three-dimensional assembly of nano/biomaterials of various geometries and electrical properties have been comprehensively studied using the new method. Simulation of the dynamic process of electro-manipulation of individual and multiple cells agrees well with experimental data. Preliminary results for selective deposition of viruses and stretching of a DNA molecule are also presented. Copyright q 2006 John Wiley & Sons, Ltd. Received 22 July 2006; Accepted 23 October 2006 KEY WORDS: immersed finite element method; electrokinetics; dielectrophoretic assembly; biomolecules; nanowires; Euler–Lagrange mapping 1. INTRODUCTION Assembly and manipulation of nano/biomaterials are of great interest and have become feasible with the development of advanced micro/nanoscale devices. Many methods have been proposed to pattern nano/biomaterials such as nanowires (NWs), cells, and biomolecules. For example, chemical binding has been used to assemble single molecules [1]. Magnetic fields have been utilized to orient single-walled carbon nanotubes (CNTs) into a densely packed array [2]. ∗ Correspondence to: Wing Kam Liu, Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3111, U.S.A. w-liu@northwestern.edu Current address: University of Washington, Department of Mechanical Engineering, Campus Box: 352600, Seattle, U.S.A. † E-mail: ‡ Copyright q 2006 John Wiley & Sons, Ltd. 380 Y. LIU ET AL. Electro-optical forces have been proposed to align gold nano-rods [3]. Mechanical flow has also been shown to be an efficient means for the assembly of NWs onto electrodes [4]. Forces resulting from an electric field can be advantageous for the manipulation of micro/nanoscale particles because both the amplitude and the direction of forces can be controlled by the applied electric field strength and frequency [5]. These forces were investigated as micromachining technology emerged in the 1980s, because electric field strengths over 1 kV/cm could readily be obtained at the microscale. For example, Trau et al. [6] used a direct current (DC) field-induced electrohydrodynamic flow to assemble the colloidal crystals on electrode surfaces. Brisson and Tilton [7] utilized a similar electrophoretic deposition method to form dense, quasiordered two-dimensional cell clusters adjacent to the electrode surface. These patterned cell arrays are reversible, i.e. they dissipate by diffusion on removal of the electric field. Besides the massive assembly of particles, there is growing interest in the manipulation of individual flexible biopolymers or cells. Zimmermann et al. [8] used high-intensity alternating current (AC) field to manipulate mammalian cells. By applying a 2 MHz AC field of 2 kV/cm strength, they observed the electrodeformation of human erythrocytes myeloma into a stretched shape. The hybridoma cells aligned vertically between electrodes by an AC field between two cylindrical electrodes. Recently, more precise assembly methods were proposed. Washizu et al. [9] used the electro-osmotic flow to immobilize and stretch the DNA molecules. Chung et al. [10] and Lee et al. [5] assembled the multi-walled CNTs and DNA by combining AC and DC electric fields. This technology allows the control and deposit of individual nanoscale components on micro/nano systems. Despite the many emerging applications of the various assembly procedures, the fundamental mechanisms are not fully understood due to the complexity of the assembly processes at the micro/nanoscale. The underlying mechanisms of electric field-driven assembly are yet to be elucidated; for example, (1) how the electric forces contribute to the assembly processes; (2) how the experimental conditions such as the amplitude and frequency of the electric field, electrode geometry, and material properties affect the assembly process; (3) how large-scale assembly can be achieved in a precise manner by the electric field-guided methods. Patankar and Hu [11] obtained a three-dimensional, steady-state numerical solution of electroosmotic flow using the finite volume method. Legay et al. [12] have developed a fluid–structure interaction method based on level set and the extended finite element method. The assembly of colloidal particles on electrode surfaces by electro-osmotic flow was studied by Solomentsev et al. [13]. Their model is based on an analytical solution in bispherical co-ordinates and is limited to a few rigid circular particles with simple configurations. Dimaki and Boggild [14] simulated the dielectrophoretic (DEP) deposition of CNTs. However, their simulation is based on the linear effective dipole moment (EDM) approximation of the DEP force on ellipsoid particles, and the drag force in their simulation is not obtained from a direct hydrodynamic calculation. The EDM theory is valid only when the dimensions of the particle are much smaller than the characteristic length scale of the electric field gradient, and thus, it is inaccurate when the CNTs are close to the electrodes or are bundled together. Moreover, only few analytical solutions, based on the EDM theory for simple geometries such as spheres or prolate/oblate particles, are available. The EDM theory is currently not applicable when particles have complex geometries or are deformable. Currently, no methodology is available for modelling assemblies of flexible, arbitrarily shaped nano/biomaterials in fluid environments with complex electrode geometries, under an electric field. For example, the characterization of the electrodeformation and alignment of cells is hindered by the inability to model the complex coupling of the electric field, cell deformation, and fluid flow. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 381 Here, the immersed finite element method (IFEM) [15–23] is extended to incorporate electrohydrodynamics and electrokinetic forces on solids to become the immersed electrokinetic finite element method (IEFEM). The method couples the electric field with hydrodynamics to study the motion and deformation of flexible objects immersed in a fluid under an applied electric field. The electric field-induced force on an arbitrary object is calculated based on the continuum electromechanics and the Maxwell stress tensor (MST) method [24], and is coupled with the fluid/solid motion. This coupled hydro-electrokinetic formulation can capture the motion and deformation of multiple flexible, arbitrarily shaped continua in complex environments including features such as complex electrode geometries and composite DC/AC fields. The proposed direct numerical method is important for modelling electric field-guided assembly of nano/bioparticles, in which the EDM theory breaks down near the electrodes and interaction between particles has to be considered. The dynamic processes of attraction, alignment, and deposition of immersed objects between micro-electrodes are modelled by integrating the electric force, fluid flow, and solid motion/deformation. This method is applied to simulate the assembly of NWs and can also be used to explore the mechanisms of the electro-manipulation of cells and biomolecules such as DNA and viruses. We will validate our method by comparing our calculations with experimental results. 2. ELECTROMECHANICS In this section, we describe two different approaches for electric field calculation; the appropriate method depends on the current carriers in the medium. When the inherent carriers such as electrons are the major current carriers, we adopt a complex variable formulation for the electric potential. On the other hand, if ions are the major current carriers, we solve the reduced Maxwell equations together with a species conservation equation. The nomenclature used in this paper is listed in Table I. 2.1. Electric field calculation We assume that the material properties, i.e. electrical conductivity and permittivity, are constant in the fluid and the solids. The time scale of the motion of the solids is assumed to be much larger than that of the electric field. For currents and frequencies typically found in electric assembly processes, the governing equations of the electric field can be described by the reduced quasi-electrostatic form of the Maxwell equations, and thus, the respective Gauss equation and the charge conservation equation apply [25, 26]: ∇ · D =
e D
e +∇ ·J=0 Dt (1) (2) where
e is the free bulk charge density, J is the conduction current, and D is the electric flux density or the electric displacement vector. The electric field E is assumed to be irrotational E = −∇ Copyright q 2006 John Wiley & Sons, Ltd. (3) Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 382 Y. LIU ET AL. Table I. List of principal quantities. Quantities Mechanical Spatial coordinate Displacement Velocity Stress Symbol Units x u v cm cm cm/s g/(cm s2 ) r Electrical Electric field potential Electric field AC field angular frequency AC field frequency Conductivity Permittivity Complex permittivity Complex electric field Ion concentration Free charge density Electric flux density Conduction current Valence of the kth species Mobility of the kth species  V V/m rad/s Hz S/m F/m F/m V/m m−3 C/m−3 C/m2 C/m2 s E  f   ˜ (=  + /j ) Ẽ ck
e D J zk Name S: Siemens F: Farads C: Coulomb m2 /s/V k For a homogeneous linear dielectric material with permittivity  and conductivity , J = E, D = E. Assuming that the time scale of particle motion is much longer than the electric time scale, the material time derivative in (2) can be approximated by the partial time derivative D
e *
e *
e = + v · ∇
e ≈ Dt *t *t (4) where v is the particle velocity. Hence, Equations (1) and (2) become ∇ · (E) =
e *
e + ∇ · (E) = 0 *t (5) (6) When objects with different conductivities and permittivities are immersed in the medium, the following jump conditions hold at the interfaces: n · 'E( = s e n · 'E( = − (7) *s e *t (8) where s e is the free surface charge and 'A( denotes A|n+ − A|n− . Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 383 Depending on whether the major current carriers are inherent carriers or ions, we adopt two different approaches for the electric field calculation. Case 1: When ions are major current carriers. If the ions are the major current carriers and the solids are insulating materials, the ions inside the fluid move under an electric field and their concentration is influenced by the fluid flow. The free charge density in (5) is related to the ion concentration by

e = ez k ck (9) k where e is the charge on a proton and z k is the valence of the kth species whose concentration is ck . The species conservation equation is given by [27] *ck (10) + vf · ∇ck = ∇ · [−k ez k ck E + k kB T ∇ck ], k = 1, . . . , N *t where kB is the Boltzmann constant, T is the temperature, and k is the mobility of the kth species. The first term on the right-hand side describes ion migration under electric field and the second term represents ion transport through diffusion. Equation (10) is solved to obtain the concentration for the kth species ck , which can then be used to compute the free charge density in (9) and consequently the electric potential in (5). Case 2: When inherent carriers are major current carriers. When the inherent carriers are the major current carriers, the effect of ions on the current is negligible. A complex variable formulation can be employed for the electric field. The electric potential for an AC field with frequency  is assumed to be of the form ˜ (x, t) = M (x)ejwt (11) Using *
e /*t = j
e , substituting Equation (5) into Equation (6), and adopting notation of complex permittivity ˜ =  + (/j), the complex-coupled free charge/current Maxwell equation (1), (2), (7), (8) can be reduced to a simpler form in terms of complex variables ˜ =0 ∇ · (˜∇ ) (12) n · '˜Ẽ( = 0 (13) 2.2. Maxwell stress tensor The electrostatic and hydrodynamic effects can be coupled through the MST. The DEP force arises from induced dipole moments embedded in a non-uniform electric field. The force exerted by an electric field E on a dipole with dipole moment p is given by p · ∇E and the force due to a free charge is
e E. Using the Maxwell equations, the total electrostatic force per unit volume can be written as
e E + p · ∇E = (∇ · E)E − (E) · ∇E = ∇ · (E ⊗ E − 12 (E · E)I) ≡ ∇ · TM Copyright q 2006 John Wiley & Sons, Ltd. (14) Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 384 Y. LIU ET AL. Figure 1. An illustration of the order of magnitude of various factors on the motion of a spherical particle 10
m above two planar electrodes with a spacing of 20
m, density = 1 g/cm3 , viscosity = 0.78 Pa s. where the Maxwell stress has been defined as TM ≡ E ⊗ E − 21 (E · E)I (15) The electric force  on an object can be calculated through a surface integral of the Maxwell stress tensor Fe = s TM · n dA, or as a volume integral Fe = s ∇ · TM dV . Because the time scale of the motion of the particles is much larger than the time scale of the AC frequency, the electric force used in the equation of motion of the particles should be the average force over the period of the AC field. We write an AC field E(x, t) as the sum of a complex variable Ẽ(x, t) = EM (x)ejt and its conjugate Ẽ∗ as E(x, t) = Re(EM (x)ejt ) = 12 (Ẽ(x, t) + Ẽ∗ (x, t)) Then, the time-averaged MST is expressed as [24]  1 T M 1 M r = T dt = Re(˜)((Ẽ ⊗ Ẽ∗ + Ẽ∗ ⊗ Ẽ) − |Ẽ|2 I) T 0 4 (16) (17) where the period of the electric field T = 2/. 2.3. Order of magnitude analysis of various factors Factors other than electrical and viscous forces may also influence the particle motion. We plot the particle displacement per second under various factors and for different particle sizes in Figure 1. In our simulation, we focus on particles with size between 100 nm and 1
m, under which the particle motion is dominantly controlled by electric forces. 3. SIMULATION METHOD Consider a fluid domain  enclosed by a sufficiently smooth boundary f , as shown in Figure 2. Suppose there exists a submerged solid domain s enclosed by a boundary s ; the entire domain Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 385 IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD Figure 2. Configuration of the problem.  is subdivided into two regions, namely the solid region s and the fluid region \s . The interface between the solid and fluid regions is denoted by s . In the fluid domain, the Eulerian description is used while the Lagrangian description is employed in the solid domain. The electric field is calculated in the entire domain . 3.1. Strong form of electromechanics 3.1.1. Electric field. If ions in the fluid are the major charge carriers and the solids are insulators, we solve ∇ · (s ∇s ) = 0 in s (18) ck = 0 in s (19) ∇ · (f ∇f ) =
e in f

e = ez k ck , f (20) (21) k *ck + vf · ∇ck = ∇ · [−k ez k ck E + k kB T ∇ck ], *t k = 1, . . . , N in f (22) The interface between the solid and the fluid is updated along with the solid domain in the Lagrangian description. The electric field can be obtained by solving Equations (18)–(22). If the inherent carriers are the major current carriers, from Equation (12), the equations for the electric field are described in complex variables as ˜ f) = 0 ∇ · (˜f ∇  in f (23) ˜ s) = 0 ∇ · (˜s ∇  in s (24) where the superscripts s and f denote solid and fluid, respectively. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 386 Y. LIU ET AL. From Equation (13), electric field jump condition at the interface is governed by ˜s Ẽs · n = ˜f Ẽf · n on s (25) 3.1.2. Fluid equations. In this paper, we consider an incompressible viscous fluid described by the Navier–Stokes equations in \s , so the governing equations are the continuity equation and the momentum equation ∇ · vf = 0 (26)
f v̇f = ∇ · rf +
f g + ∇ · rM (27) where vf is the fluid velocity,
f is the fluid density, rf is the fluid Cauchy stress tensor, g is the gravitational force, and rM is the time-averaged MST. Brownian motion is not considered here because it has little influence on the motion of particles close to the electrodes (see the order of magnitude analysis section) (Figure 1). For a Newtonian fluid of uniform viscosity, rf can be written as rf = − p f I + [∇vf + (∇vf )T ] (28) where  is the fluid viscosity and p f is the fluid pressure. 3.1.3. Solid equations. The governing equations for a deformable solid s are
s v̇s = ∇ · rs +
s g + ∇ · rM (29) where vs is the velocity of the solid,
s is the solid density, rs is the solid stress tensor, and rM is the time-averaged Maxwell stress. For a hyper-elastic solid, we have, rs = − p s I + J −1 F(2(c1 + c2 Tr(C))I − 2c2 C)FT (30) where F is the deformation gradient, c1 and c2 are the material constants. The Cauchy–Green deformation tensor is defined as C = FT F and the Jacobian determinant is J = |F|. 3.2. Immersed formulation based on volumetric force and the Euler–Lagrange mapping In order to avoid remeshing in fluid–structure interaction problems, an extra force is introduced to enforce the interface traction and velocity conditions. This extra force, called the fluid–structure interaction force, will be added to the RHS of Equation (27) to extend the calculation domain from \s to . Here, we adopt the IFEM formulation given by Liu et al. [19]. Rewriting the solid momentum equation (29) into the form of Equation (27), and by noting that vs = vf in s , we have
s v̇s = ∇ · rs +
s g + ∇ · rM + (
s −
f ) (v̇s − g) + ∇ · (rf − rs ) + FFSI Copyright q 2006 John Wiley & Sons, Ltd. in s (31) Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 387 FFSI is the interacting force exerted on the fluid by the solid. In the fluid domain \s , we know that FFSI must vanish. Thus, the interacting force FFSI can be defined as the following:  −(
s −
f ) (v̇s − g) + ∇ · (rs − rf ) in s FSI F = (32) 0 in \s Based on Equation (32), Equation (27) can be extended to the whole domain  as
f v̇ = ∇ · rf +
f g + ∇ · rM + FFSI (33) Since the fluid and solid meshes do not generally coincide, in order to map the variables between the solid and fluid mesh, an Eulerian to Lagrangian mapping is introduced as ∀X ∈ s0 (X, t) ≡ (x(X, t), t) (34) This mapping function can be the standard finite element shape function or the Dirac-delta function approximated by the reproducing kernel particle function [28, 29] (RKPM). With this mapping function, we obtain the solid velocity and displacement from the Eulerian space. When the ions are the major current carriers and solids are insulators, we solve the following set of extended equations in the immersed formulation: ∇ · (f ∇) + S = 0 in  (35) where = S=  s in s f in f ⎧ ⎨ ∇ · ((s − f )∇s ) in s in f ⎩ −
e (36) (37) If the inherent carriers are the major current carriers, we can obtain the immersed formulation for the electric potential by extending Equation (23) into the whole domain  ˜ + Q =0 ∇ · (˜f ∇ ) (38) where ˜=  Q= ⎧ s ˜ ⎨ ⎩ ˜f   on s (39) on  \ s s ˜ ) ∇ · ((˜s − ˜f )∇  in s 0 in  \ s (40) Here, Q represents the influence of the difference in complex permittivities for fluid and solid on the electric field. The strong form of the immersed system is summarized in Table II. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 388 Y. LIU ET AL. Table II. Governing equations for immersed solid system. ˜ (x, t), ck (x, t)) Solution: (v(x, t), p(x, t), x(X, t),  Augmented fluid (in ) Continuity ∇ · v=0 (41) Momentum
f v̇ = ∇ · r +
f g + ∇ · rM + FFSI (42) Force FFSI (x, t) =  0, x ∈  \ s f s s (
−
)(v̇ − g) + ∇ · (r − r), x ∈ s (43) Stress r = − p I + [∇v + (∇v)T ] (44) Electric field and species equation when inherent carriers are major current carriers ˜) + Q =0 ∇ · (˜f ∇  Q=  in  (45) x ∈  \ s 0, (46) ˜ ), x ∈ s ∇ · ((˜s − ˜ f )∇  For ions as major current carriers case ∇ · (f ∇ ) + S = 0 in  S=  (47) ∇ · ((s − f )∇ ) in s
e = (48) in f −
e
ez k ck in f (49) k * ck + vf · ∇ck = ∇ · [−k ez k ck E + k kB T ∇ck ], *t k = 1, . . . , M in f (50) Solid equations (in s ) Solid stress rs = − p s I + *x *X −1 *x *x [2(c1 + c2 Tr(C)) I − 2c2 C] *X *X T (51) Maxwell stress rM = 14 Re(˜)((Ẽ ⊗ Ẽ∗ + Ẽ∗ ⊗ Ẽ) − |Ẽ|2 I) Copyright q 2006 John Wiley & Sons, Ltd. (52) Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 389 Interface (on s ) Normal stress r− · n = rs · n, r− ≡ r|\s (53) For inherent carriers as major current carriers case ˜s Ẽs · n = ˜f Ẽf · n (54) For ions as major current carriers case f E f · n = s E s · n (55) Initial values and boundary conditions v(x, 0) =  vf0 (x) ∀x ∈ f0 vs0 (xs ) ∀xs ∈ s0 xs (X, 0) = X ∀X ∈ s0  p0f (x) ∀x ∈ f0 p(x, 0) = p0s (x) ∀xs ∈ s0 v(x, t) = vfB (x) ∀x ∈ *  ck,0 (x) ∀x ∈ f0 ck (x, 0) = 0 ∀xs ∈ s0 (56) (57) (58) (59) (60) 3.3. Weak form For stabilization of the flow problem, the SUPG [30] (streamlined upwind/Petrov–Galerkin) formulation will be adopted. With stabilization parameters c and m , we take test functions for fluid velocity, pressure, electric potential, and ion concentration as w̃ ∈ V ≡ {w + m v · ∇w + c ∇q | w ∈ [H01 ()]d , q ∈ H 1 ()} (61) q̃ ∈ M ≡ {q + c ∇ · w | q ∈ H 1 (), w ∈ [H01 ()]d } (62) m ∈ P ≡ L 2 () (63) sk ∈ S ≡ L 2 () (64) In the IFEM, we consider the entire domain, including both the solid and fluid, and construct a single weak form (details and derivations are given in Liu et al. [18]). The weak form of the fluid momentum equation is   Rm (v, p, u, ) ≡
f (v̇ − g) · w̃ dx + (r + rM ) : ∇w dx  −
e Copyright q   (∇r + ∇rM ) · (m v · ∇w + c ∇q) dx e 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 390 Y. LIU ET AL. −  s0 (
−
s ) *v −g *t  *x (r̄ − r̄s ) : ∇X w̃ + s *X 0  *x dX *X · w̃ −1  *x dX = 0 ∀w̃ ∈ V *X (65) Here, r̄, r̄s , v̄, and w̃ are the mapped variables through the Eulerian to Lagrangian mapping defined in Equation (34). Individual finite elements are denoted with the subscript e and the stabilization supplementary terms are evaluated at the element interiors. The weak form of the continuity equation is  Rc (v, p, u) ≡ ∇ · vq̃ dx = 0 ∀q̃ ∈ M (66)  When the ions are the major current carriers, the weak forms of the electric potential and ions concentration equations are   Re () ≡ f ∇ · ∇mdx −
f m dx       −1 −1 *x *x *x s f dX=0 ∀m∈P (67) · (∇X m) + ( − ) (∇X ) s *X *X *X 0  Rion (ck ) ≡ [−k ez k ck E + k kB T ∇ck ] · ∇sk dx  −   * ck + v · ∇ck sk dx + *t  s0 ck sk *x dX = 0 *X ∀sk ∈ S (68) where we use a penalty to apply the constraint ck = 0 in s0 . When the inherent carrier is the major current carrier, the weak form of the electric potential equation is  ˜ · ∇m dx ˜ ≡ ˜f ∇  Re ()       −1 −1 *x *x *x s f ˜ dX = 0 ∀m ∈ P · (∇X m) (˜ − ˜ ) (∇X ) + *X *X *X s0 (69) 3.4. Finite element discretization There are two sets of discretization, one for the Lagrangian solid mesh and another for the Eulerian fluid mesh. Using the finite element shape functions, trial functions for v, p, m and c are as follows: In the fluid domain N N

vh = p J (t) N J (x) (70) v J (t)N J (x), ph = J =1 J =1 h = N
m J (t)N J (x), J =1 Copyright q 2006 John Wiley & Sons, Ltd. ch = N
c J (t)N J (x) (71) J =1 Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 391 In the solid domain vh = Ns
v̄ J (t) N̄ J (X), ph = J =1 h= Ns
Ns
p̄ J (t) N̄ J (X) (72) c J (t) N̄ J (X) (73) J =1 m J (t) N̄ J (X), ch = J =1 Ns
J =1 Here, N J (x) and N̄ J (X) are chosen as the linear FE-basis functions for fluid, and solid, respectively, on the independent meshes for the background flow and the initial solid. For any solid node in s , we note v̄ K (t) = N
v J (t) N J (x(X K , t)) ∀K = 1, . . . , N s (74) J =1 Consequently, the following finite element approximation of the Eulerian velocity transformed by Euler–Lagrange mapping is obtained:  s  N N

vh = v J (t) N J (x(X K , t)) N̄ K (X) (75) J =1 K =1 The linearity of the mapping function suggests the following approximation: ( N̄ J )(X, t) = Ns
N J (x(X K , t)) N̄ K (X) ∀J = 1, . . . , N (76) K =1 These basis functions will form a test function space for w in the weak formulation. These will also be applied to the pressure p(x, t), electric potential , and the ion concentration c. It should be noted that the moving interface s is automatically updated through the Lagrangian solid mesh and captured by the Euler–Lagrange mapping. Compared with the level set method [31] that is widely used to track the interface, this mapping function does not involve any extra equations or unknowns and is convenient to use. 4. RESULTS 4.1. Assembly of rigid particles Understanding the motion of a rigid particle near an electrode in the presence of a uniform/non-uniform electric field is important for the assembly or trapping of micro-particles. There have been some attempts by Dimaki and Boggild [14] to model the motion of CNTs under an AC field. In that work, the viscous effect of the flow was represented by a friction term in the equations of motion and the DEP force on CNTs was estimated by the point dipole theory. However, as stated before, the point dipole theory has limited accuracy. Alternative methods have been developed to study the motion of rigid particles in fluid [11, 32–34]. Here, we use the IEFEM method to simulate the assembly process of rigid particles of various shapes and under different electric fields. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 392 Y. LIU ET AL. Figure 3. Attraction of a rigid spherical particle towards an electrode gap. E-field contour shows that the E-field close to the sphere is highly non-uniform: (a) mesh; and (b) electric field contour. 4.2. Trapping of a spherical particle The term DEP force was introduced in 1951 by Pohl [35]. Jones [36, 37] gave a detailed account of an approximate theory called the EDM theory to model the DEP force. The DEP force on a particle depends on the strength and the frequency of the local electric field, the particle geometry, and the electrical properties of the particle and medium. A widely used expression for the time-averaged DEP force on a spherical particle is given as [36]   ˜s − ˜f 3 f FDEP = 2a  Re s ∇|E|2 (77) f ˜ + 2˜ where a is the radius of the particle, ˜f is the complex permittivity of the fluid, ˜s is the complex permittivity of the particle, and E is the magnitude of the electric field. Equation (77) is derived based on the EDM method via linear dipole approximation. However, when the dielectric objects Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 393 Table III. Various properties of the fluid and solid used in the simulation of trapping of a spherical particle. Fluid 36 797 nodes 190 621 elements (120 × 80 × 20)
f = 1 g/cm3  = 0.01 g cm/s f = 0.0056 S/cm f = 200 Solid 809 nodes 3684 elements
s = 1 g/cm3 D = 20
m s = 0.0112 S/cm s = 1000 Figure 4. The calculated DEP force on a spherical particle from the effective dipole moment theory and our numerical calculation at various frequencies. The sphere is placed at: (a) 10
m above the electrode gap; and (b) 5
m above the electrode gap, respectively. are very close to the electrodes or are clustered together, the linear dipole approximation is no longer accurate. We next show that the MST is more accurate than the EDM theory. To validate the MST method used in the simulation, we study the attraction of a rigid spherical particle in the presence of an electric field formed by two collinear but separate electrodes as shown in Figure 3. A spherical particle of diameter 2
m is placed 10
m above an 8
m-wide gap between two parallel electrodes. The applied AC voltage is 4 V and the centre of the particle is at the midpoint of the gap. The geometry of this simulation is shown in Figure 3, and the properties of the fluid and solid are listed in Table III. The DEP force on the spherical particle for a wide range of frequencies calculated by both the EDM theory and our numerical method are plotted in Figure 4. The DEP force computed by the EDM and MST agree well at high frequencies (the difference is less than 2%), while the MST yields higher value at low frequencies (the difference is around 6%). The difference between the results of the MST and EDM indicates that EDM theory becomes invalid when the imaginary part of the complex permittivity is large. A comparison in closer proximity to the electrodes (e.g. 5
m-high position), however, shows large differences because of the large variation in the gradient of the electric field. The EDM model leads to a much larger force on the sphere than that obtained by the MST theory. Therefore, the use of EDM model to obtain DEP force may lead to incorrect Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 394 Y. LIU ET AL. Figure 5. Dielectrophoretic force time history on the particle. Table IV. Various properties of the fluid and solid used in the simulation of assembly of NWs. Fluid 42 677 nodes 226 532 elements (120 × 80 × 24)
f = 1 g/cm3  = 0.01 g cm/s f = 0.0056 S/cm f = 200 Solid 1536 nodes 6253 elements 2745 surface elements t = 1 ms
s = 1 g/cm3 Length = 6
m s = 0.0112 S/cm s = 1000 Diameter = 600 nm Caps radius = 300 nm assembly pattern of the particles near the electrodes. The DEP force history of the particle during the trapping process is shown in Figure 5. It can be seen that as the particle approaches the gap, the DEP force on the particle increases quadratically. 4.3. Assembly of rod-shaped particles As demonstrated in the previous section, the DEP calculation near the electrode by the EDM is inaccurate. Furthermore, the EDM theory is only accurate for spherical particles, whereas many assembly processes involve non-spherical particles such as NWs, viruses, and DNA (Table IV). Here, we study the assembly of two NWs between two semi-circular electrodes to demonstrate the capability of the MST formulation to model the DEP on non-spherical geometries. We consider a pair of semi-circular electrodes of radius R = 5.0
m with a minimum gap of 5
m, as shown in Figure 6. Our objective is to simulate the assembly of two NWs using DEP. Two NWs are randomly positioned and oriented at a height of 2
m above the electrode surface. Figure 6 shows the dynamic process of assembly: (1) rotation of the particle to align its long-axis along the electric field direction; (2) trapping of the particle towards the gap center; (3) final deposition on the bottom of the electrode. Various assembly patterns can be achieved for electrodes of different shapes or in different electric fields. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 395 Figure 6. Simulation of the assembly of two rod-shaped NWs by the application of an AC field. Two pairs of semi-circular-shaped electrodes are used here with a gap size of 5
m. The AC field has a frequency of 5 MHz and amplitude of 0.5 V/
m. Figure 7. Geometry and mesh used for electrodeformation simulation. 4.4. Electrodeformation of cell The transport and manipulation of cells in micro-fluidic channels are essential elements of the lab-on-a-chip bio-devices. Using such devices, Riske and Dimova studied the electrodeformation of giant vesicles [38] with a fast digital imaging technique. A DC electric field of 3 kV/cm, with a duration of tp = 200
s was applied to a giant vesicle with a diameter of 20
m. We study the same electrodeformation process of the experiment in [38] and compare our simulated results with the experiments. Specifically, the equilibrium cell shape, the assembly Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 396 Y. LIU ET AL. Table V. Various properties of the fluid and solid used in the simulation of electrodeformation. Fluid 24 046 nodes 124 556 elements (40 × 90 × 20)
f = 1 g/cm3  = 0.012 g cm/s f out = 6
S/cm f in = 4.5
S/cm Solid 3027 nodes 8989 elements t = 4
s for the first 250 steps
s = 1 g/cm3 Outer diameter = 20
m 0.5 ms for the next 100 steps E Y = 5.7 × 105 g/(cm s2 ) Inner diameter = 19.8
m 10 ms for the last 50 steps Figure 8. A snapshot sequence of a vesicle subjected to a pulse, E = 3 kV/cm, tp = 200
s. The top ones are experimental images from [38] (Riske KA and Dimova R, Electro-deformation and poration of giant vesicles viewed with high temporal resolution, 2005, Copyright Biophysical Society, Reproduced with permission) and the bottom ones are simulation snapshots. pattern of the cells, and the dynamic manipulation process are simulated. The properties of the fluid and solid used in this simulation are listed in Table V and the geometry is shown in Figure 7. The membrane is assumed to be 100 nm thick and is represented by two layers of tetrahedral elements. The modulus of the membrane is measured to be E Y = 5.7 × 105 g/(cm s2 ) [38]. The fluid inside the membrane is a 0.2 M (M is molar, 1 molar = 1 mole per litre) sucrose solution with a conductivity of 6
S/cm and a viscosity of 0.012 g cm/s. The fluid outside the membrane is a 0.2 M glucose solution with a conductivity of 4.5
S/cm and the same viscosity as the inner fluid. The dynamic process of stretching and relaxation is illustrated in Figure 8. The vesicle is found to be deformed into a prolate shape once the electric field is applied and relaxes gradually into the original spherical shape after the electric field is shut. The dynamic deformation process and the equilibrium pattern obtained through numerical simulation are compared with the experimental images in Figure 8, which shows good agreement. As a quantitative measure of the deformation of the cell, we compare the degree of deformation as the ratio a/b, where a and b are the semiaxes along and perpendicular to the field, respectively, to the experimental values. The time history of the ratio is plotted in Figure 9. It is found that the cell is gradually stretched into a prolate shape when the electric field is applied. After the electric field is removed, the cell relaxes slowly back to its original shape under the membrane stress. As can be seen in Figure 9, the results agree with the experimental observations. 4.5. Assembly of multiple cells When multiple cells are close to each other, interaction of cells under an electric field leads to various assembly patterns. Our numerical method can be used to explore the assembly patterns Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 397 Figure 9. The time history for the degree of deformation of the vesicle. The circular marks are experimental data reproduced from [38]. Figure 10. Alignment of cells under the electric field: (a) cells aligned parallel to the electric field direction will be attracted towards each other and will be deformed; and (b) the pearl-chain formation of multiple cells aligned between cylindrical electrodes. Experimental images are from [8] (Zimmermann U et al., Electromanipulation of mammalian cells: fundamentals and application, 2000, Copyright IEEE, Reproduced with permission). under different electric field conditions. As a demonstration example, two and three cells are placed randomly between two parallel electrodes. These cells deform against each other and form a stretched chain as shown in Figure 10(a). Another example is the alignment of cells between two cylindrical electrodes. Thirty-two cells are placed randomly between two parallel electrodes. The Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 398 Y. LIU ET AL. Table VI. Problem parameters for the simulation of multiple cell assembly. Fluid 73 713 nodes 398 058 elements (200 × 40 × 40)
f = 1 g/cm3  = 0.01 g cm/s f = 200 Solid 809 nodes per cell 3684 elements per cell
s = 1 g/cm3 D = 20
m s = 1000 Table VII. Various properties of the three viruses used in the simulation of virus detection. Size (nm) Permittivity (0 ) Conductivity (mS/m) Inovirus Influenza Bacteriophage P22 200 (length) 20 (diameter) 70 8 100 (diameter) 65 (length) 3 80 30 3 non-uniform electric field near the cylindrical electrode attracts the cells towards the electrodes in sequence and forms a chain aligned perpendicular to the electrode surface, as shown in Figure 10(b). The properties of the fluid and solid used in this simulation are listed in Table VI. 4.6. Virus detection Manipulation of individual viruses is important for the lab-on-a-chip disease diagnosis system, because a diagnostic lab sample usually contains several types of viruses that need to be analysed separately. Centrifugation has been widely used to separate viruses, but it is cumbersome. Electric field-induced forces have attracted considerable interest for the separation and selective detection of viruses since they can be applied and controlled easily. Dielectrophoresis has been used to isolate viral particles between electrodes [39]. Since the sign of dielectrophoresis depends on the dielectric properties of the particles and the medium and the electric field frequency, a specific type of viral particles can be selectively attracted by applying a frequency to induce a positive dielectrophoresis. The trapped viral particles can finally be identified by a fluorescent method. Suppose the sample solution contains the Inovirus, Influenza, and Bacteriophage P22, whose physical and electrical properties are given in Table VII. Due to the limited knowledge of the elastic and electrical properties of viruses, the entire virus is assumed to be a homogeneous rigid body. We emphasize that the modelling of the viruses can be easily modified in our proposed method if more detailed information regarding their properties becomes available. For a single virus of each kind, we obtain the DEP force versus frequency curve as shown in Figure 11. The DEP force on a virus is calculated through the MST described in Equation (17). An example of a strategy to filter out the three kinds of viruses is as follows. First, the frequency indicated by the solid arrow (∼5 MHz) in Figure 11 is chosen to filter out Influenza in the first array of electrodes in our proposed device since only Influenza particles will experience positive DEP force. Then, the frequency indicated by the dashed arrow (∼1 MHz) in Figure 11 is chosen to filter out Inovirus in the second array of electrodes. To demonstrate the capability of the modelling method, a preliminary simulation is performed for the deposition of three viruses on a pair of parallel-rectangular-shaped electrodes as shown in Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 399 Figure 11. The frequency-dependent polarization factor of three viruses. The polarization factor here is defined as Re{K f } from the expression FDEP = f Re{K f }∇|Ẽ|2 , where  is 32 the volume of the virus and f is the medium permittivity. The polarization factor given by the effective dipole theory for a spherical particle is Re{K f } = Re{(˜s − ˜f )/(˜s + 2˜f )}. The solid arrow indicates the frequency at which Influenza can be selectively deposited by positive dielectrophoresis while Inovirus and Influenza remain in the medium due to negative DEP force; the dashed arrow indicates the frequency at which Inovirus can be selective deposited out of the medium after Influenza is excluded. Figure 12. Simulation of the deposition of three viruses under an AC field. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 400 Y. LIU ET AL. Figure 13. Mechanisms of electro-osmosis flow. Figure 12. The three viruses are initially aligned at the same height of 0.2
m above the electrode surface. When an AC field of 5 MHz is applied, the Inovirus and the Bacteriophage are under negative DEP, while the Influenza is under positive DEP. Finally, only Influenza gets trapped exactly in the gap and the other two are washed away by the flow. The above simulation demonstrates the capability of the IEFEM method in modelling complex electrokinetic conditions. 4.7. CNT rotation due to electro-osmotic flow In this section, we demonstrate the application of the IEFEM formulation to the study of effects of electro-osmotic flow on rotation of CNTs immersed in ethanol suspension. The electro-osmotic flow is driven by the electrostatic force applied to the charged double layer on the surfaces of the electrodes. In the actual implementation, the electro-osmotic flow is treated as a slip boundary condition for the fluid. Also, since the Debye layer (the layer close to the wall where the velocity is varying) is only a few nanometers, only the steady velocity is taken into account v=− f 0E  (78) where the zeta potential 0 = 135 mV, viscosity of ethanol  = 1.2 mN s/m2 , and the permittivity of ethanol f = 70.8 nF/m. The charge relaxation time of a liquid is defined as = f /f where = 0.01 s for ethanol. An AC electro-osmotic flow is induced only when the applied electric field period is less than the relaxation time f < 1/(2 ). In such case, the sign of the charge of the electrodes and the double layer alternates according to the sign change of the potential. The charge experiences an electrostatic force tangential to the surface, which acts on the electric double layer and induces the AC electro-osmotic flow, as illustrated in Figure 13. The forces are maximal at the electrode edge and decrease rapidly with increasing distance from the edge; such non-uniform forces may form a vortex. In our simulation, an AC field of 0.5 V/
m and 100 Hz was applied on a pair of parallel rectangular-shaped electrodes with a 5-
m gap, as shown in Figure 14. A local electro-osmotic flow near the edges of electrodes is found to induce vortices that lead to the rotation of the NWs. Similar rotation of NWs has been observed in experiments shown in Figure 14(a). Note that Figure 14(b) reproduces vortex flow in the bottom view cross-section of the experimental snapshots shown in Figure 14(a). This phenomenon of vortex flow of CNT suspensions induces local circulative flows that transport NWs towards the electrode gap. 4.8. DNA stretching Stretching and orientation of a single DNA molecule on a surface is essential in the study of DNA sequencing [40], DNA–protein interactions [41], and DNA-based lab-on-a-chip devices [42]. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 401 Figure 14. Rotation of carbon nanotubes induced by electro-osmotic flow; black arrows show direction and magnitude of the electro-osmotic flow at the surface of the electrodes. Figure 15. Illustration of the experimental set-up of DNA stretching by an applied electric field. A DNA molecule is a coil in a solution and may be adsorbed onto a surface due to non-specific adhesion. Many attempts have been made to stretch a single DNA molecule into a straightline using electric fields, mechanical flow, and other techniques. For example, DNA has been successfully stretched by the DEP force in the presence of an applied AC electric field [9, 43]. However, the dynamic process and the underlying mechanisms of DNA stretching are still unclear. Here, the IEFEM method is employed to study the stretching process of a ‘toy’ DNA molecule in the presence of an applied electric field. In experiment performed by Lee et al. [5], a pair of semi-circular-shaped aluminium electrodes (gap width = 6
m) were microfabricated on an oxide-coated (5000 Å) silicon substrate. The experimental set-up is illustrated in Figure 15. -DNA dissolved in TE buffer Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 402 Y. LIU ET AL. Figure 16. The bead-spring model of DNA. Table VIII. Various properties of the fluid and DNA used in the simulation of electrodeformation. Fluid DNA 31 754 nodes 156 703 elements Kuhn length = 1.6 nm t = 0.05 ms
f = 1 g/cm3  = 0.012 g cm/s f = 6
S/cm f = 4.5
S/cm E Y = 120 kB T /nm2 Bead mass = 1500 amu (Tris–ethylenediamine–tetraacetic–acid buffer, 400 g/mL, 12.7 nM) was diluted to ∼10 nM using ultrapure water. A 150-
L drop of the diluted solution (including about 100 DNA molecules) was added onto the electrode gap of the silicon substrate that was heated to a temperature of 100◦ C. Note that the double strand DNA (dsDNA) is denatured to a single strand DNA (ssDNA) at this temperature. When an AC field (0.82 V/
m at 5 MHz) was applied, an aggregated form of DNA was observed in the vicinity of electrodes. On the other hand, under a composite field combining an AC with a DC field (0.82 V/
m at 5 MHz, DC/AC amplitude ratio = 0.39), a single DNA molecule was stretched across electrodes [Figure 17(b)]. The manipulation of a DNA molecule involves two processes: the attraction of the DNA from far field towards the gap and the stretching of the DNA between the two electrodes. The first process is induced by the DEP force, while the second process is controlled by the electro-osmotic flow. In this simulation, we focus on the stretching process. There are various atomistic-scale models and coarse-scale models for a DNA. For simulation of DNA dynamics in nano-fluidic devices, coarse-scale models are usually employed. The widely used coarse-scale models can be categorized into discrete models and continuum models [44]. The continuum models of DNA usually treat the DNA as a continuum rod [45] that can bend and twist. As a first attempt, a highly simplified discrete model called the bead-spring model [46, 47] is adopted here to model a DNA molecule (Figure 16). In this model, each segment of the DNA is represented by a bead and the connections between the segments are represented by springs. The Brownian dynamics simulation technique we employ for bead-spring chains is taken from Doyle et al. [48, 49] and we refer the readers to these papers for more details. The Kuhn length (the length of hypothetical segments that the chain can be considered as freely joined) of a ssDNA is 1.6 nm and the Young’s modulus of the spring is 120 kB T /nm2 [48, 49] (Table VIII). In the simulation, a curved ssDNA chain with a contour length of 6
m is placed in the middle of a 3
m gap between two semi-circular electrodes. Thirty beads and 29 springs are used in the simulation for a single ssDNA chain. The equilibrium length Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 403 Figure 17. The stretching of a DNA chain between two semi-circular-shaped electrodes: (a) simulation of the stretching process; and (b) atomic force microscope (AFM) image of a stretched ssDNA. of the springs are set to be 200 nm, covering 125 Kuhn lengths. A DC field of 4 V is applied on the two electrodes. The electro-osmotic flow is modelled as a slip boundary on the electrode surface as described in Equation (78). The curved ssDNA chain is stretched into a straight line under the electro-osmotic flow [Figure 17(a)]. The dynamic stretching process and the final stretched pattern from the experiments by Lee et al. [10] are illustrated in Figure 17(b). 5. CONCLUSION A three-dimensional immersed electrokinetic finite element method (IEFEM) has been proposed and applied to study the electric field-driven assembly of particles in a liquid, the electrodeformation of cells, and the manipulation of viruses and DNA. In IEFEM, independent solid meshes move in a fixed background field mesh that models the fluid and electric field, and hence, the mesh generation for the fluid and moving deformable structure interaction problems is greatly simplified. The general continuum electromechanics equations are solved concurrently with the Navier–Stokes fluid equations. The dielectrophoretic force acting on individual particles depends on the electrical properties of the particles and medium, the geometries of the microelectrodes and the particles, the electric field strength and frequency. The assembly of both spherical and nonspherical rigid particles between opposing electrodes has been demonstrated. The dynamic process of electrodeformation of a cell has been studied and compared with experimental measurements. The pearl-chain formation of multiple cells has also been illustrated. The IEFEM method presented here is a powerful modelling and simulation tool for developing a fundamental understanding of electrokinetic phenomena. ACKNOWLEDGEMENTS The support of this research by the National Science Foundation (NSF), Office of Naval Research (ONR), Army Research Office (ARO), the NSF Summer Institute on Nano Mechanics and Materials, and NSF IGERT is gratefully acknowledged. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme 404 Y. LIU ET AL. REFERENCES 1. Rao SG, Huang L, Setyawan W, Hong SH. Large-scale assembly of carbon nanotubes. Nature 2003; 425:36–37. 2. Smith BW, Benes Z, Luzzi DE, Fischer JE, Walters DA, Casavant MJ, Schmidt J, Smalley RE. Structural anisotropy of magnetically aligned single wall carbon nanotube films. Applied Physics Letters 2000; 77(5):663–665. 3. van der Zande BMI, Koper GJM, Lekkerkerker HNW. Alignment of rod-shaped gold particles by electric fields. Journal of Physical Chemistry B 1999; 103(28):5754–5760. 4. Huang Y, Duan XF, Wei QQ, Lieber CM. Directed assembly of one-dimensional nanostructures into functional networks. Science 2001; 291(5504):630–633. 5. Lee KH, Chung JH, Lee JH. Superimposed AC- and DC electric field guided deposition of a single DNA molecule along a microfabricated gap. IEEE NANO, San Francisco, CA, August 2003. 6. Trau M, Saville DA, Aksay IA. Assembly of colloidal crystals at electrode interfaces. Langmuir 1997; 13(24): 6375–6381. 7. Brisson V, Tilton RD. Self-assembly and two-dimensional patterning of cell arrays by electrophoretic deposition. Biotechnology and Bioengineering 2002; 77(3):290–295. 8. Zimmermann U, Friedrich U, Mussauer H, Gessner P, Hamel K, Sukhoruhov V. Electromanipulation of mammalian cells: fundamentals and application. IEEE Transactions on Plasma Science 2000; 28(1):72–82. 9. Washizu M, Nikaido Y, Kurosawa O, Kabata H. Stretching yeast chromosomes using electroosmotic flow. Journal of Electrostatics 2003; 57(3–4):395–405. 10. Chung JH, Lee KH, Lee JH, Ruoff RS. Toward large-scale integration of carbon nanotubes. Langmuir 2004; 20(8):3011–3017. 11. Patankar NA, Hu HH. Numerical simulation of electroosmotic flow. Analytical Chemistry 1998; 70:1870–1881. 12. Legay A, Chessa J, Belytschko T. An Eulerian–Lagrangian method for fluid–structure interaction based on level sets. Computer Methods in Applied Mechanics and Engineering 2006; 195:2070–2087. 13. Solomentsev Y, Bohmer M, Anderson JL. Particle clustering and pattern formation during electrophoretic deposition: a hydrodynamic model. Langmuir 1997; 13(23):6058–6068. 14. Dimaki M, Boggild P. Dielectrophoresis of carbon nanotubes using microelectrodes: a numerical study. Nanotechnology 2004; 15(8):1095–1102. 15. Liu Y, Zhang L, Wang X, Liu WK. Coupling of Navier–Stokes equations with protein molecular dynamics and its application to hemodynamics. International Journal for Numerical Methods in Fluids 2004; 46:1237–1252. 16. Liu YL, Chung JH, Liu WK, Ruoff RS. Dielectrophoretic assembly of nanowires. Journal of Physical Chemistry B 2006; 110(29):14098–14106. 17. Wang X, Liu WK. Extended immersed boundary method using FEM and RKPM. Computer Methods in Applied Mechanics and Engineering 2004; 193(12–14):1305–1321. 18. Zhang L, Gerstenberger A, Wang X, Liu WK. Immersed finite element method. Computer Methods in Applied Mechanics and Engineering 2004; 193(21–22):2051–2067. 19. Kim DW, Tang SQ, Liu WK. Mathematical foundations of the immersed finite element method. Computational Mechanics. DOI 10.1007/s00466-005-0018-5, URL http://dx.doi.org/10.1007/s00466-005-0018-5:1–12, 2006. 20. Liu WK, Liu YL, Farrell D, Zhang L, Wang XS, Fukui Y, Patankar N, Zhang YJ, Bajaj C, Lee J, Hong JH, Chen XY, Hsu HY. Immersed finite element method and its applications to biological systems. Computer Methods in Applied Mechanics and Engineering 2006; 195(13–16):1722–1749. 21. Liu WK, Karpov EG, Zhang S, Park HS. An introduction to computational nanomechanics and materials. Computer Methods in Applied Mechanics and Engineering 2004; 193(17–20):1529–1578. 22. Liu YL, Liu WK. Rheology of red blood cell aggregation by computer simulation. Journal of Computational Physics 2006; 220(1):139–154. 23. Liu WK, Karpov EG, Park HS. Nano Mechanics and Materials: Theory, Multiple Scale Analysis, and Applications. Springer: Berlin, 2005. 24. Wang X, Wang XB, Cascoyne PRC. General expressions for dielectrophoretic force and electrostational torque derived using the Maxwell stress tensor. Journal of Electrostatics 1997; 39:277–295. 25. Green NG, Morgan H. AC Electrokinetics—Colloids and Nanoparticles. Microtechnologies and Microsystems. Research Studies Press, 2003. 26. Stratton JA. Electromagnetic Theory. McGraw-Hill: New York, 1941. 27. Saville DA. Electrohydrodynamics: the Taylor–Melcher leaky dielectric model. Annual Review of Fluid Mechanics 1997; 29:27–64. 28. Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids 1995; 20(8–9):1081–1106. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme IMMERSED ELECTROKINETIC FINITE ELEMENT METHOD 405 29. Li SF, Liu WK. Meshfree Particle Methods. Springer: Berlin, 2004. 30. Hughes TJR, Mallet M. A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective–diffusive systems. Computer Methods in Applied Mechanics and Engineering 1986; 58:305–328. 31. Quecedo M, Pastor M. Application of the level set method to the finite element solution of two-phase flows. International Journal for Numerical Methods in Engineering 2001; 50:645–663. 32. Patankar NA, Sharma N. A fast projection scheme for the direct numerical simulation of rigid particulate flows. Communications in Numerical Methods in Engineering 2005; 21(8):419–432. 33. Patankar NA, Joseph DD. Lagrangian numerical simulation of particulate flows. International Journal of Multiphase Flow 2001; 27(10):1685–1706. 34. Patankar NA, Joseph DD. Modeling and numerical simulation of particulate flows by the Eulerian–Lagrangian approach. International Journal of Multiphase Flow 2001; 27(10):1659–1684. 35. Pohl HA. Dielectrophoresis. Cambridge University Press: Cambridge, 1978. 36. Jones TB. Electromechanics of Particles. Cambridge University Press: New York, 1995. 37. Jones TB. Basic theory of dielectrophoresis and electrorotation. IEEE Engineering in Medicine and Biology Magazine 2003; 22(6):33–42. 38. Riske KA, Dimova R. Electro-deformation and poration of giant vesicles viewed with high temporal resolution. Biophysical Journal 2005; 88(2):1143–1155. 39. Morgan H, Green NG. Dielectrophoretic manipulation of rod-shaped viral particles. Journal of Electrostatics 1997; 42(3):279–293. 40. Woolley AT, Guillemette C, Cheung CL, Housman DE, Lieber CM. Direct haplotyping of kilobase-size DNA using carbon nanotube probes. Nature Biotechnology 2000; 18(7):760–763. 41. Kabata H, Kurosawa O, Arai I, Washizu M, Margarson SA, Glass RE, Shimamoto N. Visualization of single molecules of RNA-polymerase sliding along DNA. Science 1993; 262(5139):1561–1563. 42. Porath D, Bezryadin A, de Vries S, Dekker C. Direct measurement of electrical transport through DNA molecules. Nature 2000; 403(6770):635–638. 283RD Times Cited:532 Cited References Count:22. 43. Washizu M, Kurosawa O, Arai I, Suzuki S, Shimamoto N. Applications of electrostatic stretch-and-positioning of DNA. IEEE Transactions on Industry Applications 1995; 31(3):447–456. 44. Olson WK. Simulating DNA at low resolution. Current Opinion in Structural Biology 1996; 6:242–256. 45. Manning RS, Maddocks JH. A continuum rod model of sequence-dependent DNA structure. The Journal of Chemical Physics 1996; 105:5626–5646. 46. Levitt M. Computer-simulation of DNA double-helix dynamics. Cold Spring Harbor Symposia on Quantitative Biology 1982; 47:251–262. 47. Olson WK, Zhang PS. Computer-simulation of DNA supercoiling. Methods in Enzymology 1991; 203:403–432. 48. Doyle PS, Shaqfeh ESG, Gast AP. Dynamic simulation of freely-draining, flexible polymers in steady linear flows. Journal of Fluid Mechanics 1997; 334:251–291. 49. Doyle PS, Shaqfeh ESG, Gast AP. Rheology of ‘wet’ polymer brushes via Brownian dynamics simulation: steady vs oscillatory shear. Physical Review Letters 1997; 78:1182–1185. Copyright q 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2007; 71:379–405 DOI: 10.1002/nme