Two-Dimensional Electromagnetohydrodynamic (EMHD) Flows of Fractional Viscoelastic Fluids with Electrokinetic Effects

Abstract

The present study provides analytical and numerical solutions for an electromagnetohydrodynamic (EMHD) flow using a Caputo time-fractional Maxwell model. The flow is a typical rectangular channel flow. When the scale of the cross-stream is much smaller than the streamwise and spanwise scales, the model is approximated as a two-dimensional slit parallel plate flow. Moreover, the influence of the electric double layer (EDL) at the solid–liquid interface is also considered. The electro-osmotic force generated by the interaction between the electric field and the EDL will induce a flow (i.e., electro-osmotic flow). Due to the application of the electric field at the streamwise and the vertical magnetic field, the flow is driven by Lorentz force along the spanwise direction. Simultaneously, under the action of the magnetic field, the electro-osmotic flow induces a reverse Lorentz force, which inhibits the electro-osmotic flow. The result shows that resonance behavior can be found in both directions in which the flow is generated. However, compared with the classical Maxwell fluid, the slip velocity and resonance behavior of fractional Maxwell fluid are suppressed. In the spanwise direction, increasing the strength of magnetic field first promotes the slip velocity and resonance behavior, and then suppresses them, while in the streamwise direction, both the electro-osmotic flow and resonance behavior are suppressed with the magnetic field.

1. Introduction

Microfluidic devices are greatly utilized in many fields, such as biological engineering, the transport of chemicals in the body and heat transfer in electronic components [1,2]. Microfluidic transport can be realized by pressure-driven micropumps [3,4], electro-osmotic micropumps [5,6,7,8], and electromagnetohydrodynamic (EMHD) pumps [9,10,11]. The EMHD micropump is driven by the Lorenz force, which is produced by the interaction of an external electric field and a magnetic field [12]. Compared with other types of micropumps, the EMHD micropump has several advantages, namely, a simple manufacturing process, continuous flow power and two-way pumping capacity. For instance, it can be used in fluid pumping, flow control in fluidic networks and fluid stirring and mixing [10]. In addition, the EMHD micropump can reduce energy consumption and save costs in industry, so related areas of study have always attracted attention from researchers, such as rotating EMHD pumps [11], stagnation point flow [13], EMHD flow in corrugated walls [14], etc. Recently, Khan et al. [15,16] considered the impact of gyrotactic microorganisms on the nonlinear mixed convective MHD flow of thixotropic and Prandtl–Eyring nanofluids. They examined the variations of heat transfer characteristics subject to nonlinear radiative flux and heat source/sink, and found that for larger thixotropic fluid parameters, the velocity field boosts up, while for rising values of the Hartmann number, the velocity and temperature have opposite behaviors. Furthermore, at the channel walls, most solids spontaneously acquire surface electric charges when brought into contact with a polar medium, and then the number of positive and negative ions in the solution near the solid–liquid interface is inconsistent, forming an electric double layer (EDL). The electro-osmotic force generated by the interaction between the electric field and the EDL will move the fluid along the electric field. In EMHD flow, however, the influence of the electro-osmotic force is often ignored (such as Refs. [10,11,12,13,14,17,18]). Therefore, the coupling effect between EMHD flow and electro-osmotic flow needs to be further investigated, which is one of the intentions of this paper.

Viscoelastic fluids (e.g., blood and polymer solutions) are often used in microfluidic transports. Since the fluids with both elastic and viscous properties, in the flow field, deform/flow under external force, the external force is vanished, and the deformation will return to the specified threshold with time evolution [19]. Due to this characteristic, a viscoelastic fluid often experiences resonance during the flow process, i.e., it increases abruptly the oscillation amplitude of a system when an external force matches the system’s natural frequency [20]. The resonance behavior with viscoelastic fluid flow has been studied by many researchers. Yin and Zhu [21] studied the unidirectional oscillating flow produced by a periodic pressure gradient through the fractional Maxwell model. Andrienko et al. [22] investigated resonance-like phenomena in axisymmetric Poiseuille flows of viscoelastic fluids. Lambert et al. [23] analyzed the heat transfer enhancement in an oscillatory flow of a viscoelastic fluid in tubes. Tsiklauri and Beresnev [24] found sharp enhancements of the flow through researching the process of transition from a dissipative regime to an elastic regime with Maxwell fluids. Emilio Herrera [25] predicted flow enhancement on the pulsating flow in an aqueous worm-like micellar solution of cetyltrimethyl ammonium tosilate for various concentrations. In two-dimensional EMHD flow, will the viscoelastic fluid also have a similar flow enhancement phenomenon? Another purpose of this paper is to describe the resonance behavior during flow using a viscoelastic model.

In recent years, fractional calculus has been widely used in the research of abnormal diffusion, viscoelastic model, and soft materials [21]. Compared with the integer calculus, fractional calculus can more concisely and accurately describe the physical process with historical memory and spatial correlation [26]. The germination of fractional calculus can be traced back to L’ Hospital’s letter to Leibniz in 1695. After a series of work of many researchers, including Fourier, Abel, Liouville, and Caputo [27], the embryonic form of fractional calculus came out. Nowadays, fractional calculus has developed into a systematic branch of mathematics, and it is widely used in various fields [28]. Feng et al. [29] used the time-distributed Caputo fractional model to study the rotating electro-osmotic slip flow of generalized Maxwell fluids in a periodic electric field. Cao et al. [30] numerically investigated the electro-osmotic flow of double layers consisting of Newtonian fluid and fractional second-order fluid with a rotating frame in a parallel microchannel. Yang et al. [31] studied the electro-osmotic flow of Maxwell fluids with Riemann–Liouville fractional derivatives in a rectangular microchannel, and proposed the nonlinear conjugate gradient method to obtain the viscoelastic parameters. Abdulhameed et al. [32] explored unsteady pressure-driven and EMHD flow of an incompressible Maxwell fluid with time-fractional Caputo–Fabrizio derivatives through a circular tube. Abro and Solangi [33] analytically researched the heat transfer and free convection problem in second-grade fluid with porous impacts using Caputo–Fabrizoi fractional derivatives. Liu et al. [34] studied the unsteady EMHD flow of fractional Oldroyd-B fluids between parallel plates on heat transfer. Inspired by these studies, the fractional Maxwell model was chosen to explore the characteristics of viscoelastic fluid flow.

Furthermore, the velocity slip on the channel wall is also an important feature in micro/nano-scale flow [35,36,37]. In macro-fluid flow, the no-slip boundary condition is applied successfully to model some experiments. This success may not always reflect the accuracy of the boundary condition but may reflect the insensitivity of the experiment to the boundary condition [38]. At the microscale level, slip boundary conditions will become important when the length scale over which the liquid velocity changes approaches the slip length in a channel [39], where the slip length is the distance inside the walls at which the fluid velocity extrapolates to zero. Pascall and Squires [40] measured the induced charge of electro-osmosis over gold electrodes; both the magnitude and frequency dependence of the measured slip velocity are captured quantitatively by accounting for the physical capacitance and surface chemistry of the dielectric layer. Galea and Attard [41] presented a model of atomic roughness to study the effect of solid roughness on the slip boundary condition during shear flow. By using evanescent waves, Bouzigues [42] measured velocity profiles and diffusion profiles in pressure-driven and electro-osmotic flows, and determine the hydrodynamic slip lengths with 10 nm accuracy in the Debye layer for hydrophilic and hydrophobic surfaces. Khair and Squires [43] found that an enhancement occurs in the electrophoretic motion of a colloidal particle whose surface exhibits hydrodynamic slip.

As already mentioned, viscoelastic fluids have resonance phenomena in a tube/channel flow. Whether the resonance behavior occurs at the micro-nano scale is exactly what we want to know. At the micro-nano scale, electrostatic force is particularly important, so on the basis of previous research [1,12,17,44], we extended the research of Jian [44] to fractional viscoelastic fluids to study the resonance behavior in channel flow. For this purpose, current paper investigates the EMHD flow of fractional viscoelastic fluids with electro-osmosis and velocity slip effects. The special feature of this model is that the Lorentz force has a component in the direction of the electric field. The scale of the cross-stream (height) is much smaller than the streamwise (length) and spanwise (width) scales (i.e., approximately an infinite slit plate). The physical mechanism is as follows. Due to the application of the electric field at the streamwise direction and the vertical magnetic field, the flow is driven by Lorentz force, which originates from the interaction between an imposed periodic electric field and the magnetic field, and its direction is along the spanwise direction. Moreover, the influence of the electric double layer (EDL) at the solid–liquid interface is also considered. The electro-osmotic force generated by the interaction between the electric field and the EDL will also induce a flow (i.e., electro-osmotic flow). Therefore, for such a two-dimensional channel flow, the Fourier transform can be used to obtain an analytical solution. Meanwhile, a numerical algorithm can also be obtained by the finite difference method. The effects of several dimensionless numbers, such as the Hartmann number $Ha$, fractional parameters $\alpha$, and slip length L, on the velocity and resonance behavior are analyzed graphically in detail. The result shows that resonance behavior can be found in both directions in which the flow is generated. However, compared with the classical Maxwell fluid, the slip velocity and resonance behavior of fractional Maxwell fluid are suppressed. In the spanwise direction, increasing the strength of magnetic field first promotes the slip velocity and resonance behavior, and then suppresses them. In the direction of the electric field, both electro-osmotic flow and resonance behavior are suppressed with the magnetic field. The rest of this paper is arranged as follows. In Section 2, the problem is formulated, and solutions of the EMHD velocity are presented. In Section 2.1, the net charge density is given by solving the linearized Poisson–Boltzmann equation. The Caputo fractional derivative and the governing equation are presented in Section 2.2. In Section 2.3, the Fourier transform method is used to obtain velocity expression. A finite difference scheme is provided for the velocity distribution in Section 2.4. Velocity distribution and the resonance of fractional viscoelastic fluid are investigated in Section 3.1 and Section 3.2, respectively. Finally, the study is summarized and concluded in Section 4.

2. Mathematical Modeling

In this paper, we studied the slip flow of EMHD pump combined with the EDL effect in the micro-parallel channel. The model sketch is shown in Figure 1. The electrical field $E\left(t\right)$ and a steady magnetic field of strength B are applied between two negatively charged micro-parallel plates separated by a distance $2H$ in the y-axis (streamwise) and z-axis (cross-stream) directions, respectively. The flow is a typical rectangular channel flow, as indicated in Section 1, assuming that the scale of the cross-stream is much smaller than the streamwise and spanwise scales. In this case, the model can be approximated as a two-dimensional slit-parallel plate flow. It can be seen from Figure 1a that the forces acting on the fluid during the entire flow are only electro-osmotic force and Lorentz force, and the velocity in the x-axis and y-axis directions are functions of z and t. Thus, the velocity vector can be expressed as $\mathbf{U}=\left(u(z,t),v(z,t),0\right)$. The electro-osmotic force acting on the fluid is ${\rho }_{e}E$ along the y-axis direction. The Lorentz force $\mathbf{F}=\mathbf{J}\times \mathbf{B}=\sigma B\left(E-Bu,-Bv,0\right)$ has two components, one of which is perpendicular to the electric and magnetic fields, and the other is parallel to the direction of the electric field, where the local ion current density is $\mathbf{J}=\sigma \left(\mathbf{E}+\mathbf{U}\times \mathbf{B}\right)=\sigma \left(Bv,E-Bu,0\right)$. For convenience, we shall denote $F_1=\sigma B\left(E-Bu\right)$ and $F_2=-\sigma B^{2}v$, respectively.

For the walls at $z=\pm H$, the slip boundary condition will be used. The so-called slip boundary condition means that on the solid wall (on the boundary), the velocity of the fluid is different from that of the solid wall, that is, the slip velocity is generated (see Figure 1b). Here, the partial slip boundary condition is used, which means that the tangential velocity of the boundary has a certain gradient. Navier pointed out that the velocity of the fluid at the solid wall is proportional to the gradient of the fluid velocity along the normal to the boundary surface, where the proportionality factor is called the slip length. The slip length l is the distance inside the wall, where the fluid velocity would extrapolate to zero [45].

2.1. The Local Net Charge Density and the Electrical Potential

The local net charge density ${\rho }_e$ and the electrical potential $\psi \left(z\right)$ are described by the following Poisson–Boltzmann equations [19,46]:

$$\begin{array}{cc}\hfill & {\nabla }^2\psi \left(z\right)=-\frac{{\rho }_e\left(z\right)}{ϵ},\hfill \\ \hfill & {\rho }_e\left(z\right)=z_{\nu }e\left(n_{+}-n_{-}\right)=-2n_0{z}_{\nu }esinh\left[\left(z_{\nu }e\psi \right)/\left(k_{b}T\right)\right],\hfill \end{array}$$

(1)

where $n_{\pm }=n_{0}exp\left[\mp \left(e{z}_{\nu }\psi \right)/\left(k_{b}T\right)\right]$ are the ionic number concentrations for the cations and anions in the liquid, respectively, $ϵ$ is the dielectric coefficient of the electrolyte liquid, e is the electron charge, $z_{\nu }$ is the valence, $k_b$ is the Boltzmann constant, $n_0$ is the ion density of bulk liquid, and T is the absolute temperature. When the potential $\psi \left(z\right)$ is low, the Debye–Hückel approximation can be applied to obtain the linearized equation [19,46]

$$\frac{{\mathrm{d}}^2\psi }{\mathrm{d}{z}^2}=k^2\psi ,\mathrm{with}{k}^2=\frac{2n_0{z}_{\nu }^2{e}^2}{ϵk_{B}T}.$$

(2)

The boundary conditions of potential are given as follows:

$$\psi ={\psi }_0,\mathrm{at}z=\pm H,$$

(3)

where ${\psi }_0$ is the potential at walls. The solution of Equation (2) subject to the condition Equation (3) is given by

$$\psi \left(z\right)={\psi }_0\frac{cosh\left(kz\right)}{cosh\left(kh\right)},-H\le z\le H.$$

(4)

Using the above potential distribution, the net charge density can be expressed as

$${\rho }_e\left(z\right)=-ϵk^2{\psi }_0\frac{cosh\left(kz\right)}{cosh\left(kH\right)}.$$

(5)

2.2. Caputo Fractional Derivative and the Governing Equation

The Caputo fractional derivative is defined as [29]

$$D_t^{p}f\left(t\right)=\left\{\begin{array}{cc}\frac{1}{\Gamma (n-p)}{\int }_0^t\frac{f^{\left(n\right)}\left(s\right)}{{(t-s)}^{p-n+1}}\mathrm{d}s,\hfill & n-1<p<n,\hfill \\ \frac{{\mathrm{d}}^{n}f\left(t\right)}{\mathrm{d}{t}^n},\hfill & p=n\in N,\hfill \end{array}\right.$$

(6)

where p is the fractional derivative parameter, and $\Gamma (·)$ is the Gamma function.

The Navier–Stokes equation is

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\rho \left(\mathbf{u}·\nabla \right)\mathbf{u}=-\nabla P+\nabla ·\tau +\mathbf{f},$$

(7)

where $\rho$ is the fluid density, P is the pressure, $\tau$ is the stress tensor, and $\mathbf{f}$ is the external body force vector.

The fractional Maxwell constitutive equation is given by

$$\left(1+{\lambda }^{\alpha }{D}_t^{\alpha }\right)\tau =\mu {\lambda }^{\beta -1}{D}_t^{\beta -1}\dot{\gamma },\mathrm{with}0\le \alpha \le \beta \le 1,$$

(8)

where $\dot{\gamma }$ is the shear rate, $\lambda$ is the relaxation time of the fluid, $\mu$ is a dynamic viscosity, and $\alpha$ and $\beta$ are the fractional derivative parameters. When $\alpha =\beta =1$, it reduces to the ordinary Maxwell model, while for $\alpha =0$ and $\beta =1$, it is simplified as the classical Newtonian fluid [47].

In the case of low Reynolds number flow, the inertia term can be ignored. Based on the assumption that the pressure gradient between the plates is zero, the simplified governing equations are obtained:

$$\left\{\begin{array}{c}\begin{array}{c}\rho \frac{\partial u}{\partial t}=\frac{\partial {\tau }_{zx}}{\partial z}+\sigma B\left(E-Bu\right),\hfill \\ \rho \frac{\partial v}{\partial t}=\frac{\partial {\tau }_{zy}}{\partial z}+{\rho }_{e}E-\sigma B^{2}v.\hfill \end{array}\hfill \end{array}\right.$$

(9)

The simplified constitutive equation is

$$\left\{\begin{array}{c}\begin{array}{c}\left(1+{\lambda }^{\alpha }{D}_t^{\alpha }\right){\tau }_{zx}=\mu {\lambda }^{\beta -1}{D}_t^{\beta -1}\left(\frac{\partial u}{\partial z}\right),\hfill \\ \left(1+{\lambda }^{\alpha }{D}_t^{\alpha }\right){\tau }_{zy}=\mu {\lambda }^{\beta -1}{D}_t^{\beta -1}\left(\frac{\partial v}{\partial z}\right).\hfill \end{array}\hfill \end{array}\right.$$

(10)

Substituting Equation (9) into Equation (10), we obtain

$$\left\{\begin{array}{c}\begin{array}{c}\left(1+{\lambda }^{\alpha }{D}_t^{\alpha }\right)\left[\rho \frac{\partial u}{\partial t}-\sigma B\left(E-Bu\right)\right]=\mu {\lambda }^{\beta -1}{D}_t^{\beta -1}\left(\frac{{\partial }^{2}u}{\partial z^2}\right),\hfill \\ \left(1+{\lambda }^{\alpha }{D}_t^{\alpha }\right)\left[\rho \frac{\partial v}{\partial t}-\left({\rho }_{e}E-\sigma B^{2}v\right)\right]=\mu {\lambda }^{\beta -1}{D}_t^{\beta -1}\left(\frac{{\partial }^{2}v}{\partial z^2}\right).\hfill \end{array}\hfill \end{array}\right.$$

(11)

The slip boundary conditions are

$$\left\{\begin{array}{c}\begin{array}{c}{\left[u\pm l\frac{\partial u}{\partial z}\right]}_{z=\pm 1}=0,\hfill \\ {\left[v\pm l\frac{\partial v}{\partial z}\right]}_{z=\pm 1}=0.\hfill \end{array}\hfill \end{array}\right.$$

(12)

Introduce a set of dimensionless parameters as follows:

$$\begin{array}{cc}\hfill & \overline{z}=\frac{z}{H},\overline{u}=\frac{u}{U_{eo}},\overline{v}=\frac{v}{U_{eo}},K=kH,\overline{t}=\frac{t\mu }{\rho H^2},L=\frac{l}{H},G=\sqrt{\frac{\sigma }{\mu }}\frac{E_{0}H}{U_{eo}},\hfill \\ \hfill & Ha=BH\sqrt{\frac{\sigma }{\mu }},De=\frac{\lambda \mu }{\rho H^2},\overline{E}=\frac{E}{E_0}\mathrm{with}{U}_{eo}=-\frac{ϵE_0{\psi }_0}{\mu },\hfill \end{array}$$

(13)

where $E_0$ is the maximum value of the electric field E. The dimensionless forms of Equations (11) and (12) are as follows:

$$\left\{\begin{array}{c}\begin{array}{c}\left(1+D{e}^{\alpha }{D}_{\overline{t}}^{\alpha }\right)\left\{\frac{\partial \overline{u}}{\partial \overline{t}}-\left[HaG\overline{E}-H{a}^2\overline{u}\right]\right\}=D{e}^{\beta -1}{D}_{\overline{t}}^{\beta -1}\left(\frac{{\partial }^2\overline{u}}{\partial {\overline{z}}^2}\right),\hfill \\ \left(1+D{e}^{\alpha }{D}_{\overline{t}}^{\alpha }\right)\left\{\frac{\partial \overline{v}}{\partial \overline{t}}-\left[K^2\frac{cosh\left(K\overline{z}\right)}{cosh\left(K\right)}\overline{E}-H{a}^2\overline{v}\right]\right\}=D{e}^{\beta -1}{D}_t^{\beta -1}\left(\frac{{\partial }^2\overline{v}}{\partial {\overline{z}}^2}\right),\hfill \end{array}\hfill \end{array}\right.$$

(14)

$$\left\{\begin{array}{c}\begin{array}{c}{\left[\tilde{u}\pm L\frac{\partial \tilde{u}}{\partial \overline{z}}\right]}_{\overline{z}=\pm 1}=0,\hfill \\ {\left[\tilde{v}\pm L\frac{\partial \tilde{v}}{\partial \overline{z}}\right]}_{\overline{z}=\pm 1}=0.\hfill \end{array}\hfill \end{array}\right.$$

(15)

2.3. Analytical Solutions of the EMHD Velocity Field

In order to obtain the analytical solution, we introduce the Fourier transform

$$\tilde{X}\left(s\right)={\int }_{-\infty }^{+\infty }X\left(t\right)e^{-ist}\mathrm{dt},$$

(16)

where s is a real parameter. It should be mentioned here that the use of Fourier transform or Laplace transform to solve fractional differential equations is very mature. The analytical solution obtained by using the method of similarity transformation is currently the most common method for dealing with fractional differential equations [48]. Then, using the Fourier transform with respect to $\overline{t}$ for the Equation (14), we obtain the following equation:

$$\left\{\begin{array}{c}\begin{array}{c}\left[1+{\left(is\right)}^{\alpha }D{e}^{\alpha }\right]\left(is\tilde{u}-HaG\tilde{E}+H{a}^2\tilde{u}\right)={\left(is\right)}^{\beta -1}D{e}^{\beta -1}\frac{{\mathrm{d}}^2\tilde{u}}{\mathrm{d}{\overline{z}}^2},\hfill \\ \left[1+{\left(is\right)}^{\alpha }D{e}^{\alpha }\right]\left(is\tilde{v}-K^2\frac{cosh\left(K\overline{z}\right)}{cosh\left(K\right)}\tilde{E}+H{a}^2\tilde{v}\right)={\left(is\right)}^{\beta -1}D{e}^{\beta -1}\frac{{\mathrm{d}}^2\tilde{v}}{\mathrm{d}{\overline{z}}^2},\hfill \end{array}\hfill \end{array}\right.$$

(17)

The simplified form of Equation (17) is

$$\left\{\begin{array}{c}\begin{array}{c}\frac{{\mathrm{d}}^2\tilde{u}}{\mathrm{d}{\overline{z}}^2}-A\tilde{u}=-BHaG,\hfill \\ \frac{{\mathrm{d}}^2\tilde{u}}{\mathrm{d}{\overline{z}}^2}-A\tilde{v}=-B{K}^2\frac{cosh\left(K\overline{z}\right)}{cosh\left(K\right)},\hfill \end{array}\hfill \end{array}\right.$$

(18)

with the boundary conditions

$$\left\{\begin{array}{c}\begin{array}{c}A=\frac{1+{\left(isDe\right)}^{\alpha }}{{\left(isDe\right)}^{\beta -1}}\left(is+H{a}^2\right),\hfill \\ B=\frac{1+{\left(isDe\right)}^{\alpha }}{{\left(isDe\right)}^{\beta -1}}\tilde{E}.\hfill \end{array}\hfill \end{array}\right.$$

(19)

By solving Equation (18), the following equations are obtained:

$$\left\{\begin{array}{c}\begin{array}{c}\tilde{u}=\frac{BHaG}{A}\left[1-\frac{cosh\left(\sqrt{A}\overline{z}\right)}{cosh\sqrt{A}+L\sqrt{A}sinh\sqrt{A}}\right],\hfill \\ \tilde{v}=\frac{B{K}^2}{K^2-A}\left[\frac{\left(1+LKtanhK\right)cosh\left(\sqrt{A}\overline{z}\right)}{cosh\sqrt{A}+L\sqrt{A}sinh\sqrt{A}}-\frac{cosh\left(K\overline{z}\right)}{cosh\left(K\right)}\right].\hfill \end{array}\hfill \end{array}\right.$$

(20)

The velocity can be expressed by inverse of a Fourier transform as follows:

$$\left\{\begin{array}{c}\begin{array}{c}\overline{u}=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\frac{BHaG}{A}\left[1-\frac{cosh\left(\sqrt{A}\overline{z}\right)}{cosh\sqrt{A}+L\sqrt{A}sinh\sqrt{A}}\right]{\mathrm{e}}^{s\overline{t}}\mathrm{ds},\hfill \\ \overline{v}=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\frac{B{K}^2}{K^2-A}\left[\frac{\left(1+LKtanhK\right)cosh\left(\sqrt{A}\overline{z}\right)}{cosh\sqrt{A}+L\sqrt{A}sinh\sqrt{A}}-\frac{cosh\left(K\overline{z}\right)}{cosh\left(K\right)}\right]{\mathrm{e}}^{s\overline{t}}\mathrm{ds}.\hfill \end{array}\hfill \end{array}\right.$$

(21)

For a given external electric field, the velocity can be obtained from Equation (21).

2.4. Numerical Solutions of the EMHD Velocity Field

The numerical solutions can also be provided through the finite difference method. First, a uniform discretization is introduced over the rectangular computational domain $\left(D=\{(\overline{z},\overline{t})\mid \overline{z},\overline{t}\in [-1,1]\times [0,T]\}\right)$ by ${\overline{z}}_j=j\mathsf{\Delta }z$, $j=0,1,2,\cdots ,M$; ${\overline{t}}_n=n\mathsf{\Delta }t$, $n=0,1,2,\cdots ,N$, in which $\mathsf{\Delta }z=2/M$, and $\mathsf{\Delta }t=T/N$ are the spatial and temporal step sizes, respectively. For given $g(\overline{z},\overline{t})$, the following approximations are introduced

$$\frac{\partial g(\overline{z},\overline{t})}{\partial \overline{t}}{\mid }_{\overline{z}={\overline{z}}_j}^{\overline{t}={\overline{t}}_n}\approx \frac{g_{j+1}^n-g_j^n}{\mathsf{\Delta }t},$$

(22)

$$\frac{{\partial }^{2}g(\overline{z},\overline{t})}{\partial {\overline{z}}^2}{\mid }_{\overline{z}={\overline{z}}_j}^{\overline{t}={\overline{t}}_n}\approx \frac{g_{j+1}^n-2g_j^n+g_{j-1}^n}{\mathsf{\Delta }{z}^2}.$$

(23)

Based on the $L_1$ and $L_2$ approximations of the Caputo fractional derivative [27], the following approximations are obtained:

$$D_{\overline{t}}^{\gamma }g(\overline{z},\overline{t}){\mid }_{\overline{z}={\overline{z}}_j}^{\overline{t}={\overline{t}}_n}=\frac{\mathsf{\Delta }{t}^{-\gamma }}{\Gamma (2-\gamma )}\sum _{j=0}^{n-1}{b}_j^{\left(\gamma \right)}\left(g_{{\overline{z}}_j}^{n-j}-g_{{\overline{z}}_j}^{n-j-1}\right),0\le \gamma <1,$$

(24)

$$D_{\overline{t}}^{\gamma }g(\overline{z},\overline{t}){\mid }_{\overline{z}={\overline{z}}_j}^{\overline{t}={\overline{t}}_n}=\frac{\mathsf{\Delta }{t}^{-\gamma }}{\Gamma (3-\gamma )}\sum _{j=0}^{n-1}{b}_j^{(\gamma -1)}\left(g_{{\overline{z}}_j}^{n-j+1}-2g_{{\overline{z}}_j}^{n-j}+g_{{\overline{z}}_j}^{n-j-1}\right),1\le \gamma <2,$$

(25)

$$D_{\overline{t}}^{-\gamma }g(\overline{z},\overline{t}){\mid }_{\overline{z}={\overline{z}}_j}^{\overline{t}={\overline{t}}_n}=\frac{\mathsf{\Delta }{t}^{-\gamma }}{2\Gamma (1-\gamma )}\sum _{j=0}^{n-1}{b}_j^{(1+\gamma )}\left(g_{{\overline{z}}_j}^{n-j}+g_{{\overline{z}}_j}^{n-j-1}\right),\gamma >0,$$

(26)

where $b_j^{\left(\gamma \right)}={(j+1)}^{1-\gamma }-j^{1-\gamma }$$\left(b_0^{\left(\gamma \right)}=1\right)$.

The discretized finite difference scheme of Equation (14) is as follows:

$$\begin{array}{cc}\hfill & \frac{D{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha -1}}{\Gamma (2-\alpha )}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left({\overline{u}}_j^{n-i+1}-2{\overline{u}}_j^{n-i}+{\overline{u}}_j^{n-i-1}\right)+\frac{{\overline{u}}_j^{n+1}-{\overline{u}}_j^n}{\mathsf{\Delta }t}\hfill \\ \hfill & +\frac{H{a}^{2}D{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha }}{\Gamma (2-\alpha )}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left({\overline{u}}_j^{n-i+1}-{\overline{u}}_j^{n-i}\right)+H{a}^2{\overline{u}}_j^n\hfill \\ \hfill & =\frac{D{e}^{\beta -1}\mathsf{\Delta }{t}^{1-\beta }}{\Gamma (2-\beta )}\sum _{i=0}^{n-1}{b}_i^{\left(\beta \right)}\frac{{\overline{u}}_{j+1}^{n-i}-2{\overline{u}}_j^{n-i}+{\overline{u}}_{j-1}^{n-i}}{\mathsf{\Delta }{z}^2}\hfill \\ \hfill & +\frac{HaGD{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha }}{\Gamma (2-\alpha )}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left[\overline{E}\left({\overline{t}}_{n-i}\right)-\overline{E}\left({\overline{t}}_{n-i-1}\right)\right]+HaG\overline{E}\left({\overline{t}}_n\right),\hfill \\ \hfill & n=0,1,2,\cdots ,N-1,j=1,2,3,\cdots ,M-1,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & \frac{D{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha -1}}{\Gamma (2-\alpha )}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left({\overline{v}}_j^{n-i+1}-2{\overline{v}}_j^{n-i}+{\overline{v}}_j^{n-i-1}\right)+\frac{{\overline{v}}_j^{n+1}-{\overline{v}}_j^n}{\mathsf{\Delta }t}\hfill \\ \hfill & +\frac{H{a}^{2}D{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha }}{\Gamma (2-\alpha )}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left({\overline{v}}_j^{n-i+1}-{\overline{v}}_j^{n-i}\right)+H{a}^2{\overline{v}}_j^n\hfill \\ \hfill & =\frac{D{e}^{\beta -1}\mathsf{\Delta }{t}^{1-\beta }}{\Gamma (2-\beta )}\sum _{i=0}^{n-1}{b}_i^{\left(\beta \right)}\frac{{\overline{v}}_{j+1}^{n-i}-2{\overline{v}}_j^{n-i}+{\overline{v}}_{j-1}^{n-i}}{\mathsf{\Delta }{z}^2}\hfill \\ \hfill & +\frac{K^{2}D{e}^{\alpha }\mathsf{\Delta }{t}^{-\alpha }}{\Gamma (2-\alpha )}\frac{cosh\left(K{\overline{z}}_j\right)}{coshK}\sum _{i=0}^{n-1}{b}_i^{\left(\alpha \right)}\left[\overline{E}\left({\overline{t}}_{n-i}\right)-\overline{E}\left({\overline{t}}_{n-i-1}\right)\right]+K^2\frac{cosh\left(K{\overline{z}}_j\right)}{coshK}\overline{E}\left({\overline{t}}_n\right),\hfill \\ \hfill & n=0,1,2,\cdots ,N-1,j=1,2,3,\cdots ,M-1.\hfill \end{array}$$

(28)

Here, Equations (27) and (28) are the discrete schemes of the x-axis and y-axis directions, respectively.

The discretized initial conditions are

$$\left\{\begin{array}{c}\begin{array}{c}\hfill {\overline{u}}_j^0=\overline{u}\left(0\right),\frac{{\overline{u}}_j^1-{\overline{u}}_j^0}{\mathsf{\Delta }t}={\overline{u}}^{\prime }\left(0\right),j=0,1,2,\cdots ,M,\\ \hfill {\overline{v}}_j^0=\overline{v}\left(0\right),\frac{{\overline{v}}_j^1-{\overline{v}}_j^0}{\mathsf{\Delta }t}={\overline{v}}^{\prime }\left(0\right),j=0,1,2,\cdots ,M,\end{array}\hfill \end{array}\right.$$

(29)

where $\overline{u}\left(0\right)$ and $\overline{v}\left(0\right)$ are the values of Equation (21) at $\overline{t}=0$, ${\overline{u}}^{\prime }\left(0\right)$ and ${\overline{v}}^{\prime }\left(0\right)$ are the values at $\overline{t}=0$ after the derivative of Equation (21) with respect to $\overline{t}$. The discretized boundary conditions are

$$\left\{\begin{array}{c}\begin{array}{c}\hfill {\overline{u}}_0^n-L\frac{{\overline{u}}_1^n-{\overline{u}}_0^n}{\mathsf{\Delta }z}=0,{\overline{u}}_M^n+L\frac{{\overline{u}}_{M-1}^n-{\overline{u}}_M^n}{\mathsf{\Delta }z}=0,n=0,1,2,\cdots ,N,\\ \hfill {\overline{v}}_0^n-L\frac{{\overline{v}}_1^n-{\overline{v}}_0^n}{\mathsf{\Delta }z}=0,{\overline{v}}_M^n+L\frac{{\overline{v}}_{M-1}^n-{\overline{v}}_M^n}{\mathsf{\Delta }z}=0,n=0,1,2,\cdots ,N.\end{array}\hfill \end{array}\right.$$

(30)

3. Results and Discussion

In the previous section, through the method of Fourier transform, the analytical expression of the velocity was derived for the EMHD flow of fractional Maxwell fluids between micro-parallel plates. There are three dimensionless numbers involved in Equation (13), namely, Hartmann number $Ha$, dimensionless parameter G and Deborah number $De$. The Hartmann number gives an estimate of the magnetic forces compared to the viscous forces, and G is a non-dimensional parameter representing the strength of the x direction electric field [44]. The Deborah number $De$ is a dimensionless quantity in rheology, which is used to describe the fluidity of materials under certain conditions. It can be used as one of the parameters to measure the viscoelasticity of the fluid. It is defined as the ratio between the relaxation time and the observation time of the mechanical response of the material under the observation conditions. A definition of the Deborah number is given in the research of Ding [19,49], and this definition is also adopted in this paper to study the resonance behavior in the flow process.

In this section, an applied periodic electric field $E=E_{0}cos\left(\omega t\right)$ with the dimensionless frequency $\overline{\omega }=\rho H^2\omega /\mu$ is considered. The ‘ifourier’ command in Matlab symbolic calculation is used to obtain the analytical solution of the velocity. The comparison between the numerical and the analytical solutions is shown in the Figure 2, which shows that they have a good agreement. The fractional derivative parameters $\alpha$ and $\beta$ also play a significant role in the flow, where $0\le \alpha \le \beta \le 1$. Friedrich [50] proved that the rheological constitutive equation exhibits fluid-like behavior only in the case where $\beta =1$, so the fractional Maxwell fluid model with single fractional parameter $\alpha$ is considered. Based on Ref. [44], the electrical conductivity is $\sigma$∼$2.2\times {10}^{-4}$–4$\times {10}^3\mathrm{S}/\mathrm{m}$, viscosity is $\mu$∼${10}^{-3}$–1.5$\times {10}^{-3}\mathrm{kg}/\left(\mathrm{ms}\right)$, the imposed electrical field is $E_0$∼0–30 $\mathrm{V}/\mathrm{m}$, the imposed magnetic field is B∼0.01–5 $\mathrm{T}$, and the Hartmann number is from 0 to 10. Thus, the order of G can be evaluated, and its range is changed from 0 to $6\times {10}^3$; here, it is fixed as 20. This means that the physical configuration can be uniquely determined given the values of these parameters that satisfy the above conditions. For example, when $U_{eo}$∼$2\times {10}^{-4}$, $\sigma$∼$4\times {10}^2$, $\mu$∼$1\times {10}^2$, $E_0$∼$20\mathrm{V}/\mathrm{m}$, h∼$1\times {10}^{-4}$, it corresponds to a physical study of the influence of magnetic field strength and fluid elasticity on the flow.

3.1. Velocity Distribution

Velocity profiles for different fractional parameters $\alpha$ at Hartmann number $Ha=0.5$ are shown in Figure 3. From Figure 3, it can be observed that the amplitude of velocity increases as the $\alpha$ gradually increases. An increased $\alpha$ increases the slip velocity at the wall. From Equation (15), the slip velocity is the product of the fluid shear rate at the wall and the slip length. Therefore, increasing the value of $\alpha$ will increase the fluid shear rate at the wall since the slip length is fixed, thereby promoting slip velocity. Figure 4 exhibits the velocity profiles for different dimensionless slip lengths at dimensionless frequency $\overline{\omega }=0$. It can be seen that as the dimensionless slip length increases, the slip velocity at the wall increases. In addition, when L increases, we can see that the velocity amplitude also has the same elevated trend as the slip velocity. From Figure 4b, the concave part in the middle is caused by the Lorentz force, opposite to the electro-osmotic force.

The variations of velocity profiles with half channel width at different values of the electrokinetic width K for $\alpha =0.7$ and $\beta =1$ are plotted in Figure 5. There is no change in the velocity in the x-axis, which means that the flow perpendicular to the direction of the electric and magnetic fields is induced by the Lorentz force and has nothing to do with the electric double layer; this can be found from Equation (9). In the direction of flow induced by electro-osmotic force, both the slip velocity and amplitude of velocity will increase with the increase in K. Figure 6 shows the effects of the dimensionless Deborah number $De$ on the velocity profiles of EMHD flow for fractional Maxwell fluids. From the figure, an increase in $De$ makes the slip velocity and amplitude of velocity increase. In a physical sense, the relaxation time (its magnitude is characterized by $De$) refers to the time required for a viscoelastic fluid to return to its normal state after deformation. Its magnitude can reflect the elasticity of the fluid. The elastic effect of the fluid can enhance the flow, which has been verified by experiments [51].

Figure 7 illustrates the relationship between velocity profiles and Hartmann number $Ha$ for $\alpha =0.7$ and $\beta =1$. In order to simplify the discussion, the steady case ($\overline{\omega }=0$) is considered, and for this case, $u>0$, $v>0$. In the direction perpendicular to the electric and magnetic fields, the critical value $H{a}_c\simeq 2$ can be found in the figure. When $Ha<H{a}_c$, the wall slip velocity and velocity amplitude increase with the increase in $Ha$, and when $Ha>H{a}_c$, they decrease with the increase in $Ha$. Along the direction of the electric field, the slip velocity and the amplitude of the velocity will decrease with increasing $Ha$. Moreover, at the center of the microchannel, it can be found that a larger value of $Ha$ will have a strong inhibitory effect on the velocity. As $Ha$ increases, the velocity at the center gradually tends to 0. This result can be explained by the force on the fluid in the microchannel. From Equation (9), the body force of the fluid in the x-axis direction and the y-axis direction are $\sigma B\left(E-Bu\right)$ and ${\rho }_{e}E-\sigma B^{2}v$, respectively, and the dimensionless forms are $HaG-H{a}^2\overline{u}$ and $K^{2}cosh\left(K\overline{z}\right)/cosh\left(K\right)-H{a}^2\overline{v}$. In the x-axis direction, the body force on the fluid first increases with the increase in $Ha$, reaches the maximum value at $G/\left(2\overline{u}\right)\simeq H{a}_c$, and then decreases with the increase in $Ha$. Along the direction of the electric field, the body force on the fluid always decreases as $Ha$ increases. Increasing the external force acting on the fluid will increase the velocity amplitude and slip velocity at the wall.

3.2. Resonance Behavior of Fractional Viscoelastic Fluids

In this part, we discuss the flow behavior of fractional viscoelastic fluids through the volumetric flow rate of the fluid. The dimensionless volumetric flow rate is defined as

$$\overline{Q}={\int }_{-1}^1\mathbf{U}(\overline{z},\overline{t})\mathrm{d}\overline{z}.$$

(31)

The variations of the volumetric flow rate for fractional Maxwell model are given graphically at four parameters (fractional parameter $\alpha$, dimensionless slip length L, Deborah number $De$, and Hartmann number $Ha$) in detail. It can be seen from these figures that drastic enhancements of the volumetric flow rate occur at certain frequencies, which reflect resonance phenomena. When resonance occurs, it is obvious that the amplitude (peak) of the enhancement decreases as the frequency increases, and the largest enhancement occurs at the smallest resonance frequency. In addition, in the x-axis direction, the resonance peak at the minimum frequency is large, and then the peaks drops sharply. Therefore, in the x-axis direction, the resonance decay is fast.

Figure 8 illustrates the effects of fractional parameters $\alpha$ on the volumetric flow rate of the x-axis and y-axis, respectively. From the figure, it can be observed that as the $\alpha$ increases, the resonance peak increases. When $\alpha$ decreases, the number of resonance peaks will decrease. This indicates that the resonance behavior of the fractional viscoelastic fluid is suppressed compared with the classical Maxwell fluid. Figure 9 and Figure 10 depicts the volumetric flow rate of the x-axis and y-axis with $\overline{\omega }$. It is noticed that the increase in slip length L or Deborah number $De$ can not only enhance resonance, but also reduce the frequency at which resonance occurs.

Figure 11 shows the effects of Hartmann number $Ha$ on volumetric flow rate of the x-axis and y-axis, respectively. From Figure 11a, a very small Hartmann number will result in a very small volumetric flow rate. When the Hartmann number increases, the volumetric flow rate and resonance peak will increase first, and then decrease. The volumetric flow rate reaches the maximum when $Ha\simeq 1$. The volumetric flow rate represents the net flow of the fluid in the flow process. A slight change in some positions will not affect the value of $Ha$, which can maximize the volumetric flow rate. However, in Section 3.1, the critical number $H{a}_c$ is the value that makes all positions reach the maximum velocity. Therefore, $Ha$ that makes the volumetric flow rate reach the maximum is different from the critical number $H{a}_c$. From Figure 11b, the peaks of the resonance will decrease as $Ha$ increases. It should be mentioned that a larger $Ha$ always suppresses resonance. This means that when a suitable external magnetic field is applied, the resonance behavior will be eliminated.

4. Conclusions

The analytical solutions of combined unsteady electro-osmotic and magnetohydrodynamic velocity are obtained through using the Fourier transform method for the fractional Maxwell model in a slit parallel plate microchannel. Similarly, a numerical solution can also be provided by the finite difference method. The velocity and resonance behavior of periodic EMHD flow are shown graphically and discussed in detail. The resonance greatly enhances the relevant quantities, and the appearance of resonance behavior depends on the value of $\alpha ,L,De,Ha$. The main conclusions are as follows: (i) The fractional derivative parameter $\alpha$, the dimensionless slip length L, and the Deborah number $De$ promote the velocity amplitude and slip velocity, and enhance the resonance behavior. (ii) The electrokinetic width K promotes the flow along the direction of the electric field, and the change of K does not affect the flow in the direction perpendicular to the electric and magnetic fields. (iii) In the direction perpendicular to the electric and magnetic fields, when $Ha$ increases, the flow will be promoted first and then suppressed; in the direction of the electric field, $Ha$ has an inhibitory effect on the flow.