Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Electrokinetic Effects on Heterogeneous Ionic Transport Phenomenon based on
NonEquilibrium Molecular Dynamics Method in Microchannel Flow
Yanling Wang, Yueshe Wang, Hao Yu and Qin Yu
Xi'an Jiaotong University, State Key laboratory of Multiphase Flow in Power Engineering,
28' Xianning west road, Xi'an, 710049, China
Email:wangys@mail.xjtu.edu.cn
Keywords: molecular dynamics, electrokinetic effects, external sinusoidal force, LJ potential
Abstract
Electrokinetic effects have a great impact on the flow characteristics of the solution streaming in microchannels. The
influences of the ionic particles in the solution which result from electrokinetic effects bring instabilities of the flow behavior
in it. In this paper, a NonEquilibrium molecular dynamics (NEMD) method is used to simulate flow characteristic of the
heterogeneous multi ionic particles, a concentration of 2.08mol/L NaCl solution, in a 5.84x5.84x5.84 nm cubic microchannel.
The solution is composed of Na+, Cl and water molecules. In order to illustrate the interactive force between any two
individual particles, LennardJones (LJ) potential energy, Coulomb force and the iondipole interaction between any two
particles are taken into account. A periodic disturbed force model is proposed, in which an external perturbing sinusoidal force
has been applied on the flow passing through the microchannel. The simulation results show the existence of the resistance
arose from the electrokinetic effects, which is in action mainly near the channel walls but ignorable in the bulk region.
Additionally, the velocity profile predicted by molecular dynamics also conforms to the continuum predictions of the
analytical solution based on the NavierStokes equations especially in the bulk region. All of these results provide a probable
reason why the flow characteristic through a microchannel is particularly different from that in a macro duct.
Introduction
The liquid flow hydrodynamics in microchannels has
gained increasing attention in recent years. This micro
fluidics has a wide range of applications in many fields,
such as microelectromechanical systems (MEMS),biochip,
labonachip, microrobot and micro air vehicle (MAV)
Daejoong Kim (2006). Fundamental understanding of
liquid flow in micro channels is critical to the design and
process control of various MEMS and modern onchip
instruments used in chemical analysis and biomedical
diagnostics. Since the surfacetovolume in micro channels
is larger than that in the conventional channels with the
size reduction, the flows in micro channels are strongly
affected by the channels surface properties. So many
phenomena of liquid flow in microchannels, such as
unusually high flow resistance, cannot be explained by the
conventional theories of fluid mechanics.
One of the significant influences of such interfacial
phenomena is electroviscous effects in the microscale,
which is induced by the interaction between the charged
channel walls and the charged particles in the fluids. To
understand the fluidic phenomenon, many experimental
and computational research efforts have been taken.
Experimentally, Shankar Devasenathipathy & Juan G.
Santiago (2002) had applied particle tracking techniques
to obtain spatially resolved velocity measurements in
electrokinetic flow devices, and measured the bulk
velocity field for flow in the region of a borosilicate
rectangular capillary (40 by 400 y m) using current
monitoring. Gh.Mohiuddin Mala & Dongqing Li
(1999)investigated the water flowing through microtubes
with radius ranging from 50 to 254 p m, the experimental
results are roughly agreement with the conventional theory.
For a low Re, the required pressure drop is approximately
the same as predicted by the poiseuille flow velocity, but,
as Re increases, there is a significant increase in the
pressure gradient compared to that predicted by poiseuille
flow theory. On the other hand, various types of simulation
techniques have also been used, ranging from continuum
models and atomistic simulation. Jonathan B.Freund used
an atomistic simulation of an electroosmotic flow in a
65.3Awide channel to study its physical details and
evaluate continuum models. He found that the ion
distribution in the atomistic simulation was in general
agreement with the PoissonBoltzmann theory away from
the wall, but ions were more attracted to the wall than this
theory predicted. A.S.Ziarani and A.A.Mohamad (2005)
had explored a nonequilibrium molecular dynamics
(NEMD) simulation and analytical solution of Poiseuille
flow through a nanochannel. The argon fluid was studied
with a total number of 2000 particles.
As clearly reported in some review papers(A.S.Ziarani, et
al. 2005), fluid flow in micro channels has been
Paper No
extensively studied, while, on the other hand, only a
limited number of investigations deal with electrokinetic
flow in micro channels. To the authors' knowledge, few
data are available for the flow velocity and the distribution
characteristic of different charged particles in fluid flow in
a microchannel with an external perturbing force field.
The present paper will report the effect of an external
perturbation, sinusoidal force, on the flow passing through
a microchannel. In the follow section, the details and
results of the NEMD simulation will be presented. The
analytical solution of the NS equation for the velocity
profile will be also discussed.
Nomenclature
Re Reynolds number
n Concentration of fluid (mol.L1)
z, The valence of typei ions
e unit charge (1.60x1019C)
kB Boltzmann constant (J/K)
T temperature (K)
Pe Volume charge density.(number/m3)
U Potential energy (V)
r, Position of molecule i relative toj (m)
re The cutoff radius (m)
q Quantity of charge (C)
1 Distance between the centre of the ion and the
centre of the water molecule (m)
u Velocity of particles (m/s)
p Fluid density (kg/m3)
Fo External applied force on particles (N)
At Time steps (s)
Greek letters
0 Electric potential (V)
e Dielectric constant (J/K)
o Characteristic length (m)
i Shearviscosity (Pa.s)
Mathematical model and external force
I Electrokinetic theory
It is wellknown that most solid surfaces carry electrostatic
charges as well as the liquid contains a certain amount of
ions because there must be some impurities in the fluid.
The interaction between the charged channel walls and
ions in the fluid will have a great influence on the flow
characteristic in the microchannel flow.
Assuming equilibrium Boltzmann distribution is applicable
to describe the distribution of the ions, which coincides
with the derivation based on the statistical mechanics.
ze y/
n, = n,, exp( Z ) (1)
kBT
Where n,o and z, are the bulk ionic concentration and the
valence of typei ions, respectively, 0 is the electrical
potential, e is the electron charge (i.e., 1.6x1019C). kB is
the Boltzmann constant, and T is the absolute temperature .
According to the theory of electrostatics, the relationship
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
between the electrical potential and the local net charge
density per unit volume at any point in the solution is
described by the Poisson equation:
V2= P (2)
Where E is the dielectric constant of the solution. pe is the
volume charge density.
For symmetric cations and anions
Pe = ze(n n)
ze y
2zeno sinh( )
kBT
Substituting eq.(3) into the Poisson equation (2) leads to
the wellknown PoissonBoltzmann equation:
2zen, ze i/
S 2zeno sinh(e V) (4)
E kT
Generally, solving this equation with appropriate boundary
conditions, we can obtain the electrical potential
distribution ,and the local charge density distribution Pe
can then be determined from Eq.(3).
II Model in rectangular microchannel
In this context, we use NEMD method to simulate the
distribution of all the flowing particles in a rectangular
microchannel. Figure 1 shows a schematic diagram of the
system under investigation. The model is a cube of
5.84x5.84x5.84 nm. The two channel walls are
symmetrical with respect to the channel center line. Each
wall is made up with two layers of silicon atoms. The
channel height is 5.256 nm. The sodium chloride dilute
solution of 2.08mol/L is used to simulation the charged
particles flow in the rectangular channel. We studied two
different cases in this simulation. The number and the
radius of the particles are listed in the table 1. There are
6600 particles in all of the model systems.
In this paper, we assume that the positive charges are
discrete in case and uniform distribution among the wall
atoms in case2, and the total charge quality of the fluid is
neutral in case. However, there is only Cl in the water
molecules in case2, and the number of charges in which
are equal to the charges on the walls. The innermost wall
layers take positive charges, corresponding to surface
charge of 0.938C/m2. All of the molecules are spherical
and rigid.
All particles of our microchannel flow model interact via
the LennardJones potential
U/12 1 6 12 !6
U,(r)= 4 r1 + ] (5)
Here U ind r, are the potential and the distance between
particle i and particle j, respectively. e is a parameter
characterizing the strength of interactions, and a is the
characteristic length.
In addition, there are coulomb interactions between
charged particles. The potential is
U, (r) q (6)
S 4xe r
Whereq, q, is the charges magnitude of the ion i and ionj
measured in units of e.
Moreover, while the water molecule is polar molecule, the
Paper No
charges distribution is asymmetrical though the total
charges are zero. The iondipole interaction between
charged particles and water molecular is not neglectable.
z
x_ X
Figure 1: A
investigation
schematic of the channel system under
0.3 0.35 0.4 0.45
distance between two ions a
Figure 2: LJ potential energy between different particles
The iondipole potential energy between charged particles
and water molecules is
U, ) qz,el cos 0
U, (r)= 4zr r2 (7)
47TS r
Where q is the charges magnitude of the ion, E o is
vacuum dielectric constant, z, is the ionic charges of water
molecule measured in units of e, I is the distance between
the centre of the ion and the centre of the water molecule,
and 0 is the angle between the line links the position of
ion and the mass center of the water molecule and the line
links the anion and the cation of the water molecule.
Generally, u =ql is called as the dipole distance.
Table 1 : Summary of the simulations performed
Si Cl Na+ H20 Total
radius (nm) 0.146 0.181 0.102 0.14
Particle
Casel number 1600 250 250 4500 6600
Particle
Case2 number 1600 150 4850 6600
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
performed to remove the local contacts. To start the
simulation, an initial velocity sample from a Maxwell
distribution at 300 K is assigned to each molecule in the
system. The particles are bounded by horizontal walls in Z
direction and periodic boundary conditions are imposed in
the other two direction. We use velocity Verlet algorithm to
integrate the equation of motion of the interacting particles
and follow their trajectory.
Table 2 : Parameters for the LennardJones
potential
interaction e (kJ/mol) a (nm)
00 0.6502 0.317
0Na 0.5216 0.286
OC1 0.5216 0.375
0Si 0.3278 0.327
Na+ Na 0.4184 0.257
Na Cl 0.4184 0.339
Na+Si 1.0118 0.295
C1l C1 0.4184 0.445
ClSi 1.0118 0.388
III. External perturbing force field
The fluid molecules have been assured to be constitutes of
quantities of spherical molecules. A sinusoidal force has
been applied to the system in the X direction, which is
Fe=Fo[l+sin( ct)] In this work, the amplitudes of the
applied force range from 1.17 x10 1Nto 1.17x10ON, and
the <=50 in the two cases. The continuum analytical
solution of the sinusoidal forced flow in a channel has been
derived. The governing impressible NS equation (White,
2003) with the aforementioned form of force function can
be written in the following form:
au LI, a2
a az2 + F [1 + sin( ot)] (8)
3t p z2
Where p is the fluid density, u velocity in the X direction,
and i is shear viscosity.
*1
11 l J 1 ilJ
Az
Figure 3 : Node of the difference method
To solve the equation (8), we divide the channel height and
simulation time uniformly into m and n parts, and use the
finite difference method to discrete the equation (8) in the
following form:
When setting up the simulation, the flow particles are
positioned uniformly. An energy minimization is
Paper No
At
1u,)= 
pAz
2u, +u,,)
+ Fo [1+ sin(o x At x j)]
Uj+1 = ru_1,j + (1 2r)ul + ru +l,j
Eq. (10)
E + AtF [+ sin( x At x j)(10)
Wherei=1,2,, m1andj=0,1, n1.
Numerical results and analysis
The varieties of energy profile are presented in Fig.4 and
Fig.5. Fig.4 is the simulation results about energy varieties
without the action of an external force, while the Fig.5
shows the energy varieties with the action of a sinusoidal
force. According to Fig.4, the total energy is roughly
equilibrium between kinetic energy and potential energy
when there is no action of external forces. It should be
mentioned that these profiles are obtained by statistical
averaging of each particles.
Ix 1021
0 0.5 1 1.5 2
time steps x104
Figure 4 : Varieties of the energy without external force
x1021
1000 2000 3000 4000 5000
Timesteps
Varieties of the energy with sinusoidal force
Fo=1.17xl10N.
Fig.5 shows the kinetic energy profiles under the action of
a sinusoidal force. The kinetic energy plays an important
role in the varieties of the total energy which is affected by
the force amplitude, as expected. The total energy varies
with the changes of the external force within a certain
range.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The distribution of all the particles flow in the
microchannel with the NEMD simulation is shown in the
Fig.6. An external sinusoidal force Fo=1.17x 10N is used
to simulated the characteristic of the particles flow. The Cl
number density as a function of distance from the wall is
shown in Fig.6. As expected, the profile is sharply peaked
near the walls, and fall to a low value by the middle of the
channel. This is mainly because the coulomb attraction
interactions between the charged walls and the Cl.The
small bumps near the peak at both walls result from
molecular stacking. They are roughly two water molecules
width away from the peak. There is a peak value of the
number density of Na+ near the channel wall as well,
however, the peaks do not cling to the channel walls and
there is an obviously decline about the number density of
Na+ between the peak and the channel wall as there is a
coulomb repulsion interactions between the charged walls
and the Na+. The total particles number density is also
depicted in Fig.6. According to this graph, the number
density profiles in the bulk region reach a fairly constant
profile. The number density is affected by the action of
electrokinetic effects near the channels, but the effects of
electrokinetic are ignorable in bulk region.
2.5
2
S1.5
0
I 1
E
c
0.5
total particles
 Na+
Cl
0 0.2 0.4 0.6 0.8 1
position (z)
Figure 6 : Number density of Na+, Cl and total particles
with sinusoidal force for case 1, Fo=1.17 x101N.
11i
0 08
a
) 07
.o
5)
mean velocity
 Na+ velocity
 Cl velocity
04
0 01
02 03 04 05 06 07 08 09
position (z)
Figure 7 : Related velocity versus channel depth for
sinusoidal force for case 1, Fo= 1.17 x10N.
The simulation results of the mean velocity profile of fluid
for the amplitude of 1.17 x1011 for the sinusoidal force is
shown in Fig.7. We see that they are approximately
parabola in the middle of the microchannel, but do not
continue as a parabola all the way as one would expect for
6
4
2
0
2 kinetic energy
potential energy
4 total energy
6
kinrtic energy
 potential energy
total energy
Figure 5 :
m
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
constant g Poiseuille flow. Instead, there is considerably
more resistance close to the wall, and the velocity drops
continuously toward the wall. Fig.7 also shows the
varieties of velocity about Na+ and Cl versus the depth of
the channel. We can see that either the velocity of Na+ or
Cl drops near the wall. However, the mean velocities of
Na+ and Cl are both lower than the total mean velocity of
the fluid. Besides, there are two velocity peaks of Cl
towards to the walls, and the peak value is about 20%
higher than the mean velocity of C1. In interpreting these
results we should note that the strength of the external
force field necessary to force a flow velocity large enough
to converges velocity statistics might preclude the
existence of the electrokinetic. Furthermore, we can see
that there is a velocity slip near the channel walls.
0.4 0.(
position (z)
Figure 8 : Comparison of velocity versus channel depth for
sinusoidal force for case 1, F=5.85 x1011N.
Fig. 8 shows the comparison of analytical solution of the
equation (8) and molecular simulation result for the
velocity profile of the sinusoidal force flow. The sinusoidal
force is Fo=5.85x10"11N. The result investigated by
B.Freund about electroosmosis in a nanometerscale
channel with an external electric field is also shown in Fig.
8, and the applied electric field acted on each particle in
the flow direction is 2.42x10"N in B.Freund's work. This
force is large but necessary to drive the flow at high
enough speed to converge velocity statistics.
It is clear that the overall velocity profile produced by
NEMD shows good agreement with the continuum
model's predictions of parabolic behavior and the result of
B.Freund beyond the verynear wall region. However, the
velocity of the result of NEMD simulation is much larger
than the other two results near the walls, which means that,
for a channel with a height of 5.256 nm, the action of the
external force plays an important role in the flow
characteristics. With the increase of the external forces, the
continuum hypothesis is invalid and there is a velocity slip
near the channel walls.
The mean number density for case 1 under the external
sinusoidal force Fo=5.85x1011N is depicted in Fig.9. The
graph also shows the result simulated by B.Freund who
investigated the fluid flow in a channel under an external
electric field. The mean density simulated by NEMD
methods agrees well with the result simulated by B.Freund.
But the NEMD result has a larger fluctuation than the
result of B.Freund. Clearly the magnitude of the external
force is very important.
200
 result of MD
 result of B Freund
0 0.2 0.4 0.6 0.
position (z)
Figure 9 : Comparison of mean number
channel depth for sinusoidal force
Fo=5.85x10"N.
density versus
for case 1,
B.Freund compared electroosmotic flows for uniformly
charged walls as well as for discretely charged walls and
reported no significant difference for the two cases. In this
work, we also simulated uniformly charged walls in case 1
as well as for discretely charged walls in case2. The fluid
is electric neutral as the number of Na+ is equal to the
number of Cl in case, and the charges distribution at the
walls is discrete. Whereas there is only Cl in the fluid in
case2, as a whole, the fluid is charged. The number of
charges in fluid is equal to the number of the charges on
the walls in case2, and the charges distribution on the walls
is uniformity.
NEMD simulation of particles distribution is displayed in
the Fig.10. The coulomb interaction between Cl and
charged walls is obvious in Fig. 10, and consequently there
is a peak of number density of Cl near to each wall. While
the mean velocity close to the walls is lower than that in
the bulk region. For the velocities of both the Cl and the
mean fluid are smooth during the bulk area.
0.8
S0.6 emean fluid number density b
  CI number density
0.4
0.2
0 0.2 0.4 0.6 0.8 1
position (z)
Figure 10: Number density versus channel depth for
sinusoidal force for case 2, Fo=5.85 x 10 N.
We finally study the NEMD simulation of fluid flow in
microchannel with a larger external force of
Fo=1.17xl010N. The results of the velocity distribution
and number density are shown in Fig.11 and Fig.12,
respectively. From the Fig. 11 it can be seen that the mean
Paper No
Paper No
fluid velocity is almost a fixed value across the channel
depth, while the velocity of Cl has a drop close to the
walls. This shows that the external force play an important
role in the varieties of fluid flow velocity, and the
electrokinetic effects are mainly in action near the
surfaces. It is also shown that the larger the external force
is, the greater the velocity slip will be along the wall
surface.
Fig.12 depicts the number density of the Cl and mean
fluid for case2. From this graph, we can see that the
number density of Cl verus channel depth has little
fluctuation under a larger external force Fo=1.17xl01'N.
This is because the electrokinetic effect is very little
related to the external force, and there is almost no obvious
difference about the forces applied on the Cl and water
molecule.
a
S0.6
0.4
 0.4
.me umainmumauamInuauauuI

d
 mean fluid velocity
 Cl velocity
0 0.2 0.4 0.6 0.8 1
position (z)
Figure 11 : Comparison of velocity versus channel depth
for sinusoidal force for case 2, Fo= 1.17 x10ON.
1.4
1.2
1
S0.8
0
S0.6
E
S0.4
0.2
0
e mean fluid number density
SC number density
I I
aI
0 0.2 0.4 0.6 0.8 1
position (z)
Figure 12 : Number density versus channel depth for
sinusoidal force for case 2, Fo=1.17 x101N.
Conclusions
In this content, a nonequilibrium molecular dynamics
simulation is used to study the characteristic of fluid flow
in microchannel. A total number of 6600 particles
including the fluid and the walls are investigated. The
simulation results show that the charged surfaces of the
walls have a great influence of the distribution of the flow
particles under the action of electrokinetic effects, which
will lead to high concentration of the particles near the
wall and complex the flow characteristic. At the same time,
an external force is applied to study the varieties of flow
velocity in the microchannel. The total flow characteristic
induced by electrokinetic effects is diminishing with the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
increasing of the external force. The microchannel flows
need to be further investigated, and the NEMD simulation
is one of the most important methods to study the flows
characteristic in microchannel.
Acknowledgements
The authors would like to acknowledge the financial
support of the National natural science of Foundation of
China (Grant No.50676077) and (No. 50823002)
References
Daejoong Kim & Eric Darve. Molecular dynamics
simulation of electroosmotic flows in rough wall
nanochannels. Physic Review E 73,051203(2006).
Shankar Devasenathipathy, Juan G. Santiago. Particle
tracking techniques for electrokinetic microchannel flows.
Analytical Chemistry, vol.74, 37043713, 2002
Gh.Mohiuddin Mala, Dongqing Li. Flow characteristics of
water in microtubes. International Journal of Heat and
Flow. Vol. 20, 142148, 1999
Jonathan B.Freund. Electroosmosis in a nanometerscale
channel studied by atomistic simulation. Journal of
Chemical of Physics. 116(5)2002.
A.S.Ziarani.and A.A.Mohamad, A molecular dynamics
study of perturbed poiseuille flow in a nanochannel,
Microfluid nanofluid, vol.2,pp. 1220,2005
White FM. Fluid mechanics, 5t edn. McGraw Hill, Boston,
Massachusetts.
Liu Jing. Micro/nanochannel scales heat transfer. Science
Press, Beijing, China, 2001
Michael R.Philpott, James N.Glosli, Shengbaizhu.
Molecular dynamics simulation of adsorption in electric
double layers .Surface Science 1995, 335:422431
R,Qiao,N.R.Aluru. Ion concentrations and velocity profiles
in nanochannel electroosmotic flows .Journal of Chemical
Physics 2002, 118:46924701
Wei Zhu, Sherwin J. Singer. Electroosmotic flow of a
model electrolyte. Physical Review,2005, E 71,041501
Dongqing Li.Electroviscous effects on pressuredriven
liquid flow in microchannels .Colloids and Surfaces
2001,195:3557
George Em Kamiadakis, Ali Beskok. Micro flows
fundamentals and simulation. springerverlag, New Y
ork.2006
Li yigui, Lu jiufang. Electrolyte solution theory.Tinghua
Press, Beijing, 2005
Feynman,R. Leighton,R, Sands,M .The Feynman Lectures
on Physics,1977
