Electro-osmotic flow of a third-grade fluid past a channel having stretching walls

1 Introduction

The electroosmotic flow (EOF) is the flow induced by the application of an electric field across a channel and due to the presence of electric double layer (EDL) at channel walls. Electric double layers exist in all heterogeneous fluid based systems including biofluids. The electroosmotic flow is most significant in small channels where the EDL is much smaller than the characteristic length scale of the channel. The electroosmotic flow is considered to be an inseparable part of many scientific and engineering applications used in microfluidic devices. The various applications of microfluidic devices include micro-reactors, heat transfer in electronic circuits, biosensors, drug delivery, DNA analysis and sequencing, etc. It has been projected that microfluidic devices utilizing EOF will have great application in medical research [1, 2]. As the biofluids are electrically conducting in nature, EOF will be vital for drug delivery and separation/mixing of fluids at the atomic level. Studies on EOF of Newtonian fluids abound in contrast to non-Newtonian fluids as evident from the literature. In the recent years, to handle complex fluids in microfluidics, EOF of non-Newtonian fluids has gained attention.

The study of EOF was first reported by F.F. Reuss [3] in 1809 where he has showed that water can flow through a plug of clay by the application of suitable electric voltage. Most of the early studies [4,5,6] have focused on the theoretical analysis of electroosmotic flow of Newtonian fluids. Santiago [7] analyzed the effects of fluid inertia and pressure on the velocity and vorticity field of electroosmotic flow of Newtonian fluids. The analytical solution for electrostatic potential of double layer has been obtained by solving the complete Poisson-Boltzmann equation for arbitrary zeta potentials using Green’s function method for Newtonian fluids in [8]. A theoretical model was developed by Chakraborty [9] for studying the capillary filling dynamics of a power-law fluid in a microchannel. Park and Lee [10] numerically obtained flow field solution of full Phan-Thien and Tanner (PTT) constitutive equation in a rectangular duct under the action of an external electric field and applied pressure gradient by employing finite volume method. Berli and Olivares [11] introduced a theoretical description of the electrokinetic flow of power law fluid through slit and cylindrical microchannels. Exact solution of the velocity distribution are obtained for the EOF of power-law fluids in a slit channel analytically by Zhao et al. [12]. Siddiqui and Lakhatia [13] fomulated and solved the boundary value problem of steady, symmetric one-dimensional electroosmotic flow of a micropolar fluid in a uniform rectangular microchannel under the action of uniformly applied electric field. A numerical study of electroosmotic flow in microchannels for non-Newtonian fluid obeying power-law model was attempted for the first time in [14]. Analytical solutions for channel and pipe flows of visco-elastic fluids obeying PTT model and the FENE-P model, under the mixed influence of electrokinetic and pressure forces were obtained in [15]. Vasu and De [16] have numerically analyzed the fluid flow behavior for EOF of power-law fluids in a rectangular microchannel at high zeta potential by considering the EDL potential distribution without Debye Huckel linear approximation. A numerical investigation of flow of power-law fluid in a two dimensional microchannel with non-uniform zeta potential distributions along the channel walls by finite volume scheme was carried out by Hadigol et al. [17]. Analytical solutions were obtained for transverse distribution of velocity, and thermal transport characteristics of the steady fully developed electroosmotic flow of PTT and FENE-P model fluid in a slit microchannel in [18]. Semi-analytical solutions for transient EOF of Maxwell fluids between micro-parallel plates are presented in [19] using Laplace transform method. A numerical analysis is performed in [20] on alternating current EOF of non-Newtonian incompressible fluids based on the Ostwald-de Waele power-law model in circular microchannels. The electroosmotic flow of Maxwell fluids in a rectangular microchannel with asymmetric high zeta potentials at the walls is studied analytically using Green’s functions and separation of variables method be Jimenez et al. [21]. Xie and Jian [22] analyzed the rotating electroosmotic flow of power-law fluids at high zeta potentials in a slit microchannel using finite difference method. Electroosmotic flow through a channel between two undulated surfaces induced by an external electric field is studied in [23] using lubrication approximation theory and lattice Boltzmann method. The EOF of power-law fluids in microchannels with different geometries is discussed in [24, 25].

A theoretical study on the electroosmotic flow of a fluid past a channel having stretching walls has been investigated by Misra et al. [26] by considering the biofluid to be a second grade fluid. Although the second grade fluid model can predict the characteristics of biofluids, it does not take into account the shear thinning phenomenon which is a vital property of biofluids. This motivated us to extend their work and analyze the electroosmotic flow of a third grade fluid past a channel having stretching walls in the present work. The partial differential equation governing the flow is transformed to a nonlinear ordinary differential equation using appropriate similarity transformation which is then solved using a suitable finite difference scheme. The results computed are presented in the form of graphs for various flow parameters governing the flow.

2 Basic Equations

The basic laws of the conservation of mass, conservation of momentum, and conservation of energy for an incompressible fluid are

where V is the velocity field, ρ is the fluid density, C_p is the specific heat at constant pressure, κ is the coefficient of thermal conductivity, F is the body force per unit volume, T is the temperature, and τ is the stress tensor. The constitutive equation of Cauchy stress tensor for a third grade fluid is

where p is the pressure, μ is the coefficient of viscosity, α₁, α₂, β₁, β₂, β₃ are material constants and the Rivlin-Ericksen tensors A₁, A₂, and A₃ are given by

Fosdick and Rajagopal [27] have shown that the motion satisfies Clausius-Duhem inequality to be compatible with thermodynamics and if it is assumed that the specific Helmholtz free energy is minimum when the fluid is at rest then the following constraints must be satisfied

The F term in equation (2) is given as

where E is the applied electric field and ρ_e is the net electric charge density. For a symmetric electrolyte the anions and cations have same charge valence z⁺ = z⁻ = z and the net electric charge density ρ_e = ez(n⁺ − n⁻) where n⁺ and n⁻ are the concentrations of positive and negative ions, respectively and e is the elementary electric charge.

The electric field intensity E = −∇Φ, where Φ is the electric potential governed by

where ε is the electric permittivity of the solution.

In electroosmotic flows, two types of electric fields act on the flow. One is the applied electric field ϕ which is generated by the electrodes at the inlet and outlet of the channel and the other one ψ is due to the net charge distribution in the electric double layer which is acquired at the wall. The total electric field Φ is a linear superposition of ϕ and ψ.

At the centerline of the channel the EDL does not overlap and variation of the electrical potential distribution ψ occurs only in the normal direction to the channel walls which can be expressed in the y-direction as

According to Bruss [28], the net electric charge density in the EDL can be expressed as

where n₀ is the bulk number concentration of ions in the electrolyte solution, k_B is the Boltzmann constant and T_A is the absolute temperature. The Debye-Hückel linear approximation yields

3 Problem Formulation

We consider the fully developed electroosmotic flow of a viscous incompressible fluid, that obeys the third grade fluid model, in a channel bounded by stretching walls with equations y = ± R₀ (Fig.1). The flow is assumed to be symmetric about the centerline and hence we confine our study to the region 0 ≤ y ≤ R₀. The height of the channel is much smaller than the length of the channel, L. The effect of Joule heating, gravity, and pressure gradient have been neglected in the present study. The flow is driven by stretching of the channel walls and the electro-osmotic effect. The channel walls in contact with the electrolyte develop surface charge due to ion adsorption which results in the formation of the electric double layer. The physical extent of the EDL is represented by the Debye length. The EDL is very thin, therefore the effective fluid velocity is quantified by the Helmholtz-Smoluchowski equation

Fig. 1

Sketch of the problem

where ζ is the zeta potential. In the boundary layer, the equations that govern the flow for the present problem are [29]

and

The electric potential distribution for a symmetric electrolyte described by the Poisson-Boltzmann equation in Eq.(5) and Eq.(6) is linearized using Eq.(7), and is given by

where λD=εkBTA2n0z2e2 is the Debye length.

The boundary conditions for the velocities and electric potential for the flow problem are

Using the following similarity transformation

and introducing the non-dimensional variables

the partial differential equation (10) transforms to the following ordinary differential equation

where

The boundary conditions become

On solving the linearized Poisson-Boltzmann equation with the boundary conditions, we get

Eq. (8) now becomes

4 Numerical Methods

The resulting equation arising from the problem of EOF of a third grade fluid past a channel having stretching walls is highly non-linear in nature and does not admit an exact analytical solution. We, therefore, look for a numerical solution of the problem under consideration. For simplicity, better stability, accuracy and efficiency, the finite difference technique is used for discretization. In order to discretize the differential quotients the following suitable finite difference strategy is opted.

We have divided the domain [0, 1] into a set of non-overlapping grids such that the spacing of the grid points, Δη is uniform. The grid points are identified by an index i running in the η-direction.

At the initial point (i = 0), Euler’s forward finite difference scheme and at the end point (i = N), the backward finite difference is used.

For i = 1, the following difference scheme is used

The central difference scheme is used for discretization of the terms at the nodes i = 2, 3, …, N − 2.

For i = N − 1, the following difference scheme is used

The resulting algebraic system of equations is then expressed in the residual form as

where the residuals are given as follows

For m = 2, 3, …, N-2

For computation of the solution of non-linear system of algebraic equations (17), Newton’s iterative method that converges quadratically, is used. In order to stabilize the convergence at an early stage of the iteration and to save computing time the damped Newton’s method [30] is used which is given by

where f = [f₁, f₂, …, f_N]^T represents the column vector of unknowns, 0 < λ^k < 1 is the k^th damping parameter, and J(f^k) is the Jacobian matrix evaluated at the k^th iterate. Here the Jacobian matrix J(f) is

A good initial guess is necessary for the convergence of the damped Newton’s method and very important for the fast solution of the iterative process. The initial solution guess for the axial velocity is approximated according to the boundary conditions. For this problem the convergence of the damped Newton method is achieved when

5 Results and Discussion

The result of our investigation is discussed in this section with graphical presentation of f′(η) and f(η) which correspond to the axial velocity and vertical velocity, for various governing flow parameters. The effect of the visco-elastic parameter α, the cross viscous parameter y, the third grade material parameter β, the Reynold’s number Re and the Debye Hückel parameter K on the velocity profiles is studied through numerical comparisons among the relevant curves. Although computation for several combinations of flow parameters were performed, only some key graphs are provided in this paper.

Figures 2-6 illustrate the shape of the dimensionless axial velocity profile for different flow parameters. It is depicted from Fig. 2 that the axial velocity increases for increasing viscoelastic parameter, α. The velocity profile maintains a constant value upto a certain distance from the centreline and thereafter decreases monotonically as approaches the wall. Similar behavior is observed for the case of varying β, the third grade material parameter and y, the cross viscous parameter as shown in Fig. 3 and Fig. 4, respectively. Variation of the axial velocity for various Re is provided in Fig. 5. The velocity increases for increasing value of Re upto η = 0.6 and for η > 0.6 the behavior of the velocity gets reversed for increasing Re. Further Re does not have much effect on f′ near the boundary for Re > 5.0. Fig. 6 shows the variation of horizontal component of velocity for different values of the Debye-Hückel parameter, K. As K increases f′ increases but a deviation is observed very close to the channel wall which is displayed in the inset figure.

Fig. 2

Variation of f′ with α at β = 0.005, y = 0.01, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 3

Variation of f′ with β at α = 0.01, y = 0.01, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 4

Variation of f′ with y at α = 0.01, β = 0.005, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 5

Variation of f′ with Re at α = 0.01, β = 0.005, y = 0.01, K = 10.0 and Λ = 0.001.

Fig. 6

Variation of f′ with K at α = 0.01, β = 0.005, y = 0.01, Re = 10.0 and Λ = 0.001.

Figures 7-11 represent vertical velocity profiles for different flow parameters. It appears to increase from zero on the axis while moving away from it and finally approaches zero again at the boundary. The velocity increases with increasing values of α, β, y, K and Re. The figures show that these flow parameters do not affect the velocity near the channel wall. The shape of the vertical velocity curve changes for Re > 5.0. It can be estimated quantitavely that the magnitude of the vertical velocity component is lesser than the axial velocity component by comparing the results presented in both flow cases.

Fig. 7

Variation of f with α at β = 0.005, y = 0.01, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 8

Variation of f with β at α = 0.01, y = 0.01, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 9

Variation of f with y at α = 0.01, β = 0.005, Re = 10.0, K = 10.0 and Λ = 0.001.

Fig. 10

Variation of f with K at α = 0.01, β = 0.005, y = 0.01, Re = 10.0 and Λ = 0.001.

Fig. 11

Variation of f with Re at α = 0.01, β = 0.005, y = 0.01, K = 10.0 and Λ = 0.001.

The figures presenting the velocity profile reveal the fact that velocity increases with the rise in non-Newtonian parameters α, y and β which is physically relevant. As the visco-elastic parameter α and cross viscous parameter y increases, the shear rate increases leading to increase in the deformation, rearrangement and aggregation of fluid particles, thereby increasing the fluid velocity in the channel. Similarly, with increasing shear thinning parameter β the viscosity of the fluid decreases, resulting in the increase of velocity. The figures associated with the Reynold’s number Re show that the inertial forces dominate the flow in the middle of the channel whereas near the wall viscous forces dominate and this is due to the presence of the electric double layer. The Debye Hückel parameter is directly proportional to the width of the channel and inversely proportional to the Debye length, which describes the physical extent of the EDL. It is observed that the velocity increases with increasing width of the channel but simultaneously, the effect of EDL is noticed at the walls of the channel.

The potential distribution function ψ has been plotted against η for different values of K in Fig. 12. The value of ψ is zero at the middle of the channel and for K > 20.0 it increases sharply to attain maximum on the boundary. It can also be inferred that as K increases, the potential distribution function decreases. For K > 100.0 variation among the curves becomes negligible. This result depicts the influence of Debye length on the potential distribution. As the EDL thickness reduces, the potential distribution function decreases and beyond a certain λ_D the potential distribution function ψ starts behaving similarly with negligible variation.

Fig. 12

Variation of electric potential distribution ψ with K.

In addition, we have computed the local skin-friction coefficient as

which in terms of dimensionless quantities become

where ReL=bR02ρμ is the local Reynold’s number. It is observed that the skin-friction coefficient increases with increasing viscoelastic parameter, α, and increases linearly with increasing shear-thinning parameter, β from Fig. 13 and Fig. 14, respectively. Further, b, the stretching parameter is directly proportional to the local Reynold’s number and it is evident from the figures that increase in the stretching parameter results in increasing value of coefficient of skin friction.

Fig. 13

Variation of C_f with α for β = 0.005.

Fig. 14

Variation of C_f with β for α = 0.001.

6 Conclusion

In the present paper, the electro-osmotic flow of a third grade fluid in a channel having stretching walls is investigated with a view to apply it in the study of blood flow dynamics. The non-linear differential equation is discretized by finite difference method and the resulting system of non-linear algebraic equations is solved by the damped Newton’s method. The results computed are presented graphically for various parameters of interest. The results obtained from this study are summarized below

1. Axial velocity profiles are found to be diminishing from their maxima at the middle of the channel to achieve the minimum value on the wall surface.

2. Axial velocity and vertical velocity increases with increasing visco-elastic parameter, α, cross viscous parameter, y and increasing third grade fluid parameter, β.

3. Axial velocity increases with the increasing Debye-Hückel parameter, K, but near the wall reverse of it is observed.

4. Potential distribution function, ψ is minimum at the middle of the channel and gradually increases as it approaches the wall. It decreases with increasing Debye-Hückel parameter, K.

5. Value of local skin-friction coefficient, C_f increases with increasing values of visco-elastic parameter, α, third grade fluid material parameter, β, and the stretching parameter b.