# Unsteady Electroosmotic Flow of Jeffrey Fluid in a Circular Microchannel under the Combined Action of Vertical Magnetic Field, External Electric Field, and Pressure at High Zeta Potential

Authors: Meirong Ren; Tiange Zhang; Jifeng Cui; Xiaogang Chen; Bixia Wu
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

The unsteady electroosmotic flow (EOF) for one kind of linear viscoelastic fluid, which is Jeffrey type fluid, is investigated under the common impact of vertical magnetic field, external electric field, and pressure at high Zeta potential in a circular microchannel. The numerical solutions of the potential and velocity distributions are obtained by solving the nonlinear Poisson-Boltzmann equation, the constitutive equation of the Jeffrey fluid, and the Cauchy momentum equation applying the Chebyshev spectral method and the finite difference method. By contrast, the Chebyshev spectral method has higher accuracy and less computation. The flow characteristics of Jeffrey fluid at high Zeta potential are analyzed with the numerical solution obtained by the Chebyshev spectral method. The results show that the velocity of Jeffrey fluid increases with the increase of the wall Zeta potential and electric width. The oscillation amplitude of velocity distribution increases with the increase of relaxation time but decreases with the increase of retardation time. When the Hartmann number is smaller, the increase of relaxation time leads to the increase of velocity; when the Hartmann number is larger, the increase of relaxation time leads to the decrease of velocity. No matter what the Hartmann number is, the velocity always decreases with the increase of the retardation time. The velocity tends to be stable gradually with the increase of time.


---

## Body

## 1. Introduction

The research on micro-electromechanical systems (MEMS), which began in the 1980s and 1990s, is recently an investigative hotspot in the world [1]. Although there are multifarious microfluidic driving and controlling technologies based on different theories, the micropumps based on electroosmotic mechanism are still the dominant micropump in MEMS. They are widely used not only in MEMS but also in drug delivery of medical treatment and aerospace technology. Therefore, scholars around the world have done a lot of researches on the electroosmotic flow of micro-nanofluids in micropipes with diverse geometric shapes [2–6]. However, under the action of external electric field, Joule heat will be generated inevitably with the increase of voltage. In recent years, in order to cut down the Joule heat effect efficaciously, people have studied the transport process of micro-nanofluids by applying a transverse magnetic field, hoping to elevate the flow velocity of micro-nanofluids through the effect of electromagnetic hydrodynamics [7–9]. For example, the effect of MHD and electroosmosis on the radiative tangent hyperbolic nanofluid flow through a porous medium is investigated by Ramesh et al. [10]. Yang et al. [11] considered theoretically the heat transfer traits of incomprehensible magnetohydrodynamic electroosmotic flow in a two-dimensional rectangular microchannel, acquired the analytical solution of dimensionless velocity, and deduced the dimensionless temperature and the Nusselt number. Sridhar and Ramesh [12] analyzed the couple stress nanofluid (blood–graphene/diamond) flow in an asymmetric channel with the effect of viscous dissipation, electromagnetohydrodynamics (EMHD), Joule heating, velocity slip and convective boundary conditions, and so forth.In most practical applications, the microbial fluid equipment is mainly dealing with biological fluid, for example, polymer solution, body fluid, and blood. These fluids are composed of large and long-chain molecules, which exhibit distinct non-Newtonian fluid characteristics, such as variable viscosity, memory effect, and normal stress effect. Jeffrey fluid model is a simple linear viscoelastic model among non-Newtonian models that can be used to simulate these biological fluids, interpret experimental data on linear viscoelasticity, and investigate the molecular structure of the polymer via linear viscoelastic measurements(such as the complex viscosity). The constitutive equation of Jeffrey fluid model is established on the basis of Maxwell fluid constitutive equation, which is a simple linear viscoelastic model in non-Newtonian model, and the second-grade fluid is a nonlinear viscoelastic model. The second-order fluid model is based on the Jeffrey fluid model, adding some additional nonlinear terms to better measure the relationship between stress tensor and strain rate tensor. The linear viscoelastic models can describe the complex viscosity and the transmission of small-amplitude shearing waves. It can also be shown that the linear models can describe elastic recoil, although the results are restricted to flows with small displacement gradients. An important reason for studying Jeffrey fluid model is that some background in linear viscoelasticity helps us in the study of nonlinear viscoelasticity, where flows with large displacement gradients are discussed. This model was suggested by H. Jeffreys, The Earth, Cambridge University Press, 1st edition (1924) and 2nd edition (1929), p. 265, to describe the propagation of waves in the earth’s mantle [13]. Because of its wide application in biology, industry, and other fields, the study of Jeffrey fluid has aroused great interest of many scholars. For instance, Ramesh [14] investigated and discussed the Poiseuille flow of an incompressible magnetohydrodynamic Jeffrey fluid between parallel plates through homogeneous porous medium using slip boundary conditions under the effect of heat transfer. Sridhar and Ramesh [15] dealt with the magnetohydrodynamic and electroosmotic radiative peristaltic Jeffrey nanofluid flow with the effect of slip and convective boundary conditions in the nonsymmetric vertical channel. Yang and Jian discussed the electromagnetic flow of Jeffrey fluid in a parallel microchannel [16]. Yasmeen et al. [17] talked over peristalsis in axisymmetric tube for Jeffrey fluid in a strong magnetic field, and the important characteristics of peristaltic flow, velocity, and pressure rise are calculated and analyzed; in the presence of an applied magnetic field, a vertical cylindrical natural convection magnetohydrodynamic flow is introduced by Kumar et al. [18]. Aleem et al. [19] explored the free convective and unsteady flow of Jeffrey fluid between the two hot vertical parallel plates in porous media and gave the accurate solutions of velocity and temperature with the Laplace transform method. Gao and Jian [20] analyzed the electroosmotic flow (EOF) of Jeffrey fluid in a circular microchannel under vertical magnetic field and got the accurate solutions of the velocity and volumetric flow rate based on the separation transformation method, and so on.The researches mentioned above are mostly about the flow and heat transfer characteristics for the Jeffrey fluid in microtubes under low Zeta potential (that is, the Zeta potential is less than 25 mV). However, in practical application, the wall Zeta potential of most interfaces is higher than 25 mV. Therefore, it is also of great significance to consider the high Zeta potential case (that is, the Zeta potential is higher than 25 mV) in microchannels to provide theoretical reference for performance prediction of electroosmotic pump. For example, the electrical potential distribution and velocity distribution were given for two-dimensional electroosmotic flow of power law model in circular and elliptical microchannels at high Zeta potential [21]; the electroosmotic flow for power law fluids through a curved rectangular microchannel was studied numerically at high Zeta potentials [22]; the rotating electroosmotic flow for Maxwell fluid at high Zeta potential was investigated in a parallel plate microchannel, and the velocity distribution was given by the finite difference method, and the effect of the wall Zeta potential and other related physical parameters on the velocity distribution was analyzed [23]; the electroosmotic flow for non-Newtonian fluid was explored at high Zeta potential through porous polymer membrane, and the variation of the overall flow rate with disparate physical parameters was evaluated [24], and so on.In summary, although the flow characteristics for some types of microfluids are given in micropipes at high Zeta potential, there are few results of the unsteady electroosmotic flow for Jeffrey fluid acted by the vertical magnetic field at high Zeta potential in micropipes. Therefore, the electroosmotic flow for Jeffrey fluid at high Zeta potential will be further studied in a circular microchannel under the impact of vertical magnetic field in this paper.


## 2. Mathematical Model and Formulation of the Problem

In this paper, the unsteady EOF for incompressible Jeffrey fluid acted by the vertical magnetic field is considered in a cylindrical microchannel at high Zeta potential. Mathematical model of the issues is presented in Figure1. The channel with length L has a circular cross section with radius R, and the wall surfaces are negatively charged. Setting up the cylindrical coordinate system r,θ,x as displayed in Figure 1 and assuming that the fluid is acted by the applied electric field E=Ex,Ey,0 and the uniform magnetic field in the vertical direction B=0,0,Bz, it flows along the axial direction under the coaction of the electric field force generated by the external electric field and the Lorentz force generated by the external magnetic field and pressure gradient, and the velocity is u=u¯,0,0. The density, viscosity, and conductivity of the fluid are ρ, μ, and σe, respectively.Figure 1
Sketch for EOF of Jeffrey fluid in a circular microchannel under vertical magnetic field.


### 2.1. Potential Distribution

The electrolyte solution interacts with the channel wall to form double electrical layer (EDL). In terms of the theory of electrostatics, the distribution of electric potentialψ¯ can be depicted by the Poisson equation as follows [7]:
(1)∇2ψ¯=−ρeε,where ε is the dielectric permittivity of the electrolyte solution.The ion concentration per unit volume of electrolyte solution obeying Boltzmann distribution is defined as [7]
(2)nv=nv0exp−zveψ¯kbT,where zv is the ion valence, nv0 is the ion concentration of the electrolyte solution, e is the amount of charge carried by the electron, T is the absolute temperature, and kb is the Boltzmann constant.The net electrical charge densityρe is shown by the following form:
(3)ρe=∑nvzve=e∑zvnv0exp−zveψ¯kbT.Forz+:z−=1:1 solution, Equation (3) can be reduced as
(4)ρe=−2zvn0esinh−zveψ¯kbT.For small values of electrical potential ψ¯ of the EDL (ψ≪25mV), the Debye-Hückel linearization approximation can be applied, which means physically that the electrical potential is small compared with the thermal energy of the charged species. The analytical solution of electrical potential can be obtained. However, most application requires or provides high Zeta potential (ψ¯≥25mV), and D-H linear approximation methods are no longer suitable.By substituting Equation (4) into Equation (1), the Poisson-Boltzmann equation for potential distribution ψ¯ can be written as
(5)∇2ψ¯=2zvn0eεsinhzveψ¯kbT.The boundary conditions are as follows(6)ψ¯r¯r¯=R=ψ¯0,dψ¯r¯dr¯r¯=0=0,where ψ¯0 is the wall Zeta potential.Introducing nondimensional variables as follows,(7)r=r¯R,ψ,ψ0=zvekbTψ¯,ψ¯0,K=κR,where κ=2n0zv2e2/εkbT1/2 is the Debye-Hückel parameter; 1/κ, called the Debye length, represents the thickness of EDL; and K is electrokinetic width of the EDL denoting the ratio of radius R of microchannel to Debye length 1/κ.Inserting Equation (7) into Equations (5) and (6), the dimensionless Poisson-Boltzmann equation and boundary conditions for the potential distribution have the following forms:
(8)∂2ψ∂r2+1r∂ψ∂r=K2sinhψ,(9)ψrr=1=ψ0,dψrdrr=0=0.


### 2.2. Velocity Distribution

The fluid velocity is dominated by following continuity equation and Cauchy momentum equation:(10)∇⋅u=0,(11)ρ∂u/∂t¯+u⋅∇⋅u=−∇p−∇⋅τ¯+f,where u is fluid velocity, p is pressure, t¯ is time, τ¯ is stress tensor, and f is volume force which is equal to the sum of electric field force ρexE, and the Lorentz force J×B is as follows:
(12)f=ρexE+J×B,where electric current density J satisfies Ohm’s law.
(13)J=σeE+u×B.The constitutive equation of Jeffrey fluid model is established on the basis of Maxwell fluid constitutive equation. The partial derivative of deformation velocity tensor with respect to time is introduced, and the time parameterλ¯2 is also introduced. The constitutive equation of Jeffrey fluid [13] can be expressed as
(14)1+λ¯1∂∂t¯τ¯=−η01+λ¯2∂∂t¯γ,where η0 is the zero shear rate viscosity, λ¯1 is the relaxation time, λ¯2 is the retardation time, and γ=∂ui/∂xj+∂uj/∂xi is the deformation rate tensor. Newton fluid and Maxwell fluid are typical special cases of Jeffrey fluid model. When λ¯1=0,λ¯2=0, Jeffrey fluid model is simplified as Newtonian fluid model. When λ¯1≠0,λ¯2=0, Jeffrey fluid is a typical Maxwell fluid model.By eliminating stress tensorτ from Equations (10)–(14), the governing equation and initial-boundary value conditions of the velocity distribution can be reduced as [20]
(15)1+λ¯1∂∂t¯ρ∂u¯∂t¯+dpdx−ρexEx+σeB2u¯−σeEyBz=η01+λ¯2∂∂t¯∂2u¯∂r¯2+1r¯∂u¯∂r¯,(16)u¯t¯=0=0,∂u¯∂t¯t¯=0=0,(17)u¯r¯=R=0,∂u¯∂r¯r¯=0=0.The dimensionless parameters such as the following are introduced to simplify the governing equations:(18)u=u¯vHS,t,λ1,λ2=t¯,λ¯1,λ¯2ρR2/η0,(19)Ha=BzRσeη0,β=EyRvHSσeη0,(20)vHS=−εkbTExψ0z0eη0,where Ha is Hartmann number which is the parameter to measure the ratio of magnetic and viscous forces, β is the nondimensional parameter representing the strength of the lateral direction electric field, and vHS represents the Helmholtz-Smoluchowski electroosmotic velocity.By substituting Equations (18)–(20) into Equations (15)–(17), the nondimensional equation and initial-boundary value conditions of the velocity distribution were obtained finally.
(21)1+λ1∂∂t∂u∂t−Ω−K2sinhψ+Ha2u−Haβ=1+λ2∂∂t∂2u∂r2+1r∂u∂r,(22)ut=0=0,∂u∂tt=0=0,s,(23)ur=1=0,∂u∂rr=0=0,where dp/dx=−Δp/L, Ω=vp/vHS is the ratio of axial pressure gradient velocity to electroosmosis driving velocity, and vp is the flow velocity driven by axial pressure, which can be expressed as vp=ΔpR2/η0L.


## 2.1. Potential Distribution

The electrolyte solution interacts with the channel wall to form double electrical layer (EDL). In terms of the theory of electrostatics, the distribution of electric potentialψ¯ can be depicted by the Poisson equation as follows [7]:
(1)∇2ψ¯=−ρeε,where ε is the dielectric permittivity of the electrolyte solution.The ion concentration per unit volume of electrolyte solution obeying Boltzmann distribution is defined as [7]
(2)nv=nv0exp−zveψ¯kbT,where zv is the ion valence, nv0 is the ion concentration of the electrolyte solution, e is the amount of charge carried by the electron, T is the absolute temperature, and kb is the Boltzmann constant.The net electrical charge densityρe is shown by the following form:
(3)ρe=∑nvzve=e∑zvnv0exp−zveψ¯kbT.Forz+:z−=1:1 solution, Equation (3) can be reduced as
(4)ρe=−2zvn0esinh−zveψ¯kbT.For small values of electrical potential ψ¯ of the EDL (ψ≪25mV), the Debye-Hückel linearization approximation can be applied, which means physically that the electrical potential is small compared with the thermal energy of the charged species. The analytical solution of electrical potential can be obtained. However, most application requires or provides high Zeta potential (ψ¯≥25mV), and D-H linear approximation methods are no longer suitable.By substituting Equation (4) into Equation (1), the Poisson-Boltzmann equation for potential distribution ψ¯ can be written as
(5)∇2ψ¯=2zvn0eεsinhzveψ¯kbT.The boundary conditions are as follows(6)ψ¯r¯r¯=R=ψ¯0,dψ¯r¯dr¯r¯=0=0,where ψ¯0 is the wall Zeta potential.Introducing nondimensional variables as follows,(7)r=r¯R,ψ,ψ0=zvekbTψ¯,ψ¯0,K=κR,where κ=2n0zv2e2/εkbT1/2 is the Debye-Hückel parameter; 1/κ, called the Debye length, represents the thickness of EDL; and K is electrokinetic width of the EDL denoting the ratio of radius R of microchannel to Debye length 1/κ.Inserting Equation (7) into Equations (5) and (6), the dimensionless Poisson-Boltzmann equation and boundary conditions for the potential distribution have the following forms:
(8)∂2ψ∂r2+1r∂ψ∂r=K2sinhψ,(9)ψrr=1=ψ0,dψrdrr=0=0.


## 2.2. Velocity Distribution

The fluid velocity is dominated by following continuity equation and Cauchy momentum equation:(10)∇⋅u=0,(11)ρ∂u/∂t¯+u⋅∇⋅u=−∇p−∇⋅τ¯+f,where u is fluid velocity, p is pressure, t¯ is time, τ¯ is stress tensor, and f is volume force which is equal to the sum of electric field force ρexE, and the Lorentz force J×B is as follows:
(12)f=ρexE+J×B,where electric current density J satisfies Ohm’s law.
(13)J=σeE+u×B.The constitutive equation of Jeffrey fluid model is established on the basis of Maxwell fluid constitutive equation. The partial derivative of deformation velocity tensor with respect to time is introduced, and the time parameterλ¯2 is also introduced. The constitutive equation of Jeffrey fluid [13] can be expressed as
(14)1+λ¯1∂∂t¯τ¯=−η01+λ¯2∂∂t¯γ,where η0 is the zero shear rate viscosity, λ¯1 is the relaxation time, λ¯2 is the retardation time, and γ=∂ui/∂xj+∂uj/∂xi is the deformation rate tensor. Newton fluid and Maxwell fluid are typical special cases of Jeffrey fluid model. When λ¯1=0,λ¯2=0, Jeffrey fluid model is simplified as Newtonian fluid model. When λ¯1≠0,λ¯2=0, Jeffrey fluid is a typical Maxwell fluid model.By eliminating stress tensorτ from Equations (10)–(14), the governing equation and initial-boundary value conditions of the velocity distribution can be reduced as [20]
(15)1+λ¯1∂∂t¯ρ∂u¯∂t¯+dpdx−ρexEx+σeB2u¯−σeEyBz=η01+λ¯2∂∂t¯∂2u¯∂r¯2+1r¯∂u¯∂r¯,(16)u¯t¯=0=0,∂u¯∂t¯t¯=0=0,(17)u¯r¯=R=0,∂u¯∂r¯r¯=0=0.The dimensionless parameters such as the following are introduced to simplify the governing equations:(18)u=u¯vHS,t,λ1,λ2=t¯,λ¯1,λ¯2ρR2/η0,(19)Ha=BzRσeη0,β=EyRvHSσeη0,(20)vHS=−εkbTExψ0z0eη0,where Ha is Hartmann number which is the parameter to measure the ratio of magnetic and viscous forces, β is the nondimensional parameter representing the strength of the lateral direction electric field, and vHS represents the Helmholtz-Smoluchowski electroosmotic velocity.By substituting Equations (18)–(20) into Equations (15)–(17), the nondimensional equation and initial-boundary value conditions of the velocity distribution were obtained finally.
(21)1+λ1∂∂t∂u∂t−Ω−K2sinhψ+Ha2u−Haβ=1+λ2∂∂t∂2u∂r2+1r∂u∂r,(22)ut=0=0,∂u∂tt=0=0,s,(23)ur=1=0,∂u∂rr=0=0,where dp/dx=−Δp/L, Ω=vp/vHS is the ratio of axial pressure gradient velocity to electroosmosis driving velocity, and vp is the flow velocity driven by axial pressure, which can be expressed as vp=ΔpR2/η0L.


## 3. Numerical Methods

In this paper, equations for the potential distribution and velocity distribution are solved by the Chebyshev spectral method and the finite difference method. The finite difference method is a local method; that is, the derivative of each position (∂u/∂x,∂2u/∂x2⋯⋯) is calculated from several adjacent points, and the matrices used in the calculation process are sparse matrices containing many zeros, so its accuracy is not high [25]. The Chebyshev spectral method is a global algorithm, which uses all known points to calculate the derivative of a certain point, greatly improving the accuracy [25]. Through the image comparison of the two methods mentioned above, we can see the advantages of the Chebyshev spectral method more clearly.


### 3.1. Numerical Methods for Potential Distribution

#### 3.1.1. The Chebyshev Spectral Method

We first choose the following Chebyshev points which are defined in the interval−1,1 [26]:
(24)xp=−cos2p+1π2N+1,p=0,⋯,N.We convert them into radial coordinates r, which are defined in the interval 0,1, by the relation.
(25)rp=1+xp2.Set thatψ=ψr0,ψr1,⋯,ψrN is the undetermined vector at the Chebyshev points, and we construct a Chebyshev polynomial Pr of order N or less, which satisfies the relationships Prp=ψrp,p=0,1,⋯,N. By taking the derivative of the Chebyshev polynomial P and evaluating it at the grid point, we can transform nonlinear Poisson-Boltzmann equation into linear algebraic equations and finally obtain the numerical solution of the dimensionless potential distribution from Equations (8) and (9) by Newton’s method.


#### 3.1.2. The Finite Difference Method

LetΔr=1/N, r=iΔr, and i=1,⋯,N−1; we obtain the discretization form of Equations (8) and (9) by the finite difference method:
(26)1Δr2ψi+1−2ψi+ψi−1+12iΔr2ψi+1−ψi−1−K2sinhψi=0,(27)ψ1=ψ0,ψN=ψ0,where ψi=ψri.The numerical solution of the dimensionless potential distributionψ is obtained from Equations (26) and (27) by Newton’s method.


#### 3.1.3. Comparison of Consequences

At low Zeta potential, the results obtained by using the Chebyshev spectral method and the finite difference method are compared with the approximate analytical solution obtained by the D-H linear approximation, respectively, as shown in Figure2(a). The results show that both numerical methods are feasible, and when the number of points (8 points in Figure 2(a)) is small, the accuracy of the finite difference method is poorer, while accuracy of the Chebyshev spectral method is higher. At high Zeta potential, the results by the Chebyshev spectral method only with 30 points are well coincident with the results by the finite difference method with 100 points, as shown in Figure 2(b). Therefore, the Chebyshev spectral method has higher accuracy and less computation than the finite difference method.Figure 2
Comparison of Chebyshev spectral method for nonlinear P-B equation with finite difference numerical solution and D-H linear approximate analytical solutionK=10.
(a)
ψ=1(b)
ψ=10


### 3.2. Numerical Methods for Velocity Distribution

#### 3.2.1. The Chebyshev Spectral Method

Introducing the functionvr,t=∂ur,t/∂t, reducing ∂2/∂t2 in Equation (21) to ∂/∂t, we get
(28)∂u∂t=v,∂v∂t=1λ1λ2∂2v∂r2+1r∂v∂r+∂2u∂r2+1r∂u∂r−Ha2u−1+λ1Ha2v+Ω+βHa+K2sinhψ,ut=0=0,vt=0=0,ur=1=0,∂u∂rr=0=0,vr=1=0,∂v∂rr=0=0.We use the derivative matrix of Chebyshev spectral to compute the derivative in ther direction and use fourth-fifth-order Runge-Kutta algorithm (ode45) to calculate ∂/∂t in Equation (28).


#### 3.2.2. The Finite Difference Method

LetΔr=1/N, Δt=T/M, r=jΔr, t=kΔt, j=1,⋯,N−1, and k=1,⋯,M−1, by using the finite difference scheme for Equations (21)–(23); we give the discretization forms of the dimensionless velocity distribution:
(29)λ24jΔtΔr−λ22ΔtΔr2uj−1k+1+λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk+1−λ24jΔtΔr+λ22ΔtΔr2uj+1k+1=1Δr2−12jΔr2uj−1k+2λ1Δt2−2Δr2ujk+1Δr2+12jΔr2uj+1k+λ24jΔtΔr−λ22ΔtΔr2uj−1k−1+−λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk−1−λ24jΔtΔr+λ22ΔtΔr2uj+1k−1+Ω+βHa+K2sinhψj,(30)uj0=uj1=0,(31)uNk=0,u0k=u1k,where ujk=urj,tk.The numerical solution of the dimensionless velocity distributionu is obtained from Equations (29)–(31) by Newton’s method.


#### 3.2.3. Comparison of Consequences

At low Zeta potential, the results acquired by using Chebyshev spectral method and the finite difference method are compared with the velocity distribution acquired by the D-H linear approximation (see Figure3(a)). By comparison, both numerical methods are viable, but we can see that taking a small number of points (6 points in Figure 3(a)), the Chebyshev spectral method has higher accuracy than the finite difference method. At high Zeta potential(Figure 3(b)), it also can be noted that when enough points (100 points in Figure 3(b)) are taken for the finite difference method, the velocity curve solved by the finite difference method can be relatively consistent with that solved by the Chebyshev spectral method only with 20 points in Figure 3(b), so compared with the finite difference method, the results obtained by the Chebyshev spectral method has higher accuracy and less computation. Therefore, it is feasible to extend Chebyshev spectral method to solve the fluid velocity distribution.Figure 3
Comparison of dimensionless velocity distribution solved by Chebyshev spectral method with those solved by D-H linear approximation and finite difference methodK=10,Ha=1,λ1=0.7,λ2=0.08,β=1,Ω=1,t=1.
(a)
ψ=1(b)
ψ=5


## 3.1. Numerical Methods for Potential Distribution

### 3.1.1. The Chebyshev Spectral Method

We first choose the following Chebyshev points which are defined in the interval−1,1 [26]:
(24)xp=−cos2p+1π2N+1,p=0,⋯,N.We convert them into radial coordinates r, which are defined in the interval 0,1, by the relation.
(25)rp=1+xp2.Set thatψ=ψr0,ψr1,⋯,ψrN is the undetermined vector at the Chebyshev points, and we construct a Chebyshev polynomial Pr of order N or less, which satisfies the relationships Prp=ψrp,p=0,1,⋯,N. By taking the derivative of the Chebyshev polynomial P and evaluating it at the grid point, we can transform nonlinear Poisson-Boltzmann equation into linear algebraic equations and finally obtain the numerical solution of the dimensionless potential distribution from Equations (8) and (9) by Newton’s method.


### 3.1.2. The Finite Difference Method

LetΔr=1/N, r=iΔr, and i=1,⋯,N−1; we obtain the discretization form of Equations (8) and (9) by the finite difference method:
(26)1Δr2ψi+1−2ψi+ψi−1+12iΔr2ψi+1−ψi−1−K2sinhψi=0,(27)ψ1=ψ0,ψN=ψ0,where ψi=ψri.The numerical solution of the dimensionless potential distributionψ is obtained from Equations (26) and (27) by Newton’s method.


### 3.1.3. Comparison of Consequences

At low Zeta potential, the results obtained by using the Chebyshev spectral method and the finite difference method are compared with the approximate analytical solution obtained by the D-H linear approximation, respectively, as shown in Figure2(a). The results show that both numerical methods are feasible, and when the number of points (8 points in Figure 2(a)) is small, the accuracy of the finite difference method is poorer, while accuracy of the Chebyshev spectral method is higher. At high Zeta potential, the results by the Chebyshev spectral method only with 30 points are well coincident with the results by the finite difference method with 100 points, as shown in Figure 2(b). Therefore, the Chebyshev spectral method has higher accuracy and less computation than the finite difference method.Figure 2
Comparison of Chebyshev spectral method for nonlinear P-B equation with finite difference numerical solution and D-H linear approximate analytical solutionK=10.
(a)
ψ=1(b)
ψ=10


## 3.1.1. The Chebyshev Spectral Method

We first choose the following Chebyshev points which are defined in the interval−1,1 [26]:
(24)xp=−cos2p+1π2N+1,p=0,⋯,N.We convert them into radial coordinates r, which are defined in the interval 0,1, by the relation.
(25)rp=1+xp2.Set thatψ=ψr0,ψr1,⋯,ψrN is the undetermined vector at the Chebyshev points, and we construct a Chebyshev polynomial Pr of order N or less, which satisfies the relationships Prp=ψrp,p=0,1,⋯,N. By taking the derivative of the Chebyshev polynomial P and evaluating it at the grid point, we can transform nonlinear Poisson-Boltzmann equation into linear algebraic equations and finally obtain the numerical solution of the dimensionless potential distribution from Equations (8) and (9) by Newton’s method.


## 3.1.2. The Finite Difference Method

LetΔr=1/N, r=iΔr, and i=1,⋯,N−1; we obtain the discretization form of Equations (8) and (9) by the finite difference method:
(26)1Δr2ψi+1−2ψi+ψi−1+12iΔr2ψi+1−ψi−1−K2sinhψi=0,(27)ψ1=ψ0,ψN=ψ0,where ψi=ψri.The numerical solution of the dimensionless potential distributionψ is obtained from Equations (26) and (27) by Newton’s method.


## 3.1.3. Comparison of Consequences

At low Zeta potential, the results obtained by using the Chebyshev spectral method and the finite difference method are compared with the approximate analytical solution obtained by the D-H linear approximation, respectively, as shown in Figure2(a). The results show that both numerical methods are feasible, and when the number of points (8 points in Figure 2(a)) is small, the accuracy of the finite difference method is poorer, while accuracy of the Chebyshev spectral method is higher. At high Zeta potential, the results by the Chebyshev spectral method only with 30 points are well coincident with the results by the finite difference method with 100 points, as shown in Figure 2(b). Therefore, the Chebyshev spectral method has higher accuracy and less computation than the finite difference method.Figure 2
Comparison of Chebyshev spectral method for nonlinear P-B equation with finite difference numerical solution and D-H linear approximate analytical solutionK=10.
(a)
ψ=1(b)
ψ=10


## 3.2. Numerical Methods for Velocity Distribution

### 3.2.1. The Chebyshev Spectral Method

Introducing the functionvr,t=∂ur,t/∂t, reducing ∂2/∂t2 in Equation (21) to ∂/∂t, we get
(28)∂u∂t=v,∂v∂t=1λ1λ2∂2v∂r2+1r∂v∂r+∂2u∂r2+1r∂u∂r−Ha2u−1+λ1Ha2v+Ω+βHa+K2sinhψ,ut=0=0,vt=0=0,ur=1=0,∂u∂rr=0=0,vr=1=0,∂v∂rr=0=0.We use the derivative matrix of Chebyshev spectral to compute the derivative in ther direction and use fourth-fifth-order Runge-Kutta algorithm (ode45) to calculate ∂/∂t in Equation (28).


### 3.2.2. The Finite Difference Method

LetΔr=1/N, Δt=T/M, r=jΔr, t=kΔt, j=1,⋯,N−1, and k=1,⋯,M−1, by using the finite difference scheme for Equations (21)–(23); we give the discretization forms of the dimensionless velocity distribution:
(29)λ24jΔtΔr−λ22ΔtΔr2uj−1k+1+λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk+1−λ24jΔtΔr+λ22ΔtΔr2uj+1k+1=1Δr2−12jΔr2uj−1k+2λ1Δt2−2Δr2ujk+1Δr2+12jΔr2uj+1k+λ24jΔtΔr−λ22ΔtΔr2uj−1k−1+−λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk−1−λ24jΔtΔr+λ22ΔtΔr2uj+1k−1+Ω+βHa+K2sinhψj,(30)uj0=uj1=0,(31)uNk=0,u0k=u1k,where ujk=urj,tk.The numerical solution of the dimensionless velocity distributionu is obtained from Equations (29)–(31) by Newton’s method.


### 3.2.3. Comparison of Consequences

At low Zeta potential, the results acquired by using Chebyshev spectral method and the finite difference method are compared with the velocity distribution acquired by the D-H linear approximation (see Figure3(a)). By comparison, both numerical methods are viable, but we can see that taking a small number of points (6 points in Figure 3(a)), the Chebyshev spectral method has higher accuracy than the finite difference method. At high Zeta potential(Figure 3(b)), it also can be noted that when enough points (100 points in Figure 3(b)) are taken for the finite difference method, the velocity curve solved by the finite difference method can be relatively consistent with that solved by the Chebyshev spectral method only with 20 points in Figure 3(b), so compared with the finite difference method, the results obtained by the Chebyshev spectral method has higher accuracy and less computation. Therefore, it is feasible to extend Chebyshev spectral method to solve the fluid velocity distribution.Figure 3
Comparison of dimensionless velocity distribution solved by Chebyshev spectral method with those solved by D-H linear approximation and finite difference methodK=10,Ha=1,λ1=0.7,λ2=0.08,β=1,Ω=1,t=1.
(a)
ψ=1(b)
ψ=5


## 3.2.1. The Chebyshev Spectral Method

Introducing the functionvr,t=∂ur,t/∂t, reducing ∂2/∂t2 in Equation (21) to ∂/∂t, we get
(28)∂u∂t=v,∂v∂t=1λ1λ2∂2v∂r2+1r∂v∂r+∂2u∂r2+1r∂u∂r−Ha2u−1+λ1Ha2v+Ω+βHa+K2sinhψ,ut=0=0,vt=0=0,ur=1=0,∂u∂rr=0=0,vr=1=0,∂v∂rr=0=0.We use the derivative matrix of Chebyshev spectral to compute the derivative in ther direction and use fourth-fifth-order Runge-Kutta algorithm (ode45) to calculate ∂/∂t in Equation (28).


## 3.2.2. The Finite Difference Method

LetΔr=1/N, Δt=T/M, r=jΔr, t=kΔt, j=1,⋯,N−1, and k=1,⋯,M−1, by using the finite difference scheme for Equations (21)–(23); we give the discretization forms of the dimensionless velocity distribution:
(29)λ24jΔtΔr−λ22ΔtΔr2uj−1k+1+λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk+1−λ24jΔtΔr+λ22ΔtΔr2uj+1k+1=1Δr2−12jΔr2uj−1k+2λ1Δt2−2Δr2ujk+1Δr2+12jΔr2uj+1k+λ24jΔtΔr−λ22ΔtΔr2uj−1k−1+−λ1Δt2+1+λ1Ha22Δt+λ2ΔtΔr2ujk−1−λ24jΔtΔr+λ22ΔtΔr2uj+1k−1+Ω+βHa+K2sinhψj,(30)uj0=uj1=0,(31)uNk=0,u0k=u1k,where ujk=urj,tk.The numerical solution of the dimensionless velocity distributionu is obtained from Equations (29)–(31) by Newton’s method.


## 3.2.3. Comparison of Consequences

At low Zeta potential, the results acquired by using Chebyshev spectral method and the finite difference method are compared with the velocity distribution acquired by the D-H linear approximation (see Figure3(a)). By comparison, both numerical methods are viable, but we can see that taking a small number of points (6 points in Figure 3(a)), the Chebyshev spectral method has higher accuracy than the finite difference method. At high Zeta potential(Figure 3(b)), it also can be noted that when enough points (100 points in Figure 3(b)) are taken for the finite difference method, the velocity curve solved by the finite difference method can be relatively consistent with that solved by the Chebyshev spectral method only with 20 points in Figure 3(b), so compared with the finite difference method, the results obtained by the Chebyshev spectral method has higher accuracy and less computation. Therefore, it is feasible to extend Chebyshev spectral method to solve the fluid velocity distribution.Figure 3
Comparison of dimensionless velocity distribution solved by Chebyshev spectral method with those solved by D-H linear approximation and finite difference methodK=10,Ha=1,λ1=0.7,λ2=0.08,β=1,Ω=1,t=1.
(a)
ψ=1(b)
ψ=5


## 4. Results and Discussion

The main intention of this study is to investigate the unsteady electroosmotic flow for the Jeffrey fluid at high Zeta potential under the action of vertical magnetic field in a circular microchannel. Consequently, the impression of different parameters on the velocity distributionu will be discussed in this section. Computations have been carried out for a given set of parameter values [20]: the characteristic scale of the circular micropump is R=100μm, the density of the fluid is ρ=1.06×103kg×m−3, the viscous shear rate is η0=3×10−3kg×m−1s−1, the conductivity is σe=2.2×10−4−106S/m [27], and the lower value 2.2×10−4S/m of the conductivity is the deionized water, while the upper value 106S/m of the conductivity corresponds to the liquid metal such as mercury. If the vertical magnetic field B=0.018−0.44T [28], the range of Ha can be obtained as 0.0001≤Ha≤5. If the intensity of the alternating electric field is 0≤Ey≤1V/m and vHS=100μm/s, then 0≤β≤6×104. The parameter value of the relaxation time λ1 is wide [13, 29] and can be changed from 10−4s≤λ1≤103s. In addition, because the retardation time λ2 is normally smaller than the relaxation time λ1, the retardation time λ2 and the relaxation time λ1 satisfy λ2<λ1.Figure4 depicts the variation of velocity u at different relaxation times λ1 and retardation time λ2. It is visible from Figure 4 that the amplitude of velocity enhances with increase of relaxation time λ1. This is because the greater the relaxation time is, the greater the elastic effect of the Jeffery fluid is, resulting in the increase in velocity (see Figures 4(a)–4(c)). It is also noted that the increase in the retardation time λ2 leads to decrease of fluid velocity u, and the longer the retardation time, the more obvious the retardation effect. That is because longer retardation time will cause greater flow resistance (see Figures 4(d)–4(f)). Moreover, the velocity tends to be stable with the increase of time. Figures 4(g) and 4(h) illustrate the effects of different relaxation time λ1 and retardation time λ2 on the velocity distribution u at the center of the microchannel. Here, Figure 4(g) is the velocity diagram at r=0 in Figures 4(a)–4(c), and Figure 4(h) is the velocity diagram at r=0 in Figures 4(d)–4(f). It can be seen obviously from Figures 4(g) and 4(h) that the velocity increases with increasing the relaxation time λ1, whereas it decreases with increasing the retardation time λ2.Figure 4
A dimensionless velocity distribution in a circular microchannel (λ1=0.7,λ2=0.08,β=1,Ω=1,K=10,ψ0=5,Ha=1).
(a)
λ1=0.1(b)
λ1=0.3(c)
λ1=0.5(d)
λ2=0.2(e)
λ2=0.4(f)
λ2=0.6(g)(h)Figure5 delineates the changing trend of velocity with Hartmann number at the center of microchannel for relaxation time λ1 and retardation time λ2. As shown in Figure 5, the velocity u first increases with increase of Hartmann number Ha and then decreases with increase of Hartmann number Ha, and there is a threshold of Hartmann number Ha that maximizes the velocity. It is mainly because the Lorentz force in Equation (15) consists of two parts: one is the driving force σeEyB and the other is the hindering force σeB2u. For the smaller value Ha, the driving force of Lorentz force is greater than its hindrance force, so the velocity increase with the enhance of Ha, and with the increase of Ha, the hindering force of Lorentz force exceeds its driving force, so the velocity diminishes with the enhance of Ha. In addition, when the Ha is comparatively small, the velocity increases with the increase of the relaxation time λ1, but when Ha gets bigger, the velocity decays with the increase of relaxation time λ1 (see Figure 5(a)). No matter what the value Ha is, the velocity always goes down for the increase of retardation time λ2 (see Figure 5(b)).Figure 5
Dimensionless velocity distribution atr=0 in a circular microchannel (λ1=0.7b,λ2=0.08a,β=1,Ω=1,K=10,ψ0=5,t=1).
(a)(b)In Figure6, the changing trend of the velocity for different electric width K, pressure gradient Ω, and wall Zeta potential ψ0 in the microchannel is exhibited. Figure 6 shows that the positive pressure gradient (Ω>0) promotes the flow of the fluid, whereas the negative pressure gradient (Ω<0) impedes the flow of the fluid for constant electric width K and wall Zeta potential ψ0. They also display that the velocity increases with the increase of electric width K when wall Zeta potential ψ0 and pressure gradient Ω are constant. Furthermore, high Zeta potential ψ0 results in the large magnitude of velocity for constant electric width K and pressure gradient Ω, especially within the EDL. The main reason is that there is no sufficient time for the flow diffusing far into the midplane of the wall, and the velocity variations are restricted only within the EDL.Figure 6
Dimensionless velocity distribution in a circular microchannel (λ1=0.7,λ2=0.08,β=1,t=1,K=10,Ha=1).
(a)
Ω=−1(b)
Ω=0(c)
Ω=1(d)
ψ0=5The variations of the dimensionless velocity distribution with the nondimensional parameterβ and wall Zeta potential ψ0 is drawn in Figure 7. It is noticed that the velocity u increases with the increase of β. This is due to the fact that the increase in β leads to the enhancement of the Lorentz force, and it in turn increases the velocity. It is also noted in Figure 7 that the maximum variation area of the velocity distribution is limited to the vicinity of the wall surface of the microchannel, that is, the EDL. In the EDL, the velocity increases rapidly from zero to the maximum because of the interaction between the external driving force and the EDL. Outside the EDL, the change of the velocity is small.Figure 7
Dimensionless velocity distribution in a circular microchannel (λ1=0.7,λ2=0.08,Ω=1,K=10,t=0.5,Ha=1).


## 5. Conclusion

In this article, the unsteady EOF flow of Jeffrey fluid coacted by the vertical magnetic field, external electric field, and pressure is investigated in circular microchannel at high Zeta potential. Through the use of the Chebyshev spectral method and finite difference method, the numerical solutions of the electric potential distribution and the velocity distribution are solved. By comparing the above two methods, we find that the Chebyshev spectral method is more accurate and has less computation than finite difference method. Detailed study reveals that the velocity distribution hinges on Hartmann numberHa, wall Zeta potential ψ0, relaxation time λ1, retardation time λ2, the ratio of axial pressure gradient to electroosmosis driving velocity Ω, strength of the lateral direction electric field β, and electric width K. The main conclusions are as follows:
(1)
Finite difference method and Chebyshev spectral method are effective methods to solve the potential and velocity distribution of microfluid under high Zeta potential, and the Chebyshev spectral method has higher accuracy and less calculation(2)
Due to the elastic action of Jeffery fluid, the greater the relaxation timeλ1 is, the greater the velocity profile is; the longer the retardation time λ2, the greater the flow resistance and the smaller the velocity, and the speed u tends to be stable gradually with the increase of time t(3)
WhenHa is smaller, an enhancement in the relaxation time λ1 results in the addition of the velocity u, but when Ha is larger, an increase in the relaxation time λ1 results in the decline in the velocity u. Whatever the value Ha is, the velocity always goes down for the increase of retardation time λ2(4)
The positive pressure gradient (Ω>0) promotes the flow of the fluid, whereas the negative pressure gradient (Ω<0) impedes the flow of the fluid(5)
The velocity increases with the increase of wall Zeta potentialψ0 and electric width K, and the maximum variation area of velocity u is restricted within the EDL(6)
The velocityu increases with the increase of the strength of the lateral direction electric field βBased on the research in this paper, the following issues can be further considered in the future:(1)
The flow and heat transfer characteristics of the two-layer fluid in different micropipes under high Zeta potential can be further studied(2)
It can be further expanded to study the flow and heat transfer characteristics of nanofluids in different micropipes such as T-shaped micropipes and semicylindrical micropipes at high Zeta potential under the action of external vertical magnetic field



---

*Source: 6247746-2022-07-18.xml*