Comput. Methods Appl. Mech. Engrg.

198 (2009) 3852–3868

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg.

journal homepage: www.elsevier.com/locate/cma

Lattice Boltzmann method for electrowetting modeling and simulation

H. Aminfar, M. Mohammadpourfard *
School of Mechanical Engineering, University of Tabriz, Tabriz, Iran

a r t i c l e i n f o a b s t r a c t

Article history: In this paper a free energy-based lattice Boltzmann (LB) approach for modeling and simulation of elec-
Received 11 December 2008 trowetting is presented. First, based on the solution of the Poisson–Boltzmann (PB) equation, it is shown
Received in revised form 28 June 2009 that the potential distribution inside a microsized droplet can always be assumed one-dimensional. Then,
Accepted 28 August 2009
spreading and transporting of a microsized fluid droplet on a flat substrate, due to the electrowetting, are
Available online 2 September 2009
studied. By applying a voltage to the substrate, liquid–solid surface tension is locally decreased. To esti-
mate the contact angle after applying the potential, new relations for surface tensions based on the free
Keywords:
energy functional Minimization are proposed. The results are compared to the analytical and experimen-
Lattice Boltzmann method
Free energy
tal results reported in the literature. A good agreement is found.
Poisson–Boltzmann equation Ó 2009 Elsevier B.V. All rights reserved.
Potential distribution
Electrowetting
Microsize droplet

1. Introduction CV 2
cos h ¼ cos hY þ ; ð1Þ
2rgl
Nowadays, the electrical control of wettability of liquid on a
dielectric solid, which is so-called electrowetting, has attracted where h is the droplet contact angle, hY is the contact angle with no
the attention of scientists. Electric forces, generated because of applied voltage, C is the capacity per unit surfaces, V is electric po-
interactions between electric field and surface charges, can lead tential, and rgl is liquid–gas surface tension. It is the fundamental
to modify the surface tension or surface free energy. In fact, by basis for the electrowetting. Considering Eq. (1), V equals
rffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffi
applying the voltage to the electrode which droplet has been
placed upon it, double-layer charge distributions occur in the li-
rgl drgl
V  ; ð2Þ
quid drop in close proximity to the solid surface, thereby the sur- C ee0
face energy of the liquid–solid interface is reduced [1]. here, d and e are the insulator thickness and permittivity, respec-
In connection with the Young equation which determines the tively. With the capacitance C being 100 pF/cm2 for electrowetting
wetting characteristics of the surface, this reduction in the surface on insulator coated electrodes (i.e., Fig. 2c) while in conventional
energy leads to a decrease in the equilibrium or Young’s contact investigations on bare metal surfaces (i.e., Fig. 2a), C is at least
angle (the angle at which the liquid drop joins the solid, Fig. 1a), 104 times larger [4]. Thus, instead of the 100 mV necessary to drive
causing the drop to further wet the surface, Fig. 1b. This has an a bare metal-electrolyte, the electrowetting on insulator coated
important role on the behavior of liquid layers such as spreading, electrodes needs hundreds of volts, which is not prohibitive as long
film instability, adhesion and subsequent spreading of biological as the electric field remains under the breakdown threshold in the
cells [2]. Three typical cases of the electrical effect on wetting have insulator.
been presented in Fig. 2, [3]. Recently, electrowetting has been proposed as a mechanism for
The above-mentioned modification of surface tension modifies transporting, mixing and dispensing droplets with volumes in the
the contact angle which can be estimated from the principle of range of nano-liters or less [5–11], ‘‘lab-on-a-chip” systems for
minimum energy, leading to the Lippmann–Young formula [1] applications such as DNA and protein analysis, and biomedical
diagnostics [12,13].
For the optimization of the electrowetting performance in the
microfluidic systems, it is necessary to understand the basic fluidic
issues occurring on the microscales. In such cases, computer simu-
* Corresponding author. Tel.: +98 411 3392461; fax: +98 411 3354153.
E-mail address: Mohammadpour@tabrizu.ac.ir (M. Mohammadpourfard). lations can be useful to understand the driving forces leading to the

0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.
doi:10.1016/j.cma.2009.08.021
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3853

Fig. 1. Schematic diagrams of the contact angle: (a) before applying voltage, hY and (b) after applying voltage, h.

More recently, Li and Fang [22] have used a version of LB meth-

od (Shan–Chen multiphase model [20]) to simulate the electrow-
etting of an electrolyte droplet in two dimensions. They have
proposed a LB method to investigate the behavior of the contact
angle with applied potential.
This paper aims to extend the free energy-based LB approach
for modeling and simulation of the dynamics of a microsized fluid
droplet in the presence of an electric field (i.e., electrowetting phe-
nomenon) in three dimensions. Several LB algorithms for a liquid–
Fig. 2. (a) Droplet on an electrode. (b) Droplet on a charged substrate. (c) Droplet on gas system have been reported in the literature. Among there are
a dielectric substrate. the free energy-based LB approach, which first introduced by Swift
et al. [23,24], has advantages for wetting problems. Of course, the
original version of the method suffered from the lack of Galilean
motion of a droplet under the influence of an electric field and to invariance, a serious drawback when hydrodynamic transport be-
determine the optimized potential, thereby to increase the speed comes a relevant issue. This problem was solved by Holdych
of droplet motion. Of course, experiments can address the above et al. [25] who proposed a modified expression for the relation be-
questions, but experimental work on microsized and nanosized tween the pressure tensor and the second moments of the popula-
droplets can be expensive and difficult because of the length and tion densities in the lattice Boltzmann model.
time scales involved and also fabricating the electrowetting sys- Recently, this method has been successfully employed to model
tems usually takes more than eight months per device [14]. There- microscale droplets on the flat and superhydrophobic surfaces
fore the numerical modeling is preferable to study of behavior of [26–28] and to investigate the behavior of the droplets on chemi-
these systems. cally patterned substrates [29]. The significant advantage of this
It should be noted that, the numerical methods should be able model is the ease to controlling the wetting properties of substrate
to predict the droplet behavior as a function of liquid properties such as wetting potential [26]. Therefore, this version of LB method
and surface chemistry. This kind of prediction is quite useful for was selected to model and simulate the behavior of microscale
a parametric investigation of the physical behavior and also for de- droplet shape and motion due to the electrowetting.
sign and interpretation of experiments. This paper is organised as follows: in the next section the basic
In the literature, there are a few three-dimensional investiga- concept of the free energy LB model is briefly reviewed. Section 3
tions of droplet base electrowetting. Recently, Lienemann et al. gives the definition of the total free energy of the system in the
[15] have studied the droplet base electrowetting via quasi-static presence of electric field, to obtain new relations for the surface
method (i.e., no internal fluid effects). Bahadur and Garimella tensions using the variational calculus method. The PB equation
[16] have assumed the droplet as a rigid body to simplify the used and used numerical method to solve the discretized PB equation,
three-dimensional analyses. In addition, Decamps and Coninck to find the droplets in which potential distribution is one-dimen-
[17] have proposed a dynamic model of the contact angle variation sional, are presented in Section 4. The boundary conditions used
for an axisymmetricly spreading drop and Lu et al. [18] used a dif- in the simulations are briefly described in Section 5. Results are
fuse interface model to simulate the droplet motion. They have presented and discussed in Section 6. Finally, Section 7 contains
compared their results to experiments on a scaled-up version of some conclusions.
the electrowetting device. Also, Walker and Shapiro [19] have pre-
sented a computational model for electrowetting, but because of
using two different techniques, it is difficult to understand and 2. Lattice Boltzmann model
implement. For more details, the reader is referred to [14].
In the past 10–15 years, the LB method has been successfully In LB algorithm, the Navier–Stokes equations are solved via fol-
applied to simulate fluid flows and transport phenomena. Unlike lowing the evolution of a set of distribution functions fi(r, t). Here,
the conventional computational fluid dynamics (CFD) methods, fi(r, t) represents the mass density at time t and position r with
the LB method has an intrinsic kinetic nature [20] which leads to velocity ci. In the present study, a 3D lattice which has 15 velocity
the some advantages: reduction from second order to first order vectors, so-called D3Q15, is employed (see Fig. 3). The selected lat-
partial differential equations, simplification of nonlinear modeling, tice has the following velocity vectors ci: (0, 0, 0), (±1, ±1, ±1),
computational efficiency and accuracy, simple fluid interface (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) in lattice units. Physical quantities are
boundary conditions, and a mathematical framework allowing defined as moments of fi(r, t). Thus, the particle density and
molecular level modeling, therefore it is effective for problems in momentum are obtained by
which both mesoscopic and microscopic statistics become impor- X X
tant [21].
q¼ fi ; qua ¼ fi cia ; ð3Þ
3854 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

w1 = 1/3 and w2 = 1/24 are constants. A convenient choice for bulk

pressure is found in [31]

pb ¼ pc ðmp þ 1Þ2 ; ð9Þ

here, mp = (q–qc)/qc is the reduced density and pc = 1/8 and qc = 3.5

are the critical pressure and critical density, respectively.

3. Free energy definition

To study the electric double-layer from the free energy point of

view, three energies have to be considered [32,33]: (1) The electro-
static energy of the surface charges and the bulk charges, (2) The
chemical preference of the ions forming the surface charges for
Fig. 3. Schematic view of D3Q15 lattice. the surface over the bulk or the tendency to electrolytic dissocia-
tion of groups in the surface, and (3) Mechanical free energy of
where i indicates the velocity directions and a is used to denote the system.
Cartesian directions. The electric contribution on the free energy of system, for z:z
LB collisions and particle displacements are governed by the LB electrolytes, can be defined as [32]
equation, with BGK1 [30] collision operator for the time evaluation Z Z h i
e
of distribution functions Welec ¼ rudC  ðgraduÞ2 þ P dX; ð10Þ
C X 2
1
fi ðr þ ci Dt; t þ DtÞ  fi ðr; tÞ ¼  ½fi ðr; tÞ  fieq ðr; tÞ: ð4Þ
s where X is the domain volume bounded by the surface C, u is the
In the above equation, Dt is the time step of the simulation, is the fieq potential, r is the surface charge and e is the electric permittivity of
equilibrium distribution function and is depended on local fluid prop- the droplet, P = 2nbkBT[cosh(bu)  1] is the osmotic pressure
erties and s is the relaxation parameter which is related to the kine- (where, nb is the number density of ionic species in the bulk region,
matic viscosity in D3Q15 for hydrodynamic simulations as follows k is the Boltzmann constant, T is the absolute temperature) and
[20]: b ¼ kze
BT
(here, z is the valance of ionic species, and e is the electronic
  charge).
2s  1 Dr 2 The chemical free energy can be considered as follows [32]:
t¼ ; ð5Þ
6 Dt Z
where Dr is the lattice spacing. In the present work, the power ser- Wchem ¼  rudC: ð11Þ
C
ies proposed by Swift et al. [23] is used to calculate fieq
In the constant potential case, the chemical free energy loss of the
fieq ¼ A þ Bcia ua þ Cu2 þ Dcia cib ua ub þ Gab cia cib bulk electrolyte due to charge transfer to the surface is exactly bal-
ð6Þ
f0eq ¼ A0 þ C 0 u2 ; anced by the increase of the electrostatic free energy (i.e., the sur-
face integral in Eq. (10)).
here, b will be used to denote Cartesian directions and the usual Finally, the mechanical free energy in the absence of electric
summation convention over repeated indices is assumed. field can be represented by [33]
To determine the coefficients of Eq. (6), this equation should
Z " # Z
conserve mass and satisfy the three first moments of fieq con- j ðgradqÞ2
straints, i.e., the following equations: Wmech ¼ wb ðT; qÞ þ dX þ /ðqÞdC: ð12Þ
X 2 2 C
X
fieq ¼ q;
X In the above equation, the surface integral term represents the contri-
fieq ci ¼ qua ; bution of solid–fluid interactions and wb is the free energy in the bulk
X [31]
fieq cia cib ¼ Pab þ qua ub þ t½ua @ b ðqÞ þ ub @ a ðqÞ þ uc @ c ðqÞdab ;
X c2 q wb ðT; qÞ ¼ WðT; qÞ þ /b ðTÞq  pb ðTÞ; ð13Þ
fieq cia cib cic ¼ ðua dbc þ ub dac þ uc dab Þ;
3 where, /b and pb are the chemical potential and pressure in the bulk,
ð7Þ and W is a non-negative function of q that vanishes, along with @W ,
@q

where Pab is the pressure tensor. Therefore, the coefficients of Eq. when q is equal to the liquid bulk density or to the gas bulk density
(6) for D3Q15 lattice are obtained as follows: and it can be chosen the simplest excess free energy function [31]

w j W ¼ pc ðt2  ksp Þ2 ; ð14Þ

A¼ ðp  ð@ a ðqÞ2  jq@ aa q þ tua @ a qÞ;
c2 b 2
wq wq 3wq where sp ¼ T cTT
c
is the reduced temperature, Tc = 4/7 is critical tem-
B¼ 2 ; C ¼ 2; D¼ ; perature, and k is a constant and is related to the density difference
c 2c 2c4
1 between the two phases and will be set to 0.1 in this work.
G1cc ¼ 4 ðjð@ c qÞ2 þ 2tuc @ c qÞ; ð8Þ
In the constant potential case (i.e., when a constant potential is
2c
applied to the substrate), the total free energy of the system can be
G2cc ¼ 0;
represented as follows:
1 Z Z n
G2cd ¼ ðjð@ c qÞð@ d qÞ þ tðua @ d q þ ud @ c qÞÞ; j he io
16c4 Wtu ¼ /ðq; uÞdC þ wb þ ðgradqÞ2  ðgraduÞ2 þ P dX:
C X 2 2
where j is a positive constant which is related to the surface
tension and tunes interface width, pb is the bulk pressure and
ð15Þ
It should be noted that in presence of the electric field, although ini-
1
Bhatnagar, Gross, and Krook. tially / is assumed to be a function of both the density and potential
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3855

(i.e., in Eq. (15)), but it will be shown that in the constant potential Since it will be focused on the perpendicular direction to the flat
case, this term is only dependent on the density. substrates (i.e., z-direction), hereafter the derivatives in that direc-
It should be mentioned that, the potential distribution inside tion should be handled in such a way that the wetting properties of
droplet must satisfy the nonlinear PB equation [32] the substrate can be controlled. Now, by following Cahn [35], / (qs)
  can be taken as
j2D
r2 u ¼ sin hðbuÞ; ð16Þ
b 1
/ðqs Þ ¼ /0  /1 qs þ /2 q2s þ . . . ð26Þ
b 2 2 1=2 2
where j ¼ ð2n z e =ekB TÞ
1
D denotes the Debye length.
To minimization of the obtained total free energy relation (i.e., In the above equation, /1 is called the wetting potential and the sec-
Eq. (15)) using the calculus of variations method [34], it is assumed ond order terms of qs can be neglected [35].
that the variables are changed as follows: q ? q + g, Writing Eq. (22) in one-dimensional form and comparison it
rq ? rq + rg, u ? u + b, and ru ? ru + r b, so differential of with Eq. (13), one can obtain
Eq. (15) gives  2
j dq
Z  ! !
 ¼ wb  /b q þ pb : ð27Þ
@F ~ ~ F þ b @F þ rb r
~ ~ F dX 2 dz
dWtu ¼ g
þ rg r rq ru
@q @u
Z   Now using Eqs. (13) and (27), the following relations will be found:
@/ @/  2
þ g þb dT ¼ 0; ð17Þ
j dq
@q @u
h i W¼ ð28aÞ
2 dz
where  F ¼ wb þ j2 ðgradqÞ2  2e ðgraduÞ2 þ P and  1=2 ð28Þ
~ ~  @@n~
r i þ @@@n~j þ @@@n ~
k . Using divergence theorem gives j
rn @ dz ¼  dq ð28bÞ:
Z
@x
 
@y
Z 
@z
 2wb
t @F ~ ~ @F ~ ~
dWu ¼ g  r  rr ~ q F dX þ b  r  rr~ u F dX In Eq. (28b) the plus sign of dz is taken. Similarly, Eq. (24) can be
@q @u
Z Z   written in one-dimensional form as follows:
@/
~ ~ F þ ÞdC þ b ~ ~ ~ F þ @/ dC ¼ 0;
þ gð~
ns  r rq ns  r ru  2
@q @u e du
P¼ þ const: ð29Þ
ð18Þ 2 dz

for the above functional, in two-dimensional Euler’s equations are Using the boundary condition at infinity, Eq. (29) can be rewritten
represented as follows: as follows:
8      2
>
<  @@Fq þ @x
@ @F @
þ @z @F
¼ 0 ð19aÞ
e du
@ qx @ qz P¼ ð30aÞ
    ð19Þ 2 dz ð30Þ
>
:  @@Fu þ @x
@ @F @ @F
¼ 0 ð19bÞ:  1=2
@u
þ @z @u e
x z dz ¼  du ð30bÞ:
2P
Using Eqs. (19a), (15) and (13) leads to
In Eq. (30b), the minus sign of dz should be chosen. Its reason will be
@ discussed later (i.e., after obtaining new relations for the surface
ðw Þ ¼ 2/b þ jðqxx þ qzz Þ; ð20Þ
@q b tensions). Substituting Eqs. (28a) and (30a) into Eq. (15) leads to
Z
the integral of Eq. (20) is taken as follows:
Wtu ¼ /0  /1 qs þ fW þ W  ½PðuÞ þ PðuÞgdz
Z Z
Z Z
dðwb Þ ¼ 2/b dq þ j
¼ /0  /1 qs þ 2 Wdz  2 Pdz: ð31Þ
Z     Z    
@q @q @q @q
 d þ d : ð21Þ
@x @x @z @z Employing Eqs. (28b) and (30b) and changing the variable of inte-
grations from z to q and u and substituting into Eq. (31), one can
Finally, Eq. (21) gives obtain
j Z pffiffiffiffiffiffiffiffiffiffiffiffi Z pffiffiffiffiffiffiffiffiffiffi
wb ¼ 2/b  q þ ðqx þ qz Þ2 þ const: ð22Þ
2 Wtu ¼ /0  /1 qs þ 2jW dq þ 2ePdu: ð32Þ
Now, if one considers Eq. (19b), similarly
Since the minimum value of Wtu is the equilibrium surface free en-
@P
¼ er2 u; ð23Þ ergy [35], the surface tensions at the interfaces are given by
@u
Z ql pffiffiffiffiffiffiffiffiffiffiffiffi Z ul pffiffiffiffiffiffiffiffiffiffi
integration of the above equation gives rgl ¼ 2jW dq þ 2ePdu;
qg ug
e Z qg pffiffiffiffiffiffiffiffiffiffiffiffi Z ug pffiffiffiffiffiffiffiffiffiffi
P ¼ ðruÞ2 þ const: ð24Þ
2 rsg ¼ 2jW dq þ 2ePdu; ð33Þ
q us
Imposing the following conditions: Z qs l pffiffiffiffiffiffiffiffiffiffiffiffi Z ul pffiffiffiffiffiffiffiffiffiffi
rsl ¼ 2jW dq þ 2ePdu þ /0  /1 qs :
~ ~ F þ @/ ¼ 0 ! @/ ¼ jrq ð25aÞ;
ns  r
~ qs us
rq
@q @q
ð25Þ As mentioned in Section 1, since the applied voltage leads to de-
~ @/
ns  rr
~ ~u F þ ¼ 0 ð25bÞ; crease the value of rsl and us is greater than ul, in Eq. (30b) the
@u
minus sign of dz was chosen. It should be mentioned that, because
and considering ru = 0 on the substrate (because of using constant in the present study the applied potentials will be less than 1V, the
potential on the substrate), one can obtain @@/u ¼ 0. Therefore in this linearized osmotic pressure can be used (i.e., P  12 ej2D u2 Þ. There-
case (i.e., the constant potential) / will be independent of u. fore, the electric contribution in Eq. (33) can be obtained as follows:
3856 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

Z ul pffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
pffiffiffiffiffiffi Z ul 1 1
ul To solve Eq. (40), it should first be discretized. Discretization of Eq.
2ePdu ¼ 2e ej2D u2 du ¼ ejD u2 (40) gives
us us 2 2 us
  !
us
1 1
¼ ejD u2 : ð35Þ Ui;j ¼
2 ul
2
DR2
þ i2 DR22 Dh2 þ ðjD  r c Þ2
Now the famous Young–Laplace equation will be used to relate the Uiþ1;j þ Ui1;j Uiþ1;j  Ui1;j Ui;jþ1 þ Ui;j1
 þ þ 2 : ð41Þ
contact angle to the modified surface tensions associated with a so- DR 2 i  DR2 i  DR2  Dh2
lid, liquid, and surrounding fluid
It should be noted that the Debye length is an important length in this
rsg  rsl study. Since it is a very small quantity, it must be considered in the
cos h ¼ : ð36Þ
rgl grid generation (i.e., by grid size). Using suitable boundary conditions,
Eq. (41) can be solved via Gauss–Seidel iteration procedure [37].
Finally the system of equations to study the behavior of droplet un-
der influence of electric field is closed by the following relation [31]:
5. Boundary conditions
pffiffiffiffiffiffiffiffi p   
‘
 1=2
‘
/1 ¼ 2bsp jpc sign  h cos 1  cos : ð37Þ
2 3 3 In the present study, there are two sets of boundary conditions
(BCs): hydrodynamic BCs required for LB simulation and electrical
Here ‘ = cos1(sin2h) and sign(x) gives the sign of x.
BCs required to solve the PB equation. Hydrodynamic BCs in turn
include three important conditions: periodic, no-slip and wetting
4. Poisson–Boltzmann equation discretization
conditions. Wetting and no-slip boundary conditions are applied
on the substrate (i.e., bottom side of Fig. 4) using Eq. (25a) and
As it can be seen in Eq. (35), in order to compute the electric
bounce back [38] rule, respectively. Also, the no-slip BC (via the
contribution, ul must be already estimated. It can be assumed that,
bounce back rule) is used on the top side and periodic BC [38] on
the potential distribution inside a microsized droplet varies as one-
the lateral sides (see Fig. 4).
dimensional (i.e., just in perpendicular direction to the substrate).
To investigate the validity of this assumption, the PB equation
should be solved inside the droplets for different shapes (i.e., semi-
spherical and non-spherical in the spreading and motion states,
respectively). As mentioned before, the applied potential is less
than 1, therefore the linear PB equation can be used in the compu-
tational domain (i.e., inside the droplet). It should be denoted that
ul is the potential of a position very close to the solid surface
(10 Å [36]).
Since the droplet geometry is considered as axisymmetric, two-
dimensional linear PB equation is employed

@2u 1 @u 1 @2u
þ þ ¼ j2D  u; ð38Þ
@r 2 r @r r 2 @h2
where 0 < r < rc,0 < h < p and rc is the contact radius. The variables
can be put in dimensionless forms as follows:
R ¼ r=r c ; U ¼ u=us : ð39Þ
Here, us is the applied potential. Substituting the above variables
into Eq. (38) leads to

@2U 1 @U 1 @2U
2
þ  þ  ¼ ðjD  rc Þ2  U: ð40Þ
@R R @R R2 @h2 Fig. 5. Electric boundary conditions to simulate: (a) spreading and (b) motion.

Fig. 4. Computational domain and hydrodynamics boundary conditions.

 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3857

Electrical BCs used to solve the PB equation inside the droplet This iteration procedure is repeated until the absolute values of all
have been presented in Fig. 5. It can be seen that, a constant poten- errors fall below a pre-specified stopping criteria, i.e.,
tial is applied at the substrate (i.e.,u = us). Since the ambient fluid is unew old
i;j  ui;j < 107 . The essential coefficients in the PB equation
air, the no-charge BC at air–solid and liquid–air interfaces have were considered as follows: the Debye length j1 D ¼ 9:6 nm, the
2
been considered. electric permittivity e ¼ 7:26 1010 CJm for aqueous electrolytes,
and b ¼ kze
BT
¼ 38:9 V1 for monovalent ionic species [39].
6. Results and discussion As mentioned in Section 1, the main purpose of this section
is finding droplets in which potential distribution are one-
6.1. Potential distribution inside the droplet dimensional, by solving the PB equation numerically inside
them. As in all numerical solution checking the grid indepen-
As mentioned in Section 4, the Gauss–Seidel iteration procedure dence of the results after converging of the solution (e.g., in
is employed to obtain the potential distribution inside the droplets. the present study based on the above-mentioned stopping crite-

Fig. 6. Potential variation diagram for R = 10 nm at: (a) x = 0.5R and (b) x = 0.9R.
3858 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

ria). is the final step of the numerical solution. Along with exponential and nearly one-dimensional behavior. This is the
investigating about the distribution of electric potential inside major difference between the results of case 10 nm and
different droplets, also the grid independence of the results 100 nm ones.
has been presented. For three size of droplets, potential varia- Finally, for a droplet with radius of 1000 nm, Fig. 8 indicates
tion has been studied. that if grid size approaches to 2 nm the results will be grid
For a droplet with the radius of 10 nm, Fig. 6a and b show that independent. Also it is evident from this Figure that, the poten-
when the grid size approaches to 0.2 nm in, the results become grid tial decay length compared to droplet radius, is very small.
independent. Furthermore, as it can be seen in Fig. 6, in this case Therefore, the assumption (i.e., the potential distribution is
potential variation is sharp, therefore it cannot be considered one-dimensional) is exactly valid for the semi-spherical shape
one-dimensional. micro sized droplets. To give a better description Fig. 9 shows
In a droplet with radius of 100 nm, based on Fig. 7a and b, to two-dimensional potential contours inside the three above-men-
check the grid independence of the results the grid size should tioned droplets. As it can be seen, by increasing the droplet
be decreased to 0.5 nm. In this case, potential variation has size, potential distribution inside the droplets become exactly

Fig. 7. Potential variation diagram for R = 100 nm at: (a) x = 0.5R and (b) x = 0.9R.
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3859

Fig. 8. Potential variation diagram at x = 0 and R = 1000 nm.

Fig. 9. Potential contours inside droplets with radius of 10, 100 and 1000 nm.

one-dimensional. Also, this Figure represents more clearly the The validity of one dimensionality of potential distribution for
process in which the results become grid independent in each droplets with non-spherical shape (i.e., droplets during motion)
case. was also examined. As in the case of semi-spherical case, the
3860 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

2
assumption is not valid for 10 nm droplet. Fig. 10a and b indicate d uðzÞ j2D
that for droplet with radius of 100 nm potential distribution is 2
¼ sin h½buðzÞ; ð42Þ
dz b
two-dimensional, but for droplet with radius of 1000 (i.e., micro-
size droplet), Fig. 11a and b show that the potential distribution
will be one-dimensional. Furthermore, in the different shape of the solution of the above equation gives the variation of the electric
the microsize droplets, potential variation has the same curve, potential near the substrate [40]
see Fig. 12.
" #
eu0 =2 þ 1 þ eu0 =2  1 ef
6.2. Lattice Boltzmann simulations u ¼ 2Ln ; ð43Þ
eu0 =2 þ 1  eu0 =2  1 ef
Since it was shown that the potential distribution for the micro-
sized droplets is one-dimensional (i.e., in Section 6.1), one can esti- where u ¼ bu; u0 ¼ bus , and f = jDz. Therefore, using Eqs. (33),
mate the magnitude of ul by the one-dimensional nonlinear PB (35) and (36) the contact angle in presence of an electric filed is cal-
equation as follows: culated by

Fig. 10. Potential variation diagram at different x: (a) R = 75 nm and (b) R = 50 nm.
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3861

Fig. 11. Potential variation diagram at different x: (a) R = 500 nm and (b) R = 250 nm.

cos h ¼ cos hY Here, Upzc is the potential of zero charge and dH is the distance where
  bu =2 2 ! the counter-ions are all located at this fixed distance (i.e., the Debye
ejD 2 e s þ 1 þ ðebus =2  1ÞejD z
þ u2s  Ln bus =2 : length). The obtained results were compared with this correlation
2rgl b e þ 1  ðebus =2  1ÞejD z ones for typical values of dH (i.e., j1 D Þ ¼ 10 nm; e ¼ 81 8:85
ð44Þ 1012 C 2 =Jm; rlg ¼ 0:072 N=m and also it was assumed that in the ab-
sence of an applied voltage Upzc is zero.
As mentioned in Section 1, the results have been compared with the It should be noted that, in the simulations, two-phases medium
analytical correlation presented by Mugele and Baret [1]. For an has been considered, so liquid vapor is taken as the ambient fluid
electrolyte droplet placed directly on an electrode surface, they (instead of air) to avoid complexity. In addition, a diffuse interface
found model has been used for the interface and the width of which is
calculated as [31]
e qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
cos h ¼ cos hY þ ðU  U pzc Þ2 : ð45Þ v¼ jq2c =4ksT pc : ð46Þ
2dH rlg
3862 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

Fig. 12. Potential variation diagram for three different shape of microsize droplets at x = 0.

Fig. 13. Spreading of a droplet on 60 60 40 lattice, hY = 80°, and us = 0, j = 0.001.

 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3863

Normally, in the LB simulations the dimensionless lattice units are The comparison of the results with correlation ones in this case
used. To convert these units into physical ones, one needs to define also shows good agreement, see Fig. 17.
some conversion factors. The time conversion factor can be defined Fig. 18a indicates sequential frames of the droplet profiles
as follows: during spreading due to the voltage applying. Fig. 18b and c
show the velocity vector filed inside the droplet before and after
t 2
t0 ¼ L ; ð47Þ applying voltage. An important feature of computer simulation
t 0 of fluid dynamical problems is that not only one has access to
where t* is the dimensionless kinematic viscosity, t is the kine- the system density at each point in space but also the whole
matic viscosity and L0 ¼ rr is the length factor (here, r is experi- velocity field is known. This is an interesting property since it al-
mental droplet radius and r* is the dimensionless droplet radius lows a survey of the liquid motion inside the droplet, whereby
used in the simulations). The required data to compute a suitable providing a means to estimate the viscous dissipation. Since
time conversion factor are given in the table in the corresponding the interface is considered as diffusive, in the zero potential state
section. (i.e., Fig. 18b) at some points in the interfacial region velocity
vectors are still seen.
6.2.1. Droplet spreading
Fig. 13 presents LB simulation of spreading of a droplet on a the 6.2.2. Droplet motion
substrate in the absence of electric field. To study the spreading of Using two electrodes as indicated in Fig. 19 and applying
droplet due to the electrowetting phenomenon, two cases have a voltage at one side leads to droplet motion. When the
been examined. First, a droplet with the initially prescribed contact right side electrode is made active, contact angle on that
angle hY = 80° (i.e., the hydrophilic surface) and then with hY = 120° side is reduced. With this reduction, because of mass bal-
(i.e., the hydrophobic surface), which are directly placed on the ance, droplet starts to move to find a new stable situation.
electrode surface, were considered. The predicted behavior of the This is the basic idea of pumps which work by electrowett-
droplet shape with hY = 80°, Fig. 14a and b under the variation of ing phenomenon.
the applied voltage has been presented in Fig. 14c–f. As it can be Fig. 20 shows the experimental results of Lu et al. [18]. They
seen, by increasing the applied voltage, the contact angle is re- have used 60% glycerin–water mixture as the droplet (the den-
duced. The obtained results have been compared with the analyt- sity and viscosity of the mixture are 1.1538 103 kg m3 and
ical correlation ones presented by Mugele et al [1], in Fig. 15. It can 0.02030 Pa s at 20°C, respectively). In their experiment, droplet
be seen that, there are very good agreement between the results. has been placed on the dielectric substrate, therefore the high
Fig. 16 indicate the simulation results of behavior of the droplet voltage should be applied for droplet transporting (i.e., 50 V)
shape under the variation of the applied voltage with hY = 120°. respect to the droplet directly placed on the electrode case (see

Fig. 14. Variations of the droplet shape under the influence of the applied voltage, j = 0.001: (a) and (b) hY = 80° and us = 0, (c) us = 150, (d) us = 250, (e) us = 350, and (f)
us = 450 mV.
3864 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

Fig. 15. Comparison of the variation of contact angle with analytical ones.

Fig. 16. Variations of the droplet shape under the influence of the applied voltage, j = 0.001: (a) and (b) hY = 120° and us = 0, (c) us = 150, (d) us = 250, (e) us = 350, and (f)
us = 450 mV.

Section 1). It should be noted that in the experiment the contact sults with experimental ones, presented in Fig. 20, the strength of
angle has been changed nearly from 100° to 70°due the applied the simulations will be clearer. As it can be seen, the droplet
potential. behavior (from time and shape point of view) is qualitatively sim-
Fig. 21 give the top views of the simulation results for the drop- ilar to that observed experimentally. Also, the side views of the
let motion. To calculate time in per indicated stage, the data in Ta- motion have been presented in Fig. 22. It should be mentioned that
ble 1 were used. Using the mentioned data, the time conversion since Lu et al. have sandwiched the droplet with a grounding plate
factor is equal to t0 = 3.35 106. By comparing the simulation re- on top (i.e., two plates state), their reported times are a bit more
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3865

Fig. 17. Comparison of the variation of contact angle with analytical ones.

Fig. 18. Droplet profiles and velocity field during spreading due to the applied voltage: (a) side views of sequential frames of the droplet profiles during spreading;
h:120° ? 90°. (b) Rest state and (c) velocity vector field.

than the present ones (i.e., one plate state). It is clear that in the
sandwiched case, the friction force will be greater than the single
planer surface.
Finally, sequential frames of the droplet profiles during motion
process from the start to the end of motion have been demon-
strated in Fig. 23a. Also, Fig. 23b–f show some views of the velocity
vector fields during this translation.

7. Conclusion

In this paper, first the validity range of one dimensionality

Fig. 19. Schematic diagram of a droplet transport by electrowetting. of the potential distribution inside the droplets by solving
3866 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

Fig. 20. Top views of experimental results of droplet motion, hY = 100° and h = 70°.

Fig. 21. Top views of the droplet movement due to the electrowetting, hY = 100° and h = 70°, j = 0.001.

Table 1
The required data to determine the time conversion factor for the droplet motion. modeling and simulation of the electrowetting phenomenon by
* a *
the LB method were presented. New correlations for the sur-
Lattice t = 0.133 r = 16 LB unit
Metric t = 1.76 105 m2/s r = 0.55 mm
face tensions are presented to apply the electric contribution
on the wetting potential via the contact angle. In continuation,
a 2s1 Dr 2
s ¼ 0:9 ! t ¼ 6 Dt ¼ 0:133. the droplet spreading and motion due to electrowetting were
simulated in three dimensions. Comparison of the contact
angle variation and time scale of motion with the analytical
and experimental ones shows good agreement. It is concluded
the PB equation was investigated. It was found that for the that the proposed LB method can be used as a powerful
microsized droplets, one can consider the potential distribution method to study the electrowetting phenomenon in three
exactly one-dimensional in all cases. Then a free energy based dimensions.
 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868 3867

Fig. 22. Side views of the droplet movement due to electrowetting, hY = 100° and h = 70°, j = 0.001.

Fig. 23. Droplet profiles and velocity field during motion due to the applied voltage: (a) Top views of sequential frames of the droplet motion due to the voltage applying, (b)
velocity vector field at the rest, and (c–f) velocity vector fields during motion.
3868 H. Aminfar, M. Mohammadpourfard / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3852–3868

Acknowledgement [19] S. Walker, B. Shapiro, Modeling the fluid dynamics of electrowetting on

dielectric (EWOD), J. Microelectromech. Syst. 15 (2006) 986–1000.
[20] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev.
We thank F. Varnik for useful discussions and A.Dupuis for pro- Fluid Mech. 30 (1998) 329–364.
viding us a version of his LB code. [21] H. Aminfar, M. Mohammadpourfard, Lattice Boltzmann BGK model for gas flow
in a microchannel, Proc. IMechE C: J. Mech. Engrg. Sci. 222 (2008) 1855–1860.
[22] H. Li, H. Fang, Lattice Boltzmann simulation of electrowetting, Eur. Phys. J.
References Special Top. 171 (2009) 129–133.
[23] M.R. Swift, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann simulation of
[1] F. Mugele, J.C. Baret, Electrowetting: from basics to applications, J. Phys. nonideal fluids, Phys. Rev. Lett. 75 (1995) 830–833.
Condens. Matter 17 (2005) R705–R774. [24] M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann
[2] M. Preuss, H.J. Butt, Direct measurement of particle–bubble interactions in simulations of liquid–gas and binary fluid systems, Phys. Rev. E 54 (1996)
aqueous electrolyte: dependence on surfactant, Langmuir 14 (2004) 3164–3174. 5041–5052.
[3] K.H. Kang, I.S. Kang, C.M. Lee, Wetting tension due to Coulombic interaction in [25] D.J. Holdych, D. Rovas, J.-G. Georgiadis, R.O. Buckius, An improved
charge-related wetting phenomena, Langmuir 19 (2003) 5407–5412. hydrodynamics formulation for multiphase flow lattice-Boltzmann models,
[4] C. Quilliet, B. Berge, Electrowetting: a recent outbreak, Curr. Opin. Colloid Int. J. Mod. Phys. C 9 (1998) 1393–1404.
Interface Sci. 6 (2001) 34–39. [26] A. Dupuis, J.M. Yeomans, Lattice Boltzmann modelling of droplets on chemically
[5] J. Zeng, T. Korsmeyer, Principles of droplet electrohydrodynamics for lab-on-a- heterogeneous surfaces, Fut. Gen. Comput. Sys. 20 (2004) 993–1001.
chip, Lab. Chip. 4 (2004) 265–277. [27] Q. Chang, J.I.D. Alexander, Analysis of single droplet dynamics on striped
[6] V. Srinivasan, V.K. Pamula, R.B. Fair, An integrated digital microfluidic lab-on- surface domains using a lattice Boltzmann method, Microfluid. Nanofluid. 2
a-chip for clinical diagnostics on human physiological fluids, Lab. Chip. 4 (2006) 309–326.
(2004) 310–315. [28] A. Dupuis, J.M. Yeomans, Modelling droplets on superhydrophobic surfaces:
[7] V. Cristini, Y.C. Tan, Theory and numerical simulation of droplet dynamics in equilibrium states and transitions, Langmuir 21 (2005) 2624–2629.
complex flows, Lab. Chip. 4 (2004) 257–264. [29] Huang C. Shu, Y.T. Chew, Numerical investigation of transporting droplets by
[8] Y.C. Tan, J.C. Fisher, A.I. Lee, V. Cristini, A.P. Lee, Design of microfluidic channel spatiotemporally controlling substrate wettability, J. Colloid Interface Sci. 328
geometries for the control of droplet volume, chemical, concentration and (2008) 124–133.
sorting, Lab. Chip. 4 (2004) 292–298. [30] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases I
[9] A. Wixforth, C. Strobl, C. Gauer, A. Toegl, J. Scriba, Z. von Guttenberg, Acoustic small amplitude processes in charged and neutral one-component system,
manipulation of small droplets, Annal. Bioanal. Chem. 379 (2004) 982–991. Phys. Rev. 94 (1954) 511–525.
[10] O.D. Velev, B.G. Prevo, K.H. Bhatt, On-chip manipulation of free droplets, [31] A.J. Briant, A.J. Wagner, J.M. Yeomans, Lattice Boltzmann simulations of contact
Nature 426 (2003) 515–516. line motion. I. Liquid–gas systems, Phys. Rev. E 69 (2004) 031602–031615.
[11] P.R.C. Gascoyne, J.V. Vykoukal, J.A. Schwartz, T.J. Anderson, D.M. Vykoukal, [32] J.Th.G. Overbeek, The role of energy and entropy in the electrical double layer,
K.W. Current, C. Mc Conaghy, F.F. Becker, C. Andrews, Dielectrophoresis-based Colloids Surf. 51 (1990) 61–75.
programmable fluidic processors, Lab. Chip. 4 (2004) 299–309. [33] L.D. Landau, E.M. Lifshitz, Statistical Physics, Pergamon Press, London, 1958.
[12] D.W. Langbein, Capillary Surfaces: Shape-Stability-Dynamics, in: Particular [34] A. Cherkaev, E. Cherkaev, Calculus of Variations and Applications, Lecture
under Weightlessness, Springer, Berlin, 2002. Notes, University of Utah Salt Lake City, UT, 2003.
[13] H. Gau, S. Herminghaus, P. Lenz, R. Lipowsky, Liquid morphologies on structured [35] J.W. Cahn, Critical point wetting, J. Chem. Phys. 66 (1977) 3667–3672.
surfaces: from microchannels to microchips, Science 283 (1999) 46–49. [36] J. Zhang, D.Y. Kwok, On the validity of the Cassie equation via a mean-field
[14] S.W. Walker, Modeling, Simulating, and Controlling the Fluid Dynamics of free-energy lattice Boltzmann approach, J. Colloid Interface Sci. 282 (2005)
Electro-wetting on Dielectric, PhD Thesis, University of Maryland, College 434–438.
Park, USA, 2007. [37] J.D. Hoffman, Numerical Method for Engineers and Scientists, Marcel Dekker
[15] J. Lienemann, A. Greiner, J.G. Korvink, Modeling, simulation, and optimization, Inc., USA, 2001.
IEEE Trans. Comput. Aid. Des. Integr. Circ. Syst. 25 (2006) 234–247. [38] D.A. Wolf-Gladrow, Lattice Gas Cellular Automata and Lattice Boltzmann
[16] V. Bahadur, S.V. Garimella, An energy-based model for electrowetting-induced Models: An Introduction, Springer-Verlag, Heidelberg, 2005.
droplet actuation, J. Micromech. Microengrg. 16 (2006) 1494–1503. [39] K.H. Kang, I.S. Kang, C.M. Lee, Geometry dependence of wetting tension on,
[17] C. Decamps, J.D. Coninck, Dynamics of spontaneous spreading under charge-modified surfaces, Langmuir 19 (2003) 6881–6887.
electrowetting conditions, Langmuir 16 (2000) 10150–10153. [40] E.J.W. Verwey, J.Th.G. Overbeek, Theory of the Stability of Lyophobic Colloids,
[18] H.W. Lu, K. Glasnner, A.L. Bertozzi, C.J. Kim, A diffuse interface model for Elsevier Publishing Company Inc., 1948.
electrowetting droplets in a Hele–Shaw cell, J. Fluid Mech. 590 (2007) 411–435.