# Study of a Newtonian Fluid through Circular Channels with Slip Boundary Taking into Account Electrokinetic Effect

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

---

## Abstract

We study the slip flow of fluids driven by the combined effect of electrical force and pressure gradient. The underlying boundary value problem is solved through the use of Fourier series expansion in time and Bessel function in space. The exact solutions and numerical investigations show that the slip length and electrical field parameters have significant effects on the velocity profile. By varying these system parameters, one can achieve smooth velocity profiles or wave form profiles with different wave amplitude and frequency. This opens the way for optimizing the flow by choosing the slip length, the electrical field, and electrolyte solutions.


---

## Body

## 1. Introduction

Over the past few decades, advances in nanoscience and nanotechnology have led to the development of many microelectromechanical systems and devices such as heat exchanger [1], micropump [2], lab-on-a-chip diagnostic devises [3], drug delivery systems [4], energy conversion, and biological sensing devices [5]. Most of these systems and devices involve fluid flow in microtubes and microchannels. To control the microfluidics in microchannels so as to achieve optional system performance, it is essential to study the fundamental mechanics of microflows and derive better models and understanding of the flow mechanism and flow behaviour.Flow of microfluidics may be driven by pressure gradient or by electrical forces. For electrical-driven flow, the solid surface of the microchannel is electrically charged to generate a region in the fluid with a distribution of electrical charges near the channel surface. This region is called the electrical double layer (EDL). The EDL, on one hand, has the effect of retarding liquid flow driven by external pressure gradient to form a streaming potential, whereas, on the other hand, it can induce fluid flow by applying an external electric field. In this work, we will study the flow of microfluidics in microchannels driven by the combination of pressure gradient and externally applied electric field.The equations governing the flow of microfluidics in microchannels include the Navier-Stokes equations, the incompressible continuity equation, and boundary conditions. Traditionally, the no-slip boundary condition is used [6, 7]. However, recent molecular dynamic simulations and experiments in micrometer scale show that the flow of fluids in micrometer scale is granular and slip occurs between the fluid and the solid surface [8–12]. Therefore, the no-slip condition does not work for fluid flow in microchannels. In this work, the Navier slip boundary condition will be used; namely, the tangential fluid velocity relative to the solid surface is proportional to the shear stress on the solid-fluid interface. The validity of the Navier slip boundary condition is supported by many experimental results [13–15].For many fluid flow problems, under the nonslip assumption, exact and numerical solutions have been obtained and can be found in the literature [6, 16–18]. Steady state solutions under slip conditions have also been established for flows of Newtonian fluids through pipes, channels, and annulus [19, 20]. More recently, various analytical solutions for pressure-driven time-dependent slip flows of Newtonian fluids through microtubes and microannulus were derived [16, 18]. Various attempts have also been made to study the electrically driven fluid flows. Rice and Whitehead [21] investigated the steady-state liquid flow due to an electric field in circular capillaries. Levine et al. [22] analysed the electrokinetic steady flow in a narrow parallel-plate microchannel. Yang et al. [23] studied the flow in rectangular microchannels and Mala et al. [24] in parallel-plate microchannels.Motivated by the previous work, this work aims to generalize the result in [25] for the pressure-driven slip flow to the case with the combined effect of pressure-driving forces and electrically driving forces. The rest of the paper is organized as follows. In Section 2, we present the mathematical model for the problem, consisting of the governing field equations and boundary conditions. In Section 3, we derive the exact solutions for the velocity field through the use of Bessel functions in space and Fourier series expansion in time. Based on the velocity solutions, we then establish the analytical solutions for the stress tensor and flow rate. In Section 5, we investigate the influence of the electric field and the surface slippage on the flow behaviour. Then a conclusion is given in Section 6.


## 2. Mathematical Model and Formulation

In this paper, we consider the flow of microfludics through a circular microchannel driven by both pressure gradient and external electric field. The cylindrical polar coordinate (r, θ, z), with the z-axis being in the axial direction, is used in the formulation. The governing field equations for the problem include the Navier-Stokes equations and the continuity equation. Let v_=(vr,vθ,u) be the velocity vector with vr,vθ, and u being, respectively, the components of velocity in the radial direction, the transverse direction, and axial direction. Then the continuity equation and the Navier-Stokes equation in the z-direction are
(1)∂2vr∂r2+vrr+∂vθr∂θ+∂u∂z=0,∂u∂t+vr∂u∂r+vθr∂u∂θ+u∂u∂θ=gz-1ρf∂p∂z+μρf(∂2u∂r2+∂ur∂r+1r2∂2u∂θ2+∂2u∂z2)+ρeρfEz′,
where p is the pressure, μ is the dynamic viscosity, ρf and ρe, are respectively, the fluid density and the electric charge density, and Ez′ is the externally applied electric field.Assuming that the flow is axially symmetric and the radial and transverse velocity components are negligible, then (1) admits solution of the form
(2)vr=0,vθ=0,u=u(r,t)
with u(r,t) determined by
(3)∂u∂t-μρf(∂2u∂r2+1r∂u∂r)=-1ρf∂p∂z-ρeρfEz′→.
Based on [25], the velocity boundary condition for the problem can be written as
(4)∂u∂r(0,t)=0,u(R,t)=-l∂u∂r(R,t).
In this paper, we study the flow driven by both pressure gradient ∂p/∂z and an externally applied electrical field Ez′. Without loss of generality, we express the pressure gradient by Fourier series; namely,
(5)∂p∂z=q¯(t)=a0+∑n=1∞‍[ancos(nωt)+bnsin(nωt)],
which can also be expressed by
(6)∂p∂z=Re(∑n=0∞‍cneinωt),
where Re() is to take the real part of the complex quantity. To determine the free charge density ρe, let ψ(r) be the electric potential associated with double layer at the equilibrium state. Then based on [26], we have
(7)ρe=-κ2ψ,(8)d2ψdr2+1rdψdr=κ2ψ,(9)ψ(R)=ψs,ψ(0)isfinite,
where ψs is the surface potential at the wall r=R. Now substituting (7) into (3), we have
(10)∂u∂t-μρf(∂2u∂r2+1r∂u∂r)=-1ρf∂p∂z+1ρfκ2ψEz′.
Hence the problem is to solve the boundary value problem (8)-(9) for ψ and then the partial differential equation (10) subject to boundary conditions (4) for u, which will be done in Section 3.


## 3. Solution for the Transient Velocity Field

Letκr=r¯; (8) becomes
(11)r¯2d2ψdr¯2+r¯dψdr¯-r¯2ψ=0
which gives solution
(12)ψ=AJ0(κr)+BY0(κr),
where A and B denote integration constants and J0 and Y0 are the zero-order Bessel functions of the first and second kinds, respectively. From (12) and boundary condition (9), we obtain
(13)ψ=ψsJ0(κR)J0(κr)
and (10) becomes
(14)∂u∂t-μρf(∂2u∂r2+1r∂u∂r)=-1ρfRe∑n=1∞‍cneinωt+1ρfκ2ψsJ0(κR)J0(κr)Ez′
which admits solution of the form
(15)u(r,t)=up(r,t)+ue(r,t),
where up is the solution corresponding to the first term on the right hand side and ue is the solution contributed from the second term of the right hand side. From the superposition principle, based on the work of [25], we have
(16)up(r,t)=C+a04μr2+∑n=1∞‍Re{(dnJ0(βnr)+cniρfnω)einωt},
where βn2=nβ2 with β2=(ρfω/μ)i. Now, we proceed to find the solution contributed from the second term on the right hand side of (15). Let ue(r)=C0(r)t+C1(r) and substitute it into (14); then we have
(17)-μρf[C0′′(r)+1rC0′(r)]t+C0(r)-μρf[C1′′(r)+1rC1′(r)]=1ρfκ2ψsJ0(κR)Ez′J0(κr).
For the above equation to hold for any instant of time t, we require
(18)C0′′(r)+1rC0′(r)=0,C0(r)-μρf[C1′′(r)+1rC1′(r)]=1ρfκ2ψsJ0(κR)Ez′J0(κr).
Thus,
(19)C0=D1lnr+D2.
As C0 must be bounded at r=0 but lnr has singularity at r=0, we require that D1=0. From (18), we have
(20)C1=1μψsEz′J0(κR)J0(κr)+ρfμD24r2,
and then by substituting C0 and C1 into ue(r,t), we have
(21)ue(r,t)=D2t+1μψsEz′J0(κR)J0(κr)+ρfμD24r2,
Hence, by substituting (16) and (21) into (15), we obtain
(22)u(r,t)=C+a04μr2+∑n=1∞‍Re{(dnJ0(βnr)+cniρfnω)einωt}+D2t+1μψsEz′J0(κR)J0(κr)+ρfμD24r2.Substituting (22) into boundary condition (4) yields
(23)C+a04μR2+la02μR+1μψsEz′-lκψsEz′μJ0(κR)J1(κR)+D2t+ρfD2R24μ+lρfD2R2μ+∑n=1∞‍Re[(dnJ0(βR)+cniρfnω-lβndnJ1(βnR))einωt]=0.
For the above equation to hold for any instant of time t, we require D2=0 and
(24)C+a04μR2+la02μR+1μψsEz′-lκμψsEz′J0(κR)J1(κR)=0,dnJ0(βR)+cniρfnω-lβndnJ1(βnR)=0,
which give
(25)C=-a0R24μ(1+2lR)-1μψsEz′+lκμψsEz′J0(κR)J1(κR),dn=-cniρfnω[J0(βR)-lβnJ1(βnR)].
Hence, by substituting (25) into (22), we obtain
(26)u=-a0R24μ(1+2lR-r2R2)-1μψsEz′+lκμψsEz′J0(κR)J1(κR)-∑n=1∞‍Re{cniρfnω(J0(βnr)J0(βnR)-lβnJ1(βnR)-1)einωt}+1μψsEz′J0(κR)J0(κr).
We should remark that if l=0, solution (26) reduces to the no-slip solution
(27)u=-a0R24μ(1-r2R2)-1μψsEz′-∑n=1∞‍Re{cniρfnω(J0(βnr)J0(βnR)-1)einωt}+1μψsEz′J0(κR)J0(κr).


## 4. Solutions for the Flow Rate and Stresses

From the velocity solution (26), we get the flow rate as follows:
(28)Q(t)=∫0R2πrudr=-a0πR32μ(l+R4)+lπR2κμψsEz′J0(κR)J1(κR)+∫0R1μψsEz′J0(κR)J0(κr)rdr-2πρfω×Re[{-1J0(βnr)J0(βnR)-lJ1(βnR))rd}∑n=1∞cnieinωtn∫0R(J0(βnr)J0(βnR)-lJ1(βnR)-1J0(βnr)J0(βnR)-lJ1(βnR))rdr∑n=1∞cnieinωtn]-πR2μψsEz′.
From the identity
(29)ddx[xnJn(x)]=xnJn-1(x),
we have
(30)ddr[rJ1(βnr)]=βnrJ0(βnr),
and hence
(31)∫0RrJ0(βnr)dr=1βn[rJ1(βnr)]0R=1βnRJ1(βnR).
Therefore from (28) and the above formula, we get
(32)Q(t)=-a0πR32μ(l+R4)+lπR2κμψsEz′J0(κR)J1(κR)+RμκψsEz′J0(κR)J1(κR)-2πρfω×Re[∑n=1∞‍cnieinωtn(RJ1(βnR)βnJ0(βnR)-lβnJ1(βnR)-R22)∑n=1∞‍cnieinωtn(RJ1(βnR)βn[J0(βnR)-lβnJ1(βnR)]-πR2μψsEz′.
The total amount of fluid, flowing through the tube during the period [0,T], can then be determined as follows:
(33)QT=∫0TQ(t)dt=-a0πR3T2μ(l+R4)+lπR2κμψsEz′TJ0(κR)J1(κR)+πRμκψsEz′TJ0(κR)J1(κR)-2πρfω2×Re[(RJ1(βnR)βn[J0(βnR)-lβnJ1(βnR)-R22)∑n=1∞‍cn(einωT-1)n2×(RJ1(βnR)βnJ0(βnR)-lβnJ1(βnR)-R22)∑n=1∞‍cn(einωT-1)n2]-πR2μψsEz′T.Now, we determine the stresses in the fluid. Fromv=er0+eθ0+ezu(r), we get
(34)∇v=(00∂u∂r000∂u∂r00).
From the above formulae and (26), we obtain
(35)drθ=dθz=drr=dθθ=dzz,drz=a0r2μ+∑n=1∞‍Re[cniρfnω(βnJ1(βnr)J0(βnR)-lβnJ1(βnR))einωt]-κμψsEz′J0(κR)J1(κr).
Thus from the constitutive equations for Newtonian fluid, we get
(36)σrr=σθθ=σzz=-p=q¯(t)x+p0(t),σrθ=σθz=0,σrz=a0r+2μ∑n=1∞‍Re[cniρfnω(βnJ1(βnr)J0(βnR)-lβnJ1(βnR))einωt]-2μκψsEz′J0(κR)J1(κr),
where p0(t) is an arbitrary constant which may be chosen to meet certain pressure conditions.


## 5. Numerical Investigation

In this section, we use the solutions obtained in previous sections to study the influence of the combined effect of electrical field and pressure gradient on the velocity profile and flow rate. As the pressure gradient can be expressed as a Fourier series, without loss of generality, we consider here only two cases of pressure gradient, including the constant pressure gradient and a cosine wave form pressure gradient.Case 1.
Consider(37)dpdz=a0.
In this case, we havecn=0 for all n≥1, and thus, from (26) and (32), we get the normalized velocity and flow rate as follows:
(38)u*=-4μa0R2u=(1+2lR-r2R2)-4μa0R2[lκμJ1(κR)J0(κR)+1μJ0(κr)J0(κR)-1μ]ψsEz′,Q*(t)=-2μa0πR3Q(t)=(l+R4)-2μa0πR2×[lπκRμJ1(κR)J0(κR)+1μκJ1(κR)J0(κR)-πRμ]ψsEz′.
We should address here that the velocity and flow rate are linearly propositional to the slip length l.Case 2.
Consider(39)dpdz=a1cos(ωt).
In this case, we havea0=0,c1=a1∈R,cn=0 for all n≥2. For simplification, we normalize the variables as follows:
(40)β*=βR,r*=rR∈[0,1],l*=lR,t*=ωt2π,u*=-ρfωa1u,K=κR.
Then from (26), we get
(41)u*=Re((J0(β*r*)J0(β*)-l*β*J1(β*)-1)ie2πit*)-ρfωa1[1μJ0(Kr*)J0(K)+κRl*μJ1(K)J0(K)-1μ]ψsEz′.Let(42)M=M(r*)=Re(J0(β*r*)J0(β*)-l*β*J1(β*)),N=N(r*)=Im(J0(β*r*)J0(β*)-l*β*J1(β*)),A1=-ρfωμJ0(K),A2=ρfωKμJ0(K)J1(K)R,A3=ρfωμ,
then
(43)u*=A(r*)cos(2πt*-θ)+(A1J0(Kr*)+A2l*+A3)ψsEz′a1,
where
(44)A=N2+(1-M)2,θ=arctan(M-1N).
As the M and N defined above are in terms of the complex parameter β* and the Bessel functions with complex arguments, the relationship between the velocity and the slip length l needs to be investigated. We first proceed to derive more explicit formulae relating u* and l in the real domain. As β2=-(ρω/μ)i=-(ρω/μ)e-π/2, we get
(45)β=ρω2μ(1-i)=β¯R(1-i),1β=R2β¯(1+i),β*=β¯(1-i),
where β¯=Rρω/2μ. For ∥y∥≪1, then we have the following asymptotic formulae for approximating Bessel functions [12]:
(46)Jn(x+yi)≈Jn(x)-iy2[Jn+1(x)-Jn-1(x)],Jn(y)≃1n!(y2)2.
Thus, for β¯≪1, we have
(47)J0(β*)=J0(β¯-β¯i)≈J0(β¯)+iβ¯J1(β¯)≈1+β¯22i,J1(β*)=J1(β¯-β¯i)≈β¯2-β¯2[1-β¯28]i,J0(β*r*)≈1+β¯2r*22i,J1(β*r*)=β¯r*2-β¯r*2[1-(β¯r*)28]i.Substituting the above into (48) yields
(48)M≈1-β¯4l*(1-8r*2)/161+β¯4l*((7/8)+l*),N≈-β¯2(1-r*2+2l*)2[1+β¯4l*((7/8)+l*)].
Hence we have
(49)A=A(r*)=β¯22[1+β¯4l*((7/8)+l*)]×(1-r*2+2l*)2+β¯4l*264(15-8r2+16l*)2,u*=A(r*)cos(2πt*-θ)+[A1J0(Kr*)+A2l*+A3]ψsEz′a1.Thus, it can be concluded that the velocity consists of a time-dependent part in terms of a trigonometric function and a time-independent part in terms of a Bessel function of the first kind of order zero.To demonstrate the velocity profile of the mixed pressure-driving and electrical-driving flow and the influence of the slip parameter and electric field on the flow, we carry out numerical investigation under various conditions. The typical values chosen for the model parametersR,ρf,ω,μ are R=20 mm, ρf=1000 km/m3, ω=1/2π,μ=0.86 Ns/m2. Based on [27–29], the K value varies from one electrolyte solution to another; thus in this work we use a wide range of values to demonstrate the influence of K on the flow. The slip length l is related to the smoothness of the surface of the tubes and is controllable within certain range; hence we also use a wide range of l values in the investigation to demonstrate the influence of l on the flow profile and flow rate. The ψsEz′ is set to 0.05 V in all experiments except for the experiment in which the influence of the intensity of electric field relative to the pressure-driving force is to be investigated.As the velocity profile consists of two parts due, respectively, to pressure gradient and electrical field, we first show the profile of these two velocity components in Figures1 and 2 under different electric fields ψsEz′/a1 and K values. From the figure, it is clear that it is possible to control the velocity profile in an optimal manner based on the need of application for the mixed driving flow by finding proper combination of the pressure field and electrical field. Figures 3 and 4 show typical velocity profile for mixed driving flows.Figure 1
Typical velocity profile of the velocity component due to pressure gradient under differentl* value.Figure 2
Typical profiles of the velocity component due to electric field forl*=0.5 under different ψsEz′/a1 and K values.Figure 3
Typical velocity profile for various different values of slip lengthl* at a typical time t*=θ/2π, obtained with ψsEz′/a1=10-4, K=200.Figure 4
Typical velocity profile for the mixed pressure and electric driven flow at various instants of timet*, obtained with slip length l*=1, ψsEz′/a1=10-6, and K=1500.To demonstrate the influence of the slip length on the velocity, we vary the value ofl* from 0 to 100 while maintaining other parameter values unchanged. Figure 5 shows the effect of slip length on the velocity on the wall at a typical time step t*=θ/2π, obtained with ψsEz′/a1=10-6, K=3000, and K=2000. It is clear that the relationship between the slip parameter and the velocity is nonlinear and there is an optimal value for which the slip velocity takes the maximum value.Figure 5
The influence of slip length on the slip velocity on the wall at a typical timet*=θ/2π, obtained with ψsEz′/a1=10-3, l*=1.As theK value changes across different electrolyte solutions, we investigate the influence of K on the velocity profile. Figure 6 shows the velocity profiles for different K values for ψsEz′/a1=10-3 and under the same pressure gradient and slip length l*=1. It is clear from the result that the K value influences the velocity profile significantly. For lower K, the velocity profile is smooth, while under high K values, the velocity shows wave form profile and the frequency of the waves increases as K increases. It is also noted that the amplitude of the waves increases as the ψsEz′/a1 value increases. The magnitude of the velocity is also affected significantly by the ψsEz′/a1 value as shown in Figure 7. As ψsEz′/a1 increases, the magnitude of the velocity increases.Figure 6
Influence ofK value on the velocity profile at a typical time t*=θ/2π, obtained with ψsEz′/a1=10-3,l*=1.Figure 7
Influence of magnitude of electric field on the velocity profile at a typical timet*=θ/2π, l*=0.5, K=200.


## 6. Conclusions

In this paper, an exact solution for the transient flow of an incompressible Newtonian fluid in microtubes is derived taking into account the electrokinetic effect and the boundary slip. We have shown that both boundary slip and electrokinetic field have significant effect on the flow. The results are summarized as follows.(i) Both the pressure-gradient driving force and electrical-driving force influence the velocity profile significantly, and it is possible to construct an optimal control problem to achieve the desired velocity profile based on the need of the application.(ii) For the case of constant pressure gradient, the relationship between the slip parameter and the velocity is linear. However, for the pressure gradient in wave form, the influence of the slip parameter on the velocity magnitude is nonlinear and there exists a critical slip length where the magnitude and consequently the transient flow rate attain maximum values, which indicates that there exists an optimal slip length in terms of the magnitude of velocity and flow rate.(iii) TheK value of the electrolyte solution has significant influence on the velocity profile on the cross-section of the tube. At low K value, no wave form profile exists. However, as K increases, the cross-section velocity exhibits wave form profile, and the frequency of the wave increases as K increases.(iv) The magnitude of the intensity of electric field relative to the pressure gradient, measured byψsEz′/a1, also has very significant influence on the velocity profile. As ψsEz′/a1 increases, the velocity component increases and the amplitude of the wave in the waved form velocity profile also increases.



---

*Source: 718603-2013-11-28.xml*