- Email: [email protected]
Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids
Journal Pre-proof Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids Xiaoping Wang, Huanying Xu, Haitao Qi
PII: DOI: Reference:
S0893-9659(19)30505-1 https://doi.org/10.1016/j.aml.2019.106179 AML 106179
To appear in:
Applied Mathematics Letters
Received date : 16 October 2019 Revised date : 3 December 2019 Accepted date : 3 December 2019 Please cite this article as: X. Wang, H. Xu and H. Qi, Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids, Applied Mathematics Letters (2019), doi: https://doi.org/10.1016/j.aml.2019.106179. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.
*Manuscript Click here to view linked References
Journal Pre-proof
Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids Xiaoping Wang, Huanying Xu, Haitao Qi∗ School of Mathematics and Statistics, Shandong University, Weihai 264209, PR China
Abstract
na lP repr oo f
In this paper, a rotating electro-osmotic flow of a fractional Maxwell fluid in a parallel plate microchannel with high zeta potentials is examined. The Navier’s slip law at walls is considered. The electric double layer potential distribution is derived by using the nonlinear Poisson-Boltzmann equation. Based on the L1 approximation of the Caputo derivative, a Crank-Nicolson numerical scheme is developed for obtaining the numerical solutions of the rotating electro-osmotic flow velocity profiles. With a purpose to verify the correctness of our numerical results, a comparison has been made with the analytical solutions of the Newtonian fluid given by the previous work and the excellent agreement between the solutions is clear. Finally, the influences of the fractional parameters of α and β, the slip length d and the wall zeta potential ζ on the velocity distribution are also discussed in detail. Keywords: Rotating electro-osmotic flow, Fractional calculus, Finite difference method, Slip boundary condition, High zeta potential
1. Introduction
Jo
ur
In recent years, electro-osmotic flow (EOF) has become a more noticeable topic in microfluidics due to its operational advantages, such as better flow control, high reliability and low noise and has found applications in a wide range of fields, including chemical, biological analyses and medical diagnoses [1]. On the study of EOF, many researchers have concentrated on the non-Newtonian fluid. Among of those, Sarma et al. [2] reported the EOF of the simplified Phan-Thien-Tanner model in a parallel plate microchannel with high zeta potential. Peralta et al. [3] theoretically analyzed the pulsatile EOF of a Maxwell fluid in a parallel flat plate microchannel with asymmetric wall zeta potentials. Recently, the constitutive equations with fractional derivative have been widely used to describe the rheological behaviour of non-Newtonian fluid in the EOF. Wang et al. [4] derived the analytical and numerical solutions of the electro-osmotic slip flow of the fractional second grade fluids. Wang and Zhao [5] and Yang et al. [6] analytically and numerically analyzed the EOF of fractional Maxwell fluids. Zhao and Yang [7] investigated the transient EOF of fractional Oldroyd-B fluids in a rectangular microchannel. On the other hand, microfluidic systems may be placed in a rotating environment in many practical engineering applications, such as in centrifuges for mass separation and flow control. The rotating EOF is considered as an effective technique for mixing different liquids. Thus, the effect of rotating on EOF draws the attention of many researchers. Analytical and numerical works on the rotating EOF of Newtonian and non-Newtonian fluids in a parallel plate microchannel were extensively studied, see Refs. [8–10] and references therein. Based on the nonlinear Poisson-Boltzmann equation, Xie and Jian [11] numerically studied the rotating EOF of power-law fluids in a slit microchannel at high zeta potentials. More recently, Xu and Jian [12] firstly introduced the time fractional CaputoFabrizio derivative to study the rotating EOF. But we found that most of those researches are confined to the EOF in a uniform cross-section channels with constant zeta potential and no slip velocity conditions. However, some experiments and molecular dynamic simulations have demonstrated that the fluid can slip at solid surfaces. And at the level of micro-scales, the velocity slip at the liquid-solid interface plays an important role [13]. By considering Navier’s ∗ Corresponding
author. Tel./Fax: +86 631 5687309. Email address: [email protected] (Haitao Qi)
Preprint submitted to Applied Mathematics Letters
December 3, 2019
Journal Pre-proof
slip condition, Shit et al. [14] and Wang et al. [15] have studied the effects of slip velocity on the rotating EOF in a non-uniform microchannel. But to our best knowledge, the published researches on the rotating EOF of the non-Newtonian fluids is still rare and no attempt has been made regarding the rotating EOF of the viscoelastic fluids with fractional derivative. In this paper, we study the rotating EOF with the fractional Maxwell constitutive relationship in a parallel plate microchannel at high zeta potentials [16] 1 + λα−β Dα−β σ(t) = Eλα Dαt γ(t), 0 ≤ β ≤ α ≤ 1, (1) t where σ(t) is the shear stress, γ(t) is the shear strain, λ = η/E is the relaxation time, in which η is the viscosity constant and E is the shear modulus. Further, α and β are the Caputo fractional derivative parameters. 2. Model description and mathematical formulation
na lP repr oo f
We analyze the rotating EOF of the fractional Maxwell fluids in a parallel plate microchannel. An illustration of the problem and the coordinate system is similar to Figure 1 in the Refs [9, 11]. The origin of the coordinate system is placed in the center of the microchannel. The x axis along the horizontal direction. The y axis is pointing toward the inside of the plane and the z axis is vertical to the parallel plate separating by a distance of 2H. We assume that the microchannel is filled with a liquid solution of dielectric constant ε, and the plates are uniformly charged with a zeta potential ζ that develops an interfacial region of net charge density close to the channel walls, called as the EDL. An external electric field E = (E x , 0, 0) is imposed along the x-axis direction. We consider that the fluid is initially at rest, then the entire system is rotating about the z axis with a constant angular velocity Ω. In a horizontal rotating fluid layer, all the velocity is relative velocity and is considered in non-inertial frame of reference. Because of symmetry, only the upper half domain of the slit microchannel (0 ≤ z ≤ H) is considered. The governing equations for this type incompressible fluid flow are the continuity equation ∇ · V = 0,
(2)
! ∂V + V · ∇V + 2Ω × V = −∇P + ∇ · σ + F, ρ ∂t
(3)
and the general Cauchy momentum equation
where ρ is the fluid density, V = (u, v, 0) is the velocity vector of the rotating EOF flow, Ω = (0, 0, Ω) is the angular velocity vector about z axis and P = p − ρ | Ω × r2 | /2 is the modified fluid pressure due to centrifugal force where r = (x, y, z) and p is pressure [11, 14]. Further, σ is the stress tensor and F = (ρe E x f (t), 0, 0) is the body force vector where f (t) = H(t) is the well-known Heaviside unit step function [6, 17]. From the theory of electrostatics, the net electric charge density ρe in the diffuse layer of the wall EDL and the EDL potential ψ satisfy the following Poisson-Boltzmann equation ! d2 ψ ρe z0 eψ = − , ρe = −2n0 z0 e sinh . (4) ε κb T dz2
ur
Here n0 is the bulk concentration of the ions, z0 is the valence of the ions, e is the elementary electronic charge, ε is the dielectric constant of the fluid, kB is the Boltzmann constant and T is the absolute temperature. The boundary conditions for the potential given in the following form dψ = 0, ψ(H) = ζ. (5) dz z=0
Jo
Then the analytical expression of the potential distribution in the dimensionless form can be obtained as [2] ζ ψ = 4 tanh−1 tanh eK(z−1) , 4
(6)
by introducing the following non-dimensional variables
ψ ζ z ψ¯ = , ζ¯ = , z¯ = , K = κH. ψz ψz H 2
(7)
Journal Pre-proof
Here ψz = kb T/z0 e, κ = (2n0 z20 e2 /(εkb T ))1/2 is the reciprocal of the EDL thickness (λ = κ−1 ). In view of the above-mentioned assumptions and in the absence of the modified pressure gradient, the continuity Eq. (2) is satisfied automatically. For this specific problem, the constitutive equation for a fractional Maxwell fluid (1) can be expressed as ! ! α−β α−β α α−1 ∂v α α−1 ∂u , 1 + λ D σ = Eλ D , (8) 1 + λα−β Dα−β σ = Eλ D zy zx t t t t ∂z ∂z and the general Cauchy momentum equation (3) is simplified to ! ! ∂σzy ∂σzx ∂v ∂u − 2Ωv = + ρe E x f (t), ρ + 2Ωu = , ρ ∂t ∂z ∂t ∂z
(9)
which is subject to the following initial and boundary conditions u = 0,
v = 0,
∂u ∂v = 0, = 0, at t = 0, ∂t ∂t
(10)
na lP repr oo f
∂v ∂u = 0, = 0, at z = 0, (11) ∂z ∂z ∂v ∂u = 0, v + d = 0, at z = H. (12) u+d ∂z ∂z where d is the slip length. Considering that the bottom and upper walls of the microchannel are mostly made of the same material leading to the same hydrophobic characteristics and the same slip boundary condition. So we use the same slip length. Substituting Eq. (8) into Eq. (9) and introducing the following dimensionless variables u¯ =
u v t λ z d , v¯ = , t¯ = , λ¯ = , z¯ = , d¯ = , US US 1/Ω 1/Ω H H
(13)
where US = −εζE x /η is the Helmholtz-Smoluchowski velocity and η = λE, we obtain the dimensionless governing equations as: ! ! 2 ∂u K2 α−1 α−1 ∂ u (14) − 2v = λ D + sinh ψ 1 + λα−β Dtα−β f (t), Re 1 + λα−β Dα−β t t 2 ∂t ζ ∂z ! ! 2 ∂v α−β α−β α−1 α−1 ∂ v Re 1 + λ Dt . (15) + 2u = λ Dt ∂t ∂z2 Here, Re = ρΩH 2 /η is the rotating Reynolds number, f (t) = H(t) is the well-known Heaviside unit step function, which can be approximated by 1/(1 + exp((−t + p)/q)) with p → 0 and q → 0 in the following analysis. By making some transformations, the equivalent form of Eqs. (14) and (15) can be written as
ur
2 K2 1−α α−1 ∂ u 1−α α−β 1−β Re D2−α u + λα−β D2−β v − 2λα−β D1−β + sinh ψ D + λ D f (t), t t t u − 2Dt t v =λ t ζ ∂z2
Jo
2 1−α α−1 ∂ v Re D2−α v + λα−β D2−β u + 2λα−β D1−β , t t v + 2Dt t u =λ ∂z2 which has the following dimensionless initial and boundary conditions:
u = 0,
v = 0,
∂u ∂v = 0, = 0, at t = 0, ∂t ∂t
∂u ∂v = 0, = 0, at z = 0, ∂z ∂z ∂u ∂v u+d = 0, v + d = 0, at z = 1. ∂z ∂z 3
(16) (17)
(18) (19) (20)
Journal Pre-proof
3. Numerical algorithms In this section, we are using a finite difference method to solve the Eqs. (16)-(20). We first introduce a uniform discretization over the rectangular computational domain (D = {(z, t) | z, t ∈ [0, 1] × [0, Tˆ ]}) given by zi = ih, i = 0, 1, 2, . . . M; tn = nτ, n = 0, 1, 2, . . . N, in which h = 1/M and τ = Tˆ /N are the spatial and temporal step sizes, respectively. We define uni and vni as the numerical solutions of u(z, t) and v(z, t) at the mesh point (zi , tn ) Base on the L1 approximations of the Caputo fractional derivative, the Crank-Nicolson numerical schemes of Eqs. (16) and (17) can be obtained as n−1 n−1 k−1/2 Reλα−β τβ−1 (β) n−1/2 X k−1/2 Reτα−1 (α) n−1/2 X (α) (β) (β) (α) + b δt ui b δt ui − bn−k−1 − bn−k δt ui − bn−k−1 − bn−k δt ui Γ(1 + α) 0 Γ(1 + β) 0 k=1 k=1 n−1 n−2 X Reτα−1 (α) Reτα−1 (α) n X (α) (α) (α) (α) k k n−1 b v − b v − − b − bn−k vi − bn−k−2 − bn−k−1 vi Γ(1 + α) 0 i k=1 n−k−1 Γ(1 + α) 0 i k=1 n−1 n−2 α−β β−1 α−β β−1 X Reλ τ (β) n−1 X Reλ τ (β) n (β) (β) (β) (β) k k b v − b v − − b − bn−k vi − bn−k−2 − bn−k−1 vi Γ(1 + β) 0 i k=1 n−k−1 Γ(1 + β) 0 i k=1
na lP repr oo f
n−1 n−1 uni+1 − 2uni + uni−1 + un−1 + ui−1 K 2 sinh ψ 1−α n λα−β K 2 sinh ψ 1−β n i+1 − 2ui n−1 + D Dt ( f + f n−1 ), ( f + f ) + t 2ζ 2ζ 2h2 i = 1, 2, 3, . . . , M − 1, n = 1, 2, 3, . . . , N, (21)
= λα−1
n−1 n−1 k−1/2 k−1/2 Reλα−β τβ−1 (β) n−1/2 X Reτα−1 (α) n−1/2 X (α) (β) (β) (α) + b δt v b δt vi − bn−k−1 − bn−k δt vi − bn−k−1 − bn−k δt vi i Γ(1 + α) 0 Γ(1 + β) 0 k=1 k=1 n−1 n−2 X Reτα−1 (α) Reτα−1 (α) n X (α) n−1 k b0 ui − + b u − bn−k−1 − b(α) u + b(α) − b(α) uki 0 i i n−k n−k−2 n−k−1 Γ(1 + α) Γ(1 + α) k=1 k=1 n−1 n−2 α−β β−1 α−β β−1 X Reλ τ (β) n−1 X Reλ τ (β) n (β) (β) (β) (β) k k b u − b u − bn−k−1 − bn−k ui + bn−k−2 − bn−k−1 ui + Γ(1 + β) 0 i Γ(1 + β) 0 i k=1
n−1 vn − 2vni + vni−1 + vn−1 + vn−1 i+1 − 2vi i−1 = λα−1 i+1 , 2h2
i = 1, 2, 3, . . . , M − 1,
k=1
n = 1, 2, 3, . . . , N,
(22)
where bα0 = bβ0 = 1, bαj = ( j + 1)α − jα , bβj = ( j + 1)β − jβ , j = 0, 1, 2, .... The discretized initial and boundary conditions are u0i = 0, v0i = 0,
u1i − u0i = 0, τ
v1i − v0i = 0, i = 0, 1, 2, . . . , M, τ
Jo
ur
un1 − un0 vn1 − vn0 = 0, = 0, n = 0, 1, 2, . . . , N, h h un − unM−1 vn − vnM−1 unM + d M = 0, vnM + d M = 0, n = 0, 1, 2, . . . , N. h h
(23) (24) (25)
4. Results and discussion
In the present work, all the parameters used in the governing equations are dimensionless. In order to verify the accuracy of the numerical method and the theoretical analysis in our study. We compare the obtained numerical solutions of the velocity distribution by taking α = 1, β = 0, λ = 0, d = 0, K = 10, ζ = 1, Re = 1 and t = 5 with the analytical solutions for the rotating EOF of the Newtonian fluid with the no slip boundary condition presented by Chang and Wang [8]. Results are shown in Fig. 1. It can be seen that the numerical results are agreed well with the 4
Journal Pre-proof
analytical solutions. This demonstrates that the numerical method provided in Sect. 3 is effective and feasible and can be used to compute the rotating EOF velocity profiles of the fractional Maxwell fluid. To analyze the influence of the fractional parameters α and β, the slip length d and the wall zeta potential ζ on the velocity distribution. The values of the dimensionless parameters for the EOF in a microchannel are chosen as follows: the electro-kinetic width K=10 and the Reynolds number Re = 1 [8, 14], the relaxation time λ = 0.5 [3]. Fig. 2 illustrates the variations of the rotating EOF velocity for different values of the fractional parameter α (α=0.7, 1) and the slip length d (d = 0 [8, 9, 11]: for non-slip flow and d = 0.003, 0.008 [13]: for slip flow) at t = 5. As shown in the Fig. 2(a), the axis velocity magnitude increases with α and the presence of slip velocity has an effect on the rotating EOF velocity. We can also see that the velocity profile is steep in a region near to the microchannel wall while maintains a constant value near the central axis of the microchannel. Fig. 3 presents the variations of the rotating EOF velocity for different values of the fractional parameter β (β=0.3, 0.6) and the wall zeta potential ζ (ζ = 1 [3, 8, 9]: for low zeta potential and ζ = 3, 5 [2]: for high zeta potential) at t = 5. It can be seen from the Fig. 3(a) that increasing the values of β tends to sharply enhance the velocity distribution in x direction near the wall and the peaks of velocity profiles under the high zeta potential are higher than that the case of low zeta potential. However, the velocity in y direction becomes more and more negative as the increasing α, β, d and ζ. 0
0.7
na lP repr oo f
0.8
Present numerical solution Analytical solution of Chang and Wang
-0.1
0.6
-0.2
v
u
0.5 0.4 0.3
-0.3 -0.4
0.2
Present numerical solution Analytical solution of Chang and Wang
0.1
-0.5 -0.6
0 0
0.2
0.4
z
0.6
0.8
0
1
(a) Dimensionless velocity u in x direction
0.2
0.4
0.6
0.8
1
z
(b) Dimensionless velocity v in y direction
Figure 1: Comparison of our present numerical solution with the previous analytical results presented by Chang and Wang [8].
0
1
-0.1
0.8
-0.2
0.6
α=0.7
-0.3
v
u
α=1
-0.4
0.4
-0.5
α=0.7
0.2
-0.6
α=1
-0.7
0 0
0.2
0.4
0.6
0.8
0
1
ur
z
(a) Dimensionless velocity u in x direction
0.2
0.4
0.6
0.8
1
z
(b) Dimensionless velocity v in y direction
Jo
Figure 2: Variation of the velocity distribution for different values of α and d with β = 0.3, ζ = 3 at t = 5 (d = 0 solid lines, d = 0.003 dashed lines, d = 0.008 dash-dot lines).
5. Conclusions
In this paper, we study the rotating EOF of fractional Maxwell fluid through a parallel plate microchannel. The finite difference method is used to solve the governing equations of the velocity distribution. The effectiveness of our method is proved by giving a comparison with the known analytical results of the Newtonian fluid. The effects of 5
Journal Pre-proof
0.7
0 β=0.6
0.6
β=0.3
-0.1 0.5
-0.2
u
v
0.4 0.3
-0.3
β=0.3 0.2
-0.4
β=0.6
0.1
-0.5
0 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
0.6
0.8
1
z
z
(a) Dimensionless velocity u in x direction
(b) Dimensionless velocity v in y direction
Figure 3: Variation of the velocity distribution for different values of β and ζ with α = 0.8, d = 0.008 at t = 5 (ζ = 1 solid lines, ζ = 3 dashed lines, ζ = 5 dash-dot lines).
Acknowledgements
na lP repr oo f
the fractional parameters α and β, the slip length d and the wall zeta potential ζ on the velocity profiles are analyzed graphically. Results indicate that the axis velocity profiles increase with the fractional parameters α and β. The presence of slip velocity can enhance the magnitude of the axis velocity distribution and the peaks of the axis velocity profiles under the high zeta potential are higher than that the case of low zeta potential. The velocity in y direction becomes more and more negative as the increasing α, β, d and ζ.
This work was supported by the National Natural Science Foundation of China (NSFC No. 11672163, No. 11771254), the Fundamental Research Funds for the Central Universities (No. 2019ZRJC002). References
Jo
ur
[1] H. Bruus, Theoretical Microfluidics, Oxford University Press, New York, 2008. [2] R. Sarma, N. Deka, K. Sarma, P.K. Mondal, Electroosmotic flow of Phan-Thien-Tanner fluids at high zeta potentials: An exact analytical solution, Phys. Fluids 30 (2018) 062001. [3] M. Peralta, O. Bautista, F. M´endez, E. Bautista, Pulsatile electroosmotic flow of a Maxwell fluid in a parallel flat plate microchannel with asymmetric zeta potentials, Appl. Math. Mech. -Engl. Ed. 39 (2018) 667–684. [4] X.P. Wang, H.T. Qi, B. Yu, Z. Xiong, H.Y. Xu, Analytical and numerical study of electroosmotic slip flows of fractional second grade fluids, Commun. Nonlinear. Sci. Numer. Simulat. 50 (2017) 77–87. [5] S.W. Wang, M.L. Zhao, Analytical solution of the transient electro-osmotic flow of a generalized fractional Maxwell fluid in a straight pipe with a circular cross-section, Euro. J. Mech. B/Fluids 54 (2015) 82–86. [6] X. Yang, H.T. Qi, X.Y. Jiang, Numerical analysis for electroosmotic flow of fractional Maxwell fluids, Appl. Math. Lett. 78 (2018) 1–8. [7] C.L. Zhao, C. Yang, Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels, Appl. Math. Comput. 211 (2009) 502–9. [8] C.C. Chang, C.Y. Wang, Rotating electro-osmotic flow over a plate or between two plates, Phys. Rev. E 84 (2011) 056320. [9] S.X. Li, Y.J. Jian, Z.Y. Xie, Q.S. Liu, F.Q. Li, Rotating electro-osmotic flow of third grade fluids between two microparallel plates, Colloids Surf. A 470 (2015) 240–247. [10] P. Kaushik, S. Mandal, S. Chakraborty, Transient electroosmosis of a Maxwell fluid in a rotating microchannel, Electrophoresis 38 (2017) 2741–2748. [11] Z.Y. Xie, Y.J. Jian, Rotating electroosmotic flow of power-law fluids at high zeta potentials, Colloid Surf. A 461 (2014) 231–239. [12] M.Z. Xu, Y.J. Jian, Unsteady rotating electroosmotic flow with time-fractional Caputo-Fabrizio derivative, Appl. Math. Lett. 100 (2020) 106015. [13] A. Mat´ıas, S. S´anchez, F. M´endez, O. Bautista, Influence of slip wall effect on a non-isothermal electro-osmotic flow of a viscoelastic fluid, Int. J. Thermal Sci. 98 (2015) 352–363. [14] G.C. Shit, A. Mondal, A. Sinha, P.K. Kundu, Effects of slip velocity on rotating electro-osmotic flow in a slowly varying micro-channel, Colloid Surf. A 489 (2016) 249–255. [15] S.W. Wang, N. Li, M. Zhao, M.N. Azese, Effects of slip velocity on the rotating electro-osmotic flow of the power-law fluid in a slowly varying microchannel, Z. Naturforsch 73 (2018) 825–831. [16] W. Tan, W. Pan, M. Xu, A note on unsteady flows of a viscoelastic fluid with the fractional Maxwell model between two parallel plates, Int. J. Non-Linear Mech. 38 (2003) 645–650. [17] I.C. Christov, C.I. Christov, Comment on “On a class of exact solutions of the equations of motion of a second grade fluid” by C. Fetecˇau and J. Zierep (Acta Mech. 150, 135–138, 2001), Acta Mech. 215 (2010) 25–28.
6
*Author Contributions Section
Journal Pre-proof Author Contributions Section All the authors contributed equally to this work (conceptualization, methodology, formal analysis, investigation, writing, review & editing) On behalf of myself and my co-authors.
Jo
ur
na lP repr oo f
Haitao Qi
PII: DOI: Reference:
S0893-9659(19)30505-1 https://doi.org/10.1016/j.aml.2019.106179 AML 106179
To appear in:
Applied Mathematics Letters
Received date : 16 October 2019 Revised date : 3 December 2019 Accepted date : 3 December 2019 Please cite this article as: X. Wang, H. Xu and H. Qi, Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids, Applied Mathematics Letters (2019), doi: https://doi.org/10.1016/j.aml.2019.106179. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.
*Manuscript Click here to view linked References
Journal Pre-proof
Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids Xiaoping Wang, Huanying Xu, Haitao Qi∗ School of Mathematics and Statistics, Shandong University, Weihai 264209, PR China
Abstract
na lP repr oo f
In this paper, a rotating electro-osmotic flow of a fractional Maxwell fluid in a parallel plate microchannel with high zeta potentials is examined. The Navier’s slip law at walls is considered. The electric double layer potential distribution is derived by using the nonlinear Poisson-Boltzmann equation. Based on the L1 approximation of the Caputo derivative, a Crank-Nicolson numerical scheme is developed for obtaining the numerical solutions of the rotating electro-osmotic flow velocity profiles. With a purpose to verify the correctness of our numerical results, a comparison has been made with the analytical solutions of the Newtonian fluid given by the previous work and the excellent agreement between the solutions is clear. Finally, the influences of the fractional parameters of α and β, the slip length d and the wall zeta potential ζ on the velocity distribution are also discussed in detail. Keywords: Rotating electro-osmotic flow, Fractional calculus, Finite difference method, Slip boundary condition, High zeta potential
1. Introduction
Jo
ur
In recent years, electro-osmotic flow (EOF) has become a more noticeable topic in microfluidics due to its operational advantages, such as better flow control, high reliability and low noise and has found applications in a wide range of fields, including chemical, biological analyses and medical diagnoses [1]. On the study of EOF, many researchers have concentrated on the non-Newtonian fluid. Among of those, Sarma et al. [2] reported the EOF of the simplified Phan-Thien-Tanner model in a parallel plate microchannel with high zeta potential. Peralta et al. [3] theoretically analyzed the pulsatile EOF of a Maxwell fluid in a parallel flat plate microchannel with asymmetric wall zeta potentials. Recently, the constitutive equations with fractional derivative have been widely used to describe the rheological behaviour of non-Newtonian fluid in the EOF. Wang et al. [4] derived the analytical and numerical solutions of the electro-osmotic slip flow of the fractional second grade fluids. Wang and Zhao [5] and Yang et al. [6] analytically and numerically analyzed the EOF of fractional Maxwell fluids. Zhao and Yang [7] investigated the transient EOF of fractional Oldroyd-B fluids in a rectangular microchannel. On the other hand, microfluidic systems may be placed in a rotating environment in many practical engineering applications, such as in centrifuges for mass separation and flow control. The rotating EOF is considered as an effective technique for mixing different liquids. Thus, the effect of rotating on EOF draws the attention of many researchers. Analytical and numerical works on the rotating EOF of Newtonian and non-Newtonian fluids in a parallel plate microchannel were extensively studied, see Refs. [8–10] and references therein. Based on the nonlinear Poisson-Boltzmann equation, Xie and Jian [11] numerically studied the rotating EOF of power-law fluids in a slit microchannel at high zeta potentials. More recently, Xu and Jian [12] firstly introduced the time fractional CaputoFabrizio derivative to study the rotating EOF. But we found that most of those researches are confined to the EOF in a uniform cross-section channels with constant zeta potential and no slip velocity conditions. However, some experiments and molecular dynamic simulations have demonstrated that the fluid can slip at solid surfaces. And at the level of micro-scales, the velocity slip at the liquid-solid interface plays an important role [13]. By considering Navier’s ∗ Corresponding
author. Tel./Fax: +86 631 5687309. Email address: [email protected] (Haitao Qi)
Preprint submitted to Applied Mathematics Letters
December 3, 2019
Journal Pre-proof
slip condition, Shit et al. [14] and Wang et al. [15] have studied the effects of slip velocity on the rotating EOF in a non-uniform microchannel. But to our best knowledge, the published researches on the rotating EOF of the non-Newtonian fluids is still rare and no attempt has been made regarding the rotating EOF of the viscoelastic fluids with fractional derivative. In this paper, we study the rotating EOF with the fractional Maxwell constitutive relationship in a parallel plate microchannel at high zeta potentials [16] 1 + λα−β Dα−β σ(t) = Eλα Dαt γ(t), 0 ≤ β ≤ α ≤ 1, (1) t where σ(t) is the shear stress, γ(t) is the shear strain, λ = η/E is the relaxation time, in which η is the viscosity constant and E is the shear modulus. Further, α and β are the Caputo fractional derivative parameters. 2. Model description and mathematical formulation
na lP repr oo f
We analyze the rotating EOF of the fractional Maxwell fluids in a parallel plate microchannel. An illustration of the problem and the coordinate system is similar to Figure 1 in the Refs [9, 11]. The origin of the coordinate system is placed in the center of the microchannel. The x axis along the horizontal direction. The y axis is pointing toward the inside of the plane and the z axis is vertical to the parallel plate separating by a distance of 2H. We assume that the microchannel is filled with a liquid solution of dielectric constant ε, and the plates are uniformly charged with a zeta potential ζ that develops an interfacial region of net charge density close to the channel walls, called as the EDL. An external electric field E = (E x , 0, 0) is imposed along the x-axis direction. We consider that the fluid is initially at rest, then the entire system is rotating about the z axis with a constant angular velocity Ω. In a horizontal rotating fluid layer, all the velocity is relative velocity and is considered in non-inertial frame of reference. Because of symmetry, only the upper half domain of the slit microchannel (0 ≤ z ≤ H) is considered. The governing equations for this type incompressible fluid flow are the continuity equation ∇ · V = 0,
(2)
! ∂V + V · ∇V + 2Ω × V = −∇P + ∇ · σ + F, ρ ∂t
(3)
and the general Cauchy momentum equation
where ρ is the fluid density, V = (u, v, 0) is the velocity vector of the rotating EOF flow, Ω = (0, 0, Ω) is the angular velocity vector about z axis and P = p − ρ | Ω × r2 | /2 is the modified fluid pressure due to centrifugal force where r = (x, y, z) and p is pressure [11, 14]. Further, σ is the stress tensor and F = (ρe E x f (t), 0, 0) is the body force vector where f (t) = H(t) is the well-known Heaviside unit step function [6, 17]. From the theory of electrostatics, the net electric charge density ρe in the diffuse layer of the wall EDL and the EDL potential ψ satisfy the following Poisson-Boltzmann equation ! d2 ψ ρe z0 eψ = − , ρe = −2n0 z0 e sinh . (4) ε κb T dz2
ur
Here n0 is the bulk concentration of the ions, z0 is the valence of the ions, e is the elementary electronic charge, ε is the dielectric constant of the fluid, kB is the Boltzmann constant and T is the absolute temperature. The boundary conditions for the potential given in the following form dψ = 0, ψ(H) = ζ. (5) dz z=0
Jo
Then the analytical expression of the potential distribution in the dimensionless form can be obtained as [2] ζ ψ = 4 tanh−1 tanh eK(z−1) , 4
(6)
by introducing the following non-dimensional variables
ψ ζ z ψ¯ = , ζ¯ = , z¯ = , K = κH. ψz ψz H 2
(7)
Journal Pre-proof
Here ψz = kb T/z0 e, κ = (2n0 z20 e2 /(εkb T ))1/2 is the reciprocal of the EDL thickness (λ = κ−1 ). In view of the above-mentioned assumptions and in the absence of the modified pressure gradient, the continuity Eq. (2) is satisfied automatically. For this specific problem, the constitutive equation for a fractional Maxwell fluid (1) can be expressed as ! ! α−β α−β α α−1 ∂v α α−1 ∂u , 1 + λ D σ = Eλ D , (8) 1 + λα−β Dα−β σ = Eλ D zy zx t t t t ∂z ∂z and the general Cauchy momentum equation (3) is simplified to ! ! ∂σzy ∂σzx ∂v ∂u − 2Ωv = + ρe E x f (t), ρ + 2Ωu = , ρ ∂t ∂z ∂t ∂z
(9)
which is subject to the following initial and boundary conditions u = 0,
v = 0,
∂u ∂v = 0, = 0, at t = 0, ∂t ∂t
(10)
na lP repr oo f
∂v ∂u = 0, = 0, at z = 0, (11) ∂z ∂z ∂v ∂u = 0, v + d = 0, at z = H. (12) u+d ∂z ∂z where d is the slip length. Considering that the bottom and upper walls of the microchannel are mostly made of the same material leading to the same hydrophobic characteristics and the same slip boundary condition. So we use the same slip length. Substituting Eq. (8) into Eq. (9) and introducing the following dimensionless variables u¯ =
u v t λ z d , v¯ = , t¯ = , λ¯ = , z¯ = , d¯ = , US US 1/Ω 1/Ω H H
(13)
where US = −εζE x /η is the Helmholtz-Smoluchowski velocity and η = λE, we obtain the dimensionless governing equations as: ! ! 2 ∂u K2 α−1 α−1 ∂ u (14) − 2v = λ D + sinh ψ 1 + λα−β Dtα−β f (t), Re 1 + λα−β Dα−β t t 2 ∂t ζ ∂z ! ! 2 ∂v α−β α−β α−1 α−1 ∂ v Re 1 + λ Dt . (15) + 2u = λ Dt ∂t ∂z2 Here, Re = ρΩH 2 /η is the rotating Reynolds number, f (t) = H(t) is the well-known Heaviside unit step function, which can be approximated by 1/(1 + exp((−t + p)/q)) with p → 0 and q → 0 in the following analysis. By making some transformations, the equivalent form of Eqs. (14) and (15) can be written as
ur
2 K2 1−α α−1 ∂ u 1−α α−β 1−β Re D2−α u + λα−β D2−β v − 2λα−β D1−β + sinh ψ D + λ D f (t), t t t u − 2Dt t v =λ t ζ ∂z2
Jo
2 1−α α−1 ∂ v Re D2−α v + λα−β D2−β u + 2λα−β D1−β , t t v + 2Dt t u =λ ∂z2 which has the following dimensionless initial and boundary conditions:
u = 0,
v = 0,
∂u ∂v = 0, = 0, at t = 0, ∂t ∂t
∂u ∂v = 0, = 0, at z = 0, ∂z ∂z ∂u ∂v u+d = 0, v + d = 0, at z = 1. ∂z ∂z 3
(16) (17)
(18) (19) (20)
Journal Pre-proof
3. Numerical algorithms In this section, we are using a finite difference method to solve the Eqs. (16)-(20). We first introduce a uniform discretization over the rectangular computational domain (D = {(z, t) | z, t ∈ [0, 1] × [0, Tˆ ]}) given by zi = ih, i = 0, 1, 2, . . . M; tn = nτ, n = 0, 1, 2, . . . N, in which h = 1/M and τ = Tˆ /N are the spatial and temporal step sizes, respectively. We define uni and vni as the numerical solutions of u(z, t) and v(z, t) at the mesh point (zi , tn ) Base on the L1 approximations of the Caputo fractional derivative, the Crank-Nicolson numerical schemes of Eqs. (16) and (17) can be obtained as n−1 n−1 k−1/2 Reλα−β τβ−1 (β) n−1/2 X k−1/2 Reτα−1 (α) n−1/2 X (α) (β) (β) (α) + b δt ui b δt ui − bn−k−1 − bn−k δt ui − bn−k−1 − bn−k δt ui Γ(1 + α) 0 Γ(1 + β) 0 k=1 k=1 n−1 n−2 X Reτα−1 (α) Reτα−1 (α) n X (α) (α) (α) (α) k k n−1 b v − b v − − b − bn−k vi − bn−k−2 − bn−k−1 vi Γ(1 + α) 0 i k=1 n−k−1 Γ(1 + α) 0 i k=1 n−1 n−2 α−β β−1 α−β β−1 X Reλ τ (β) n−1 X Reλ τ (β) n (β) (β) (β) (β) k k b v − b v − − b − bn−k vi − bn−k−2 − bn−k−1 vi Γ(1 + β) 0 i k=1 n−k−1 Γ(1 + β) 0 i k=1
na lP repr oo f
n−1 n−1 uni+1 − 2uni + uni−1 + un−1 + ui−1 K 2 sinh ψ 1−α n λα−β K 2 sinh ψ 1−β n i+1 − 2ui n−1 + D Dt ( f + f n−1 ), ( f + f ) + t 2ζ 2ζ 2h2 i = 1, 2, 3, . . . , M − 1, n = 1, 2, 3, . . . , N, (21)
= λα−1
n−1 n−1 k−1/2 k−1/2 Reλα−β τβ−1 (β) n−1/2 X Reτα−1 (α) n−1/2 X (α) (β) (β) (α) + b δt v b δt vi − bn−k−1 − bn−k δt vi − bn−k−1 − bn−k δt vi i Γ(1 + α) 0 Γ(1 + β) 0 k=1 k=1 n−1 n−2 X Reτα−1 (α) Reτα−1 (α) n X (α) n−1 k b0 ui − + b u − bn−k−1 − b(α) u + b(α) − b(α) uki 0 i i n−k n−k−2 n−k−1 Γ(1 + α) Γ(1 + α) k=1 k=1 n−1 n−2 α−β β−1 α−β β−1 X Reλ τ (β) n−1 X Reλ τ (β) n (β) (β) (β) (β) k k b u − b u − bn−k−1 − bn−k ui + bn−k−2 − bn−k−1 ui + Γ(1 + β) 0 i Γ(1 + β) 0 i k=1
n−1 vn − 2vni + vni−1 + vn−1 + vn−1 i+1 − 2vi i−1 = λα−1 i+1 , 2h2
i = 1, 2, 3, . . . , M − 1,
k=1
n = 1, 2, 3, . . . , N,
(22)
where bα0 = bβ0 = 1, bαj = ( j + 1)α − jα , bβj = ( j + 1)β − jβ , j = 0, 1, 2, .... The discretized initial and boundary conditions are u0i = 0, v0i = 0,
u1i − u0i = 0, τ
v1i − v0i = 0, i = 0, 1, 2, . . . , M, τ
Jo
ur
un1 − un0 vn1 − vn0 = 0, = 0, n = 0, 1, 2, . . . , N, h h un − unM−1 vn − vnM−1 unM + d M = 0, vnM + d M = 0, n = 0, 1, 2, . . . , N. h h
(23) (24) (25)
4. Results and discussion
In the present work, all the parameters used in the governing equations are dimensionless. In order to verify the accuracy of the numerical method and the theoretical analysis in our study. We compare the obtained numerical solutions of the velocity distribution by taking α = 1, β = 0, λ = 0, d = 0, K = 10, ζ = 1, Re = 1 and t = 5 with the analytical solutions for the rotating EOF of the Newtonian fluid with the no slip boundary condition presented by Chang and Wang [8]. Results are shown in Fig. 1. It can be seen that the numerical results are agreed well with the 4
Journal Pre-proof
analytical solutions. This demonstrates that the numerical method provided in Sect. 3 is effective and feasible and can be used to compute the rotating EOF velocity profiles of the fractional Maxwell fluid. To analyze the influence of the fractional parameters α and β, the slip length d and the wall zeta potential ζ on the velocity distribution. The values of the dimensionless parameters for the EOF in a microchannel are chosen as follows: the electro-kinetic width K=10 and the Reynolds number Re = 1 [8, 14], the relaxation time λ = 0.5 [3]. Fig. 2 illustrates the variations of the rotating EOF velocity for different values of the fractional parameter α (α=0.7, 1) and the slip length d (d = 0 [8, 9, 11]: for non-slip flow and d = 0.003, 0.008 [13]: for slip flow) at t = 5. As shown in the Fig. 2(a), the axis velocity magnitude increases with α and the presence of slip velocity has an effect on the rotating EOF velocity. We can also see that the velocity profile is steep in a region near to the microchannel wall while maintains a constant value near the central axis of the microchannel. Fig. 3 presents the variations of the rotating EOF velocity for different values of the fractional parameter β (β=0.3, 0.6) and the wall zeta potential ζ (ζ = 1 [3, 8, 9]: for low zeta potential and ζ = 3, 5 [2]: for high zeta potential) at t = 5. It can be seen from the Fig. 3(a) that increasing the values of β tends to sharply enhance the velocity distribution in x direction near the wall and the peaks of velocity profiles under the high zeta potential are higher than that the case of low zeta potential. However, the velocity in y direction becomes more and more negative as the increasing α, β, d and ζ. 0
0.7
na lP repr oo f
0.8
Present numerical solution Analytical solution of Chang and Wang
-0.1
0.6
-0.2
v
u
0.5 0.4 0.3
-0.3 -0.4
0.2
Present numerical solution Analytical solution of Chang and Wang
0.1
-0.5 -0.6
0 0
0.2
0.4
z
0.6
0.8
0
1
(a) Dimensionless velocity u in x direction
0.2
0.4
0.6
0.8
1
z
(b) Dimensionless velocity v in y direction
Figure 1: Comparison of our present numerical solution with the previous analytical results presented by Chang and Wang [8].
0
1
-0.1
0.8
-0.2
0.6
α=0.7
-0.3
v
u
α=1
-0.4
0.4
-0.5
α=0.7
0.2
-0.6
α=1
-0.7
0 0
0.2
0.4
0.6
0.8
0
1
ur
z
(a) Dimensionless velocity u in x direction
0.2
0.4
0.6
0.8
1
z
(b) Dimensionless velocity v in y direction
Jo
Figure 2: Variation of the velocity distribution for different values of α and d with β = 0.3, ζ = 3 at t = 5 (d = 0 solid lines, d = 0.003 dashed lines, d = 0.008 dash-dot lines).
5. Conclusions
In this paper, we study the rotating EOF of fractional Maxwell fluid through a parallel plate microchannel. The finite difference method is used to solve the governing equations of the velocity distribution. The effectiveness of our method is proved by giving a comparison with the known analytical results of the Newtonian fluid. The effects of 5
Journal Pre-proof
0.7
0 β=0.6
0.6
β=0.3
-0.1 0.5
-0.2
u
v
0.4 0.3
-0.3
β=0.3 0.2
-0.4
β=0.6
0.1
-0.5
0 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
0.6
0.8
1
z
z
(a) Dimensionless velocity u in x direction
(b) Dimensionless velocity v in y direction
Figure 3: Variation of the velocity distribution for different values of β and ζ with α = 0.8, d = 0.008 at t = 5 (ζ = 1 solid lines, ζ = 3 dashed lines, ζ = 5 dash-dot lines).
Acknowledgements
na lP repr oo f
the fractional parameters α and β, the slip length d and the wall zeta potential ζ on the velocity profiles are analyzed graphically. Results indicate that the axis velocity profiles increase with the fractional parameters α and β. The presence of slip velocity can enhance the magnitude of the axis velocity distribution and the peaks of the axis velocity profiles under the high zeta potential are higher than that the case of low zeta potential. The velocity in y direction becomes more and more negative as the increasing α, β, d and ζ.
This work was supported by the National Natural Science Foundation of China (NSFC No. 11672163, No. 11771254), the Fundamental Research Funds for the Central Universities (No. 2019ZRJC002). References
Jo
ur
[1] H. Bruus, Theoretical Microfluidics, Oxford University Press, New York, 2008. [2] R. Sarma, N. Deka, K. Sarma, P.K. Mondal, Electroosmotic flow of Phan-Thien-Tanner fluids at high zeta potentials: An exact analytical solution, Phys. Fluids 30 (2018) 062001. [3] M. Peralta, O. Bautista, F. M´endez, E. Bautista, Pulsatile electroosmotic flow of a Maxwell fluid in a parallel flat plate microchannel with asymmetric zeta potentials, Appl. Math. Mech. -Engl. Ed. 39 (2018) 667–684. [4] X.P. Wang, H.T. Qi, B. Yu, Z. Xiong, H.Y. Xu, Analytical and numerical study of electroosmotic slip flows of fractional second grade fluids, Commun. Nonlinear. Sci. Numer. Simulat. 50 (2017) 77–87. [5] S.W. Wang, M.L. Zhao, Analytical solution of the transient electro-osmotic flow of a generalized fractional Maxwell fluid in a straight pipe with a circular cross-section, Euro. J. Mech. B/Fluids 54 (2015) 82–86. [6] X. Yang, H.T. Qi, X.Y. Jiang, Numerical analysis for electroosmotic flow of fractional Maxwell fluids, Appl. Math. Lett. 78 (2018) 1–8. [7] C.L. Zhao, C. Yang, Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels, Appl. Math. Comput. 211 (2009) 502–9. [8] C.C. Chang, C.Y. Wang, Rotating electro-osmotic flow over a plate or between two plates, Phys. Rev. E 84 (2011) 056320. [9] S.X. Li, Y.J. Jian, Z.Y. Xie, Q.S. Liu, F.Q. Li, Rotating electro-osmotic flow of third grade fluids between two microparallel plates, Colloids Surf. A 470 (2015) 240–247. [10] P. Kaushik, S. Mandal, S. Chakraborty, Transient electroosmosis of a Maxwell fluid in a rotating microchannel, Electrophoresis 38 (2017) 2741–2748. [11] Z.Y. Xie, Y.J. Jian, Rotating electroosmotic flow of power-law fluids at high zeta potentials, Colloid Surf. A 461 (2014) 231–239. [12] M.Z. Xu, Y.J. Jian, Unsteady rotating electroosmotic flow with time-fractional Caputo-Fabrizio derivative, Appl. Math. Lett. 100 (2020) 106015. [13] A. Mat´ıas, S. S´anchez, F. M´endez, O. Bautista, Influence of slip wall effect on a non-isothermal electro-osmotic flow of a viscoelastic fluid, Int. J. Thermal Sci. 98 (2015) 352–363. [14] G.C. Shit, A. Mondal, A. Sinha, P.K. Kundu, Effects of slip velocity on rotating electro-osmotic flow in a slowly varying micro-channel, Colloid Surf. A 489 (2016) 249–255. [15] S.W. Wang, N. Li, M. Zhao, M.N. Azese, Effects of slip velocity on the rotating electro-osmotic flow of the power-law fluid in a slowly varying microchannel, Z. Naturforsch 73 (2018) 825–831. [16] W. Tan, W. Pan, M. Xu, A note on unsteady flows of a viscoelastic fluid with the fractional Maxwell model between two parallel plates, Int. J. Non-Linear Mech. 38 (2003) 645–650. [17] I.C. Christov, C.I. Christov, Comment on “On a class of exact solutions of the equations of motion of a second grade fluid” by C. Fetecˇau and J. Zierep (Acta Mech. 150, 135–138, 2001), Acta Mech. 215 (2010) 25–28.
6
*Author Contributions Section
Journal Pre-proof Author Contributions Section All the authors contributed equally to this work (conceptualization, methodology, formal analysis, investigation, writing, review & editing) On behalf of myself and my co-authors.
Jo
ur
na lP repr oo f
Haitao Qi