# Solution of Time Periodic Electroosmosis Flow with Slip Boundary

Authors: Qian Sun; Yonghong Wu; Lishan Liu; B. Wiwatanapataphee
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Recent research confirms that slip of a fluid on the solid surface occurs at micrometer scale. Slip on solid surface may cause the change of interior material deformation which consequently leads to the change of velocity profile and stress field. This paper concerns the time
periodic electroosmotic flow in a channel with slip boundary driven by an alternating electric field, which arises from the study of particle manipulation and separation such as flow pumping and mixing enhancement. Although exact solutions to various flow problems of electroosmotic flows under the no-slip condition have been obtained, exact solutions for problems under slip boundary conditions have seldom been addressed. In this paper, an exact solution is derived for the time periodic electroosmotic flow in two-dimensional straight channels under slip boundary conditions.


---

## Body

## 1.  Introduction

One of the major focuses in scientific research during the past decades has been the study of material behaviour at micro- and nanoscale. The subject of micro- and nanofluidics concerns mass and momentum transfer in micro- and nanoscale systems or around micro- or nanosized objects, as well as the design and application of biological and engineering devices and systems involving fluid flow through microchannels [1–3]. The study of microflows has attracted great attention from the science and engineering communities in the effort of developing a better understanding of the mechanism of microflows and consequently developing better models and control. Traditionally, in the model for fluid flows, the no-slip boundary condition is used; that is, the fluid velocity relative to the solid is assumed to be zero on the fluid-solid interface. However, it has been established that the no-slip condition is a hypothesis instead of a condition deduced from any principle. Evidence of slip of a fluid on the solid surface was reported by many researchers [4, 5].In recent years, time periodic electroosmosis has been studied intensively due to the recent development of microfluidic devices. Time periodic electroosmosis is known as AC electroosmosis, and the driving force is generated by an alternative electric field. Dose and Guiochon published numerical results for impulsive started electroosmotic flow [6]. Green et al. studied the electroosmosis on a planner microelectrode under both steady and unsteady fields [7]. Soderman and Jonsson investigated the electroosmotic flows for various geometries and derived various analytical results for the flow in two-dimensional microchannels and for the flow over a flat plate [8]. By using the green function method, Kang and coworkers [9] analysed the instantaneous velocity of the flow through a microchannel with time-dependent electric field. Jian et al. [10] studied the time periodic electroosmotic flow through circular microchannels. They discussed two different kinds of situations, which are the electroosmotic flow through the parallel plates and circular tubes.Although analytical and numerical solutions to many flow problems of Newtonian fluids have been established under the no-slip assumption [11–13], very few exact solutions for time-dependent problems have been obtained for the slip case [14, 15]. Hence, we investigate, in this paper, the time periodic electroosmotic flows in a microchannel under the Navier slip assumption. The rest of the paper is organized as follows. In Section 2, the mathematical equations for the channel flow are presented, which include the partial differential equations and the boundary conditions based on the slip assumption. In Section 3, the general solution for the velocity of the fluid in the channel is derived with some arbitrary constants. Then the slip boundary condition and symmetric condition are used to determine the arbitrary constants. Numerical investigations are then given in Section 4, followed by a concluding remark in Section 5.


## 2.  Problem Description and Mathematical Formation

The problem considered is the electrically driven flow of fluid through channels with height2
h. The geometry of the channels and the coordinates (
x
,
y
,
z
) used are shown in Figure 1.Figure 1
The flow channel and coordinate system used.Based on the principle of continuous mechanics, the governing equations for the flow of an incompressible ionized fluid are the standard Navier-Stokes equations:(1)
ρ
f
(
∂
V
→
∂
t
+
(
V
→
·
∇
)
V
→
)
=
-
∇
P
+
μ
∇
2
V
→
+
ρ
e
E
→
,
where P is the pressure, ρ
f and ρ
e are, respectively, the fluid density and the electric charge density, and V
→
=
(
u
,
v
,
w
) is the velocity vector with μ being the fluid viscosity. The last term on the right side denotes the electroosmotic body forces and E
→ denotes the externally applied electric field.We assume fully developed symmetric channel flow, and thus the swirling flow and consequently the velocity components iny and z directions can be neglected. This allows us to use the superposition principle because the resulted unsteady Navier-Stokes equations for this case are linear in velocity. This means that we can solve exactly the pressure and electroosmotically driven flows separately and then combine the two solutions to get the solution of the problem with the combined effects. The form of the electroosmotic body force, in general, depends on the form of the electric field applied on the system. In this paper, we will focus on sinusoidally driven flows.Taking into account the above assumptions, (1) reduces to
(2)
ρ
f
∂
u
∂
t
=
μ
∂
2
u
∂
y
2
+
ρ
e
E
x
sin
(
Ω
t
)
,
where E
x and Ω are, respectively, the magnitude and the frequency of the external electric field E
→. With (2), we can then determine the velocity of the flow driven by the time periodic electric field.Based on the theory of electric field [16], the ρ
e in (2) can be determined by
(3)
ρ
e
=
-
2
n
0
e
z
sinh
⁡
(
e
z
ψ
k
B
T
)
,
where n
0 denotes the average number of positive or negative ions in the buffer, z is the valance, e is the electron charge, T is the absolute temperature, ψ is the electric potential, and k
B is the Boltzmann constant. Let ψ
*
=
ψ
/
ξ be the normalized electroosmotic potential in which ξ is the zeta potential and let α
=
(
e
z
ξ
/
k
B
)
T be a parameter; then we have
(4)
ρ
e
=
-
2
n
0
e
z
sinh
⁡
(
α
ψ
*
)
.Based on the work in [16], ψ is governed by the Poisson-Boltzmann equation:
(5)
∇
2
ψ
=
-
ρ
e
ϵ
.
in which ϵ is the permittivity, and we have assumed that the effect of time fluctuation of the external electric field on the electric potential is negligible.It has been established in the literature that the effective thickness of the electric double layer (EDL) depends on the ionic energy parameterα. For example, for α
=
5, the effective electric double layer thickness is about 4.217 λ, where λ is related to other system parameters by
(6)
ω
=
1
λ
=
8
π
n
0
e
2
z
2
D
k
B
T
.In most microfluidic applications, the value ofλ is about 10–100 nm, which is generally much smaller than the height of a typical microchannel. Hence, the electroosmotic potential generally reaches zero in the channel center, and its distribution along the channel height can be expressed by [16, 17]
(7)
ψ
*
=
4
α
tanh
-
1
[
tanh
(
α
4
)
exp
⁡

(
-
χ
)
]
,
where χ
=
y
/
λ denotes the distance from the channel wall. From the definition of the Debye-Huckel parameter ω, we can write the electroosmotic body force as
(8)
-
2
n
0
e
z
E
x
sin
⁡
Ω
t
=
ω
2
μ
u
HS
α
sin
⁡
Ω
t
,
where u
HS
=
-
ϵ
ζ
E
x
/
μ is the Helmholtz-Smoluchowski velocity. Then, (2) can be written as
(9)
ρ
f
∂
u
∂
t
=
μ
∂
2
u
∂
y
2
+
ω
2
μ
u
HS
α
sinh
⁡
(
α
ψ
*
)
sin
⁡
Ω
t
.To completely determine the fluid flow in the channel, the differential equation (2) must be supplemented by a set of boundary and initial conditions. To describe the slip characteristics of a fluid on the solid surface, the Navier slip condition is used; namely, the velocity of fluid, relative to the solid surface, is assumed to be proportional to the sheer stress on the fluid-solid surface; that is,
(10)
u
(
0
,
t
)
+
l
∂
u
∂
y
(
0
,
t
)
=
0
,
where l is the so-called slip length. We should address here that, for l
=
0, condition (10) reduces to the no-slip boundary condition.In addition, due to the symmetry of velocity about the symmetric planey
=
h, we have
(11)
∂
u
∂
y
(
h
,
t
)
=
0
.Now the partial differential equation (9) together with conditions (10) and (11) constitutes a complete boundary value problem for the fluid flow in the channel.


## 3. Solution of the Boundary Value Problems

By introducing the characteristic time1
/
Ω and the characteristic length λ, (9) can be expressed in normalized form
(12)
ρ
f
u
HS
Ω
∂
U
∂
θ
=
μ
u
HS
ω
2
[
∂
2
U
∂
2
χ
+
sin
⁡
(
θ
)
α
sinh
⁡
(
α
ψ
*
)
]
,
in which θ
=
Ω
t denotes the nondimensional time, χ
=
y
/
λ denotes the nondimensional distance, and U
=
u
/
u
HS denotes the nondimensional velocity. From (12), we have
(13)
∂
U
∂
θ
=
1
κ
2
[
∂
2
U
∂
2
χ
+
sin
⁡
(
θ
)
α
sinh
⁡
(
α
ψ
*
)
]
,
where κ
=
(
Ω
λ
2
/
μ
)
1
/
2 is a normalized parameter expressed as a function of the Debye length (
λ
), the electric field excitation frequency (Ω), and the kinematic viscosity (
μ
). The parameter κ represents the ratio of the Debye length to a diffusion length scale (I
D) given by
(14)
I
D
=
λ
κ
.LetU
~ be the solution of the following equation:
(15)
∂
U
~
∂
θ
=
1
κ
2
[
∂
2
U
~
∂
2
χ
+
e
i
θ
α
sinh
⁡
(
α
ψ
*
)
]
.
Then, from the linearity and the superposition principle, the solution of (13) will be U
=
Im
⁡
(
U
~
). It is noted that (15) admits solution of the form U
~
=
e
i
θ
F
(
χ
) where we have written sin
(
θ
) as the imaginary part of exp
⁡
(
i
θ
), and thus we proceed to solve (15) first.The derivative ofU
(
χ
,
θ
) with respect to time can then be written as
(16)
∂
U
~
∂
θ
=
i
exp
⁡

(
i
θ
)
F

(
χ
)
and the second space derivatives can be written as
(17)
∂
2
U
~
∂
2
χ
=
exp
⁡

(
i
θ
)
d
2
F
d
χ
2
.
Therefore (15) becomes
(18)
i
exp
⁡

(
i
θ
)
F

(
χ
)
=
exp
⁡

⁡
(
i
θ
)
κ
2
[
d
2
F
d
χ
2
+
sinh
⁡
⁡
(
α
ψ
*
)
α
]
.
For the above to be true for all θ values, we require
(19)
d
2
F

(
χ
)
d
χ
2
-
i
κ
2
F

(
χ
)
=
-
sinh
⁡
⁡
(
α
ψ
*
)
α
.The general solution of (19) includes a complimentary solution and a particular solution. The complimentary solution can be derived by solving the associated homogenous equation, while the particular solution can be obtained by using the variation of parameter methods. Through some calculation, we have(20)
F

(
χ
)
=
A
exp
⁡

(
κ
χ
i
)
+
B
exp
⁡

(
-
κ
χ
i
)
+
1
2
κ
α
i
[
∫
0
χ

exp
⁡

(
κ
χ
i
)
sin
h
(
α
ψ
*
)
d
χ
exp
⁡

(
-
κ
χ
i
)
1111111111
×
∫
0
χ

exp
⁡

(
κ
χ
i
)
sinh
⁡
(
α
ψ
*
)
d
χ
]
-
1
2
κ
α
i

exp
⁡

(
κ
χ
i
)
×
∫
0
χ

exp
⁡

(
-
κ
χ
i
)
sinh
⁡
(
α
ψ
*
)
d
χ
,
and consequently we have
(21)
U
~
(
χ
,
θ
)
=
A
exp
⁡

(
i
θ
)
exp
⁡

(
κ
χ
i
)
+
B
exp
⁡

(
i
θ
)
exp
⁡

(
-
κ
χ
i
)
+
exp
⁡
(
i
θ
)
2
κ
α
i
[
∫
0
χ
exp
⁡

(
κ
χ
i
)
sin
h
(
α
ψ
*
)
d
χ
exp
⁡

(
-
κ
χ
i
)
1111111111
×
∫
0
χ
exp
⁡

(
κ
χ
i
)
sinh
⁡
(
α
ψ
*
)
d
χ
]
-
exp
⁡
(
i
θ
)
2
κ
α
i
[
exp
⁡

(
κ
χ
i
)
∫
0
χ
exp
⁡

(
-
κ
χ
i
)
111111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
∫
0
χ
‍
]
.Now we proceed to determine the arbitrary constants in (21). In terms of the nondimensional time θ, the nondimensional distance χ, and the nondimensional velocity U, the boundary conditions (10) and (11) become
(22)
U
(
0
,
θ
)
+
l
λ
∂
U
∂
χ
(
0
,
θ
)
=
0
,
∂
U
∂
χ
(
h
λ
,
θ
)
=
0
.AsU
(
χ
,
θ
)
=
Im
⁡
(
U
~
(
χ
,
θ
)
), (22) will be satisfied if U
~ satisfies the following boundary conditions:
(23)
U
~
(
0
,
θ
)
+
l
λ
∂
U
~
∂
χ
(
0
,
θ
)
=
0
,
∂
U
~
∂
χ
(
h
λ
,
θ
)
=
0
.
Noticing that U
~
=
e
i
θ
F
(
χ
), the above boundary conditions become
(24)
F

(
0
)
+
l
λ
F
′

(
0
)
=
0
,
F
′

(
h
λ
)
=
0
.To apply the boundary condition (24), we first need to calculate the spatial derivative of F
(
χ
). From (7) and (20), we have
(25)
F
′

(
χ
)
=
A
(
κ
i
)
exp
⁡

(
κ
χ
i
)
+
B
(
-
κ
i
)
exp
⁡

(
-
κ
χ
i
)
-
1
2
α
[
exp
⁡

(
-
κ
χ
i
)
∫
0
χ

exp
⁡

(
κ
χ
i
)
1111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
∫
0
χ
(
exp
(
κ
χ
i
)
]
-
1
2
α
[
exp
⁡

(
κ
χ
i
)
∫
0
χ

exp
⁡

(
-
κ
χ
i
)
111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
exp
⁡
(
κ
χ
i
)
∫
0
χ
‍
(
exp
(
-
κ
χ
i
)
]
.At the centre of the channel,χ
=
h
/
λ and d
F
/
d
χ
=
0 and so from (24) and (25) we have
(26)
F
′

(
h
λ
)
=
A
(
κ
i
)
exp
⁡

(
κ
h
λ
i
)
+
B
(
-
κ
i
)
exp
⁡

(
-
κ
h
λ
i
)
-
1
2
α
exp
⁡

(
-
κ
h
λ
i
)
∫
0
h
/
λ

exp
⁡

(
κ
χ
i
)
111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
-
1
2
α
exp
⁡

(
κ
h
λ
i
)
∫
0
h
/
λ

exp
⁡

(
-
κ
χ
i
)
1
1111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
=
0
.From (20), (25), and the first equation of (24), we get
(27)
(
1
+
l
λ
κ
i
)
A
+
(
1
-
l
λ
κ
i
)
B
=
0
and hence
(28)
B
=
-
(
1
-
l
λ
κ
i
)
-
1
(
1
+
l
λ
κ
i
)
A
.Substituting (28) into (26) yields
(29)
A
=
{
exp
⁡

(
κ
h
λ
i
)
+
(
1
-
l
λ
κ
i
)
-
1
11
×
(
1
+
l
λ
κ
i
)
exp
⁡
(
-
κ
h
λ
i
)
exp
⁡

(
κ
h
λ
i
)
+
(
1
-
l
λ
κ
i
)
-
1
}
-
1
×
{
1
2
α
κ
i
exp
⁡

(
-
κ
h
λ
i
)
∫
0
h
/
λ

exp
⁡

(
κ
χ
i
)
111111111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
111
+
1
2
α
κ
i
exp
⁡

(
κ
h
λ
i
)
∫
0
h
/
λ

exp
⁡

(
-
κ
χ
i
)
11111111111111111111111
×
sinh
⁡
(
α
ψ
*
)
d
χ
1
2
α
κ
i
}
.


## 4. Numerical Investigation

In this section, we investigate the characteristic of the electric driven flow. In particular, we will investigate the influence of various model parameters on the flow behaviour, including theκ value, the energy parameter α, the channel height h, and the slip parameter l. Table 1 shows the solutions of A and B for various sets of the model parameters.Table 1
A and B values for different sets of model parameters.
κ
α
A
B
0.1
1
-
1.5583
-
1.6348
i
-
4.6373
+
2.0952
i
3
-
18.0859
-
2.6534
i
-
23.1309
+
3.2856
i
5
-
23.6809
-
3.3261
i
-
30.0499
+
4.0098
i
0.01
1
23.9849
-
37.1567
i
-
34.8651
+
38.3333
i
3
38.9742
-
59.1730
i
-
57.6358
+
60.8270
i
5
10.9090
-
537.4980
i
-
104.6609
+
610.5019
i
0.001
1
372.8097
-
389.5891
i
-
385.0603
+
389.7408
i
3
593.6919
-
619.3798
i
-
611.3080
+
619.6202
i
5
732.7149
-
764.5012
i
-
754.4851
+
764.7988
iIn the first example of investigation, we set the channel height2
h to be 200
λ where the λ typically has value of hundreds of nanometres. Other parameters are chosen as α
=
1 and λ
=
1. Figures 2 and 3 show the distribution of velocity along the channel height obtained from the model with various l and κ values at θ
=
π
/
2, as functions of χ for different values of the parameters κ and l. Figure 4 presents the velocity distribution U
(
χ
,
κ
) and U
(
χ
,
l
), at θ
=
π
/
2 for α
=
1 and h
=
100
λ. The results show that as κ increases, the thickness of the layer with significant electrically driven flow decreases and essentially the velocity at the centre of the channel is zero for κ
≥
0.07 as shown in Figures 3(d) and 4(a). It is also found that the slip parameter has significant influence on the flow behaviour. At κ
=
0.001, as l increases from 0 to 0.5, the flow velocity increases from 1 to about 1.685.Effect of slip length on velocity profile obtained from the model withθ
=
π
/
2, α
=
1, and h
=
100
λ for various κ values: (a) κ
=
0.001; (b) κ
=
0.01; (c) κ
=
0.03; (d) κ
=
0.07.
(a)
κ
=
0.001
(b)
κ
=
0.01
(c)
κ
=
0.03
(d)
κ
=
0.07Effect ofκ on velocity profile obtained from the model with θ
=
π
/
2, α
=
1, and h
=
100
λ for various slip lengths: (a) l
=
0; (b) l
=
0.05; (c) l
=
0.1; (d) l
=
0.5.
(a)
l
=
0
(b)
l
=
0.05
(c)
l
=
0.1
(d)
l
=
0.53D plot ofU
(
χ
,
κ
) and U
(
χ
,
l
) obtained from the model with θ
=
π
/
2, α
=
1, and h
=
100
λ.
(a)
U
(
χ
,
κ
) when l
=
0.05
(b)
U
(
χ
,
l
) when κ
=
0.03To investigate the effect of the energy parameterα on the flow behaviour, we choose four values of α including 1, 3, 5, and 7. Figure 5 shows the influence of α on the velocity distribution for θ
=
π
/
2, h
=
100
λ, l
=
0.05 and two different κ values of 0.001 and 0.03. For a channel with fixed dimension, filled with a given buffer solution, the value of κ depends directly on the frequency of the external electric field. Hence, in practice, it is possible to obtain the desired velocity field by controlling the frequency by electric switch [7]. We found that as α increases, the flow velocity increases. At κ
=
0.001, as α increases from 1 to 7, the flow velocity at the centre increases from 1 to 1.4 as shown in Figure 5. At κ
=
0.03, as α increases from 1 to 7, the flow at the centre approaches zero.Effect of the energy parameterα on the velocity profile obtained from the model with θ
=
π
/
2, h
=
100
λ, and l
=
0.05 for two different κ values: (a) κ
=
0.001 and (b) κ
=
0.03.
(a)
κ
=
0.001
(b)
κ
=
0.03At the typical values ofα
=
1, κ
=
0.03, and l
=
0.05, the distributions of velocity along the height of the channel at various instants of time in an oscillation cycle of electric field are shown in Figure 6. The magnitude of the velocity varies with time, which is because of the fluctuation of the external electric field in time.Figure 6
2D and 3D plots of velocity profile at various instants of time obtained from the model withα
=
1, κ
=
0.03
,
l
=
0.05, and h
=
100
λ.Figure7 shows the influence of channel height on the magnitude of velocity. It can be noted that the magnitude of velocity is influenced by the channel height h, while the pattern of velocity profile along the channel is the same. Figure 8 shows the influence of κ on the first derivative of the velocity.Effect of the channel half-heighth on the velocity profile obtained from the model with θ
=
π
/
2, α
=
1, and l
=
0.05 for two different κ values: (a) κ
=
0.001 and (b) κ
=
0.03.
(a)
κ
=
0.001
(b)
κ
=
0.033D plot of(
D
U
/
D
χ
)
(
χ
,
l
) for two κ values: (a) κ
=
0.001 and (b) κ
=
0.03.
(a)
(
D
U
/
D
χ
)
(
χ
,
l
) for κ
=
0.001
(b)
(
D
U
/
D
χ
)
(
χ
,
l
) for κ
=
0.03


## 5. Concluding Results

In this paper, we study the time periodic electroosmotic flow in a two-dimensional microchannel with slip boundary. Some analytical and numerical solutions for the slip velocity field under various different conditions have been obtained.Our results show that the flow behaviour depends on the parametersκ and α as well as the slip parameter l and channel height. The parameter κ represents the ratio of the electric double layer thickness to a characteristic diffusion length scale I
D, while α is an energy parameter.We found that, depending on the value of the parameterκ, the flow pattern may change from oscillatory plug flow to uniform flat flow. Generally, as κ increases, the velocity at the centre drops to zero dramatically. We also found that the slip parameter has significant influence on the flow behaviour. As l increases, the flow rate increases.



---

*Source: 789147-2014-06-16.xml*