Transient Analysis of the Electro-Osmotic Flow of Multilayer Immiscible Maxwell Fluids in an Annular Microchannel

Abstract

This work investigates the transient multilayer electro-osmotic flow of viscoelastic fluids through an annular microchannel. The dimensionless mathematical model of multilayer flow is integrated by the linearized Poisson-Boltzmann equation, the Cauchy momentum equation, the rheological Maxwell model, initial conditions, and the electrostatic and hydrodynamic boundary conditions at liquid-liquid and solid-liquid interfaces. Although the main force that drives the movement of fluids is due to electrokinetic effects, a pressure gradient can also be added to the flow. The semi-analytical solution for the electric potential distribution and velocity profiles considers analytical techniques as the Laplace transform method, with numerical procedures using the inverse matrix method for linear algebraic equations and the concentrated matrix exponential method for the inversion of the Laplace transform. The results presented for velocity profiles and velocity tracking at the transient regime reveal an interesting oscillatory behavior that depends on elastic fluid properties via relaxation times. The time required for the flow to reach steady-state is highly dependent on the viscosity ratios and the dimensionless relaxation times. In addition, the influence of other dimensionless parameters on the flow as the electrokinetic parameters, zeta potentials at the walls, permittivity ratios, ratio of pressure forces to electro-osmotic forces, number of fluid layers, and annular thickness are investigated. The findings of this study have significant implications for the precise control of parallel fluid transport in microfluidic devices for flow-focusing applications.

1. Introduction

Microfluidics is the science and technology of systems that process or manipulate small (10${}^{-9}$ to 10${}^{-18}$ liters) amounts of fluids, using channels with dimensions of tens to hundreds of micrometers [1]; however, these channels in microfluidics can range from 1 $\mathsf{\mu }$m to 10 $\mathsf{\mu }$m called transitional microchannels, and from 10 $\mathsf{\mu }$m to 200 $\mathsf{\mu }$m called microchannels [2]. Microfluidics evolved from the convergence of technologies and principles from several pre-existing domains, such as chemistry, physics, biology, material science, fluid dynamics, and microelectronics [3]. Therefore, microfluidic devices emerged to overcome the challenges of many applications that required various assays that were not possible or very complicated to perform on the macroscale. As a result, today, great potential has been developed in personalized medicine, disease diagnosis, chemical screening, cell culture, cell separation, cell treatment, drug screening, drug delivery, and DNA sequencing [3]. For example, the lab-on-a-chip (LOC) can concentrate biological and chemical experiments in traditional laboratories on a chip of several square centimeters with high integration [4]. Hence, microfluidic devices have several benefits over conventionally sized systems because they exploit the physicochemical properties of fluids on the microscale. Some benefits of microfluidic devices are: (i) by using small volumes of samples, chemicals, and reagents, reduce overall fees in many applications, (ii) compact size, which results in a reduction of experimental times, (iii) excellent control of parameters and quality of data obtained, which allows the automation of processes, (iv) enhanced the analytical sensitivity and temperature control of samples, (v) more straightforward automation and parallelization, (vi) integration of different lab tasks or routines in one device, and (vii) cheap since they do not necessarily require the integration of various expensive equipment.

Liquid pumping in microfluidic devices is often achieved by electrokinetic effects, where the fluid movement depends on the formation of an electrical double layer on the microchannel walls [5]. Electrokinetics is a general term associated with the relative motion between two charged phases [6]. Electro-osmosis and electrophoresis are the two most widely studied electrokinetic effects, dating back to the early 1800s [7]. The electro-osmotic phenomenon was first reported as a transport technique in 1809 by Reuss [8], who demonstrated that water could flow through a plug of clay by applying an electrical voltage. Later in 1816, Porrett [9] independently reported the discovery of the electro-osmotic phenomenon; although his experiments were similar to those of Reuss, he used a bladder membrane as the porous barrier. From the above and throughout the years from 1879 to 1924, was developed the fundamental theory of the electric double layer, a basic concept for understanding an electro-osmotic flow [10,11,12,13]. With earlier knowledge of electrokinetic phenomena, the method of pumping electrolytes in narrow capillaries was first investigated by Burgreen and Nakache [14], Rice and Whitehead [15], and Levine et al. [16] in the years 1964, 1965, and 1974, respectively.

With this background, in the past three decades until today, electro-osmotic flow as a pumping mechanism has become an essential tool for fluid manipulation in microsystems. The relevance of electro-osmotic micropumps is because they have a simple and compact structure, unlike pressure-driven conventional mechanical micropumps, whose manufacturing is based on moving parts [17,18]. Therefore, considering the impact of electro-osmotic flow in miniaturized systems, analytical, numerical, and experimental investigations have been carried out on the transport of Newtonian fluids. Andreev and Lisin [19] solved numerically the electro-osmotic flow effect on the efficiency of capillary electrophoresis in a wide range of buffer concentrations and capillary diameters. Osuga et al. [20] analyze in detail the growth process and stationary flow profile of electro-osmotic flow for a large capillary-radius-to-Debye length ratio. Patankar and Hu [21] developed a numerical scheme to simulate electro-osmotic flows at the intersection of a cross-channel during the injection process for capillary electrophoresis. The numerical simulations help to determine the important parameters in controlling the shape of the injected sample. Tsao [22] studied the electro-osmotic flow in an annular channel, considering that the Helmholtz-Smoluchowski equation can describe the electro-osmotic mobility; however, this equation must be corrected by a geometry-dependent factor to evaluate the mean velocity correctly. Chang et al. [23] analytically solved the start-up electro-osmotic flow through a small pore or microchannel with an annular or rectangular cross-section under the Debye-Hückel approximation. While Yavari et al. [24] analyzed the combined electro-osmotic and pressure-driven flow in a microannulus involving the uniform Joule heating and negligible viscous dissipation. Additionally, throughout time several works about electro-osmotic flows have been developed with different conduit geometries such as parallel flat plate channels [25,26,27], rectangular channels [28,29,30], and cylindrical channels [25,31,32]. Trying to cover the wide range of applications of microfluidics, the electro-osmotic flow has been extended to handling non-Newtonian fluids to predict their complex behavior. Therefore, many studies incorporate the rheological models of Phan-Thien-Tanner [33,34], Oldroyd [35], Maxwell [36,37,38], and power-law [39,40], among others [41,42,43].

Furthermore to considering the flow of single-phase fluids such as those shown above, many applications in microfluidics have benefited from the transport of two or more immiscible fluids in microfiber generation processes [44,45,46], continuous chemical processing [47], transfer of enzyme-based processes [48], separation of antibiotics [49], viscous drag pumping [50] and flow-focusing processes [51], synthesis of nanomaterials [52], reaction [53], microextraction [54], drug analysis in blood samples [55], purification of analytes from biomedical samples [56], isolation of blood components [57], droplet formation [58] and many other applications mentioned by references [58,59,60]. In this context, parallel flows are widely used in different microsystems involving two- or three-phase fluid flows, which require developing stable liquid-liquid interfaces and maintaining a continuous parallel laminar flow pattern in the channels [58,61,62,63,64].

Following the research on parallel flows, viscous drag pumping is a commonly used method for driving a non-conducting fluid through neighboring conducting fluids. This pumping technique can be implemented in biomedical applications and chemical analysis where some fluid samples such as oil, ethanol and organic solvents have low electrical conductivity [65,66]. Gao et al. [67] presented a two-fluid electro-osmotic pump, assuming that a conducting liquid driven by electro-osmotic effects drags a non-conducting liquid by viscous forces. Their proposal was based on the fact that a fluid with low electrical conductivity cannot form an electrical double layer and, therefore, cannot be pumped alone by electro-osmotic effects. Following this idea, Haiwang et al. [68] analyzed a three-fluid electro-osmotic pumping technique, in which an electrically non-conducting fluid is delivered by the viscous interfacial force of two conducting fluids. Other studies on viscous drag pumping using combined electro-osmotic and pressure gradient effects for chemical and biomechanical processes, injection analysis, and microelectronic cooling systems can be reviewed in references [69,70,71]. The parallel flow systems mentioned above do not use special liquid-liquid boundary conditions as one of the fluids is non-conducting; then, to solve the electric potential, a specified zeta potential at liquid-liquid interfaces is enough. However, in another case, considering all the transported fluids are electrically conducting, the interfacial phenomenon include forming an electrical double layer both in solid-liquid [6] and liquid-liquid interfaces [72]. Therefore, the interfacial electrostatic effects necessary to solve the electric potential and velocity distributions should consider a potential difference, Gauss’s law, and the Maxwell electric stresses. Under this electrostatic conditions, different works about the electro-osmotic transport of two conducting fluids with Newtonian [73,74,75] and non-Newtonian [76] behavior were carried out, considering applications involving biological testing processes, biomedical diagnosis, and heat transfer. These investigations report that the discontinuity of the electrical potential and viscous stresses in the liquid-liquid interfaces produce speed jumps in these regions. Furthermore, the polarity of the concentrated electric charges within the electric double layers in the contact of these conducting fluids plays an important role in shaping the velocity profile of these flows.

Advances in microfabrication technologies have allowed the development of microfluidic devices capable of handling several fluid layers in the same conduit [45,51,77,78,79]. For this reason, theoretical research about electro-osmotic flows that handle multilayer parallel fluids has been carried out to simulate their behavior in microsystems. For Newtonian fluids, Torres and Escandón [80] in the year 2020 solved the transient electro-osmotic flow of multilayer immiscible fluids in a narrow cylindrical capillary. In the case of non-Newtonian fluids, the research dates from 2011 with the study of Li et al. [81], who analyzed the steady combined electro-osmotic and pressure-driven flow of multilayer immiscible fluids in a microchannel formed by two parallel flat plates; here, a power-law fluid is imposed in the middle layer. To extend the proposal of Li et al. [81], on the inclusion of complex fluid behavior in the multilayer scheme, recently, Escandón et al. in [82,83] studied multilayer electro-osmotic flows of Phan-Thien-Tanner and Maxwell fluids in a slit microchannel, respectively. These investigations highlighted the importance of rheological effects on the multilayer flow field in both steady and transient states.

This research intends to extend the knowledge about control techniques in multilayer flows under electrokinetic effects. Therefore, this work aims to solve the transient electro-osmotic flow of multilayer immiscible Maxwell fluids in an annular microchannel following the proposed design by [79] because this channel geometry is more realistic within multilayer flows; however, it has not been addressed yet. Also, this parametric study will focus on the interaction between the non-Newtonian behavior and the interfacial electrostatic phenomena of the conducting transported fluids. In particular, this investigation is helpful for applications handling biofluids or polymeric solutions with viscoelastic behavior, where the start-up of the flow requires a prediction with high space-time precision.

2. Mathematical Modeling

2.1. Physical Model Description

The present study theoretically analyzes an electro-osmotic pump where the multilayer flow generator device presented in Figure 1, is an adaptation of the coaxial device described by Morimoto et al. [79]. Here, the annular microchannel is formed of two concentric tubes, and the combination of the microfluidic modules allows to increase the number of fluid layers from $n=1$ to $n=i$. The cylindrical coordinate system $(r,z)$ is shown in the enlarged view with the red square into Figure 1 and marks the origin of the studied section of this electrokinetic pump. The concentric tubes that form the annular cavity have an external radius R and an internal radius $aR$, where $0<a<1$. Each n-fluid layer contains a mixture of a symmetrical electrolyte and a solute having viscoelastic characteristics. In Figure 1, the liquid-liquid interfaces are placed in $r_n$ positions within the range of $aR<r_n<R$; $n=1,2,3,\dots ,i$ are subscripts that indicate the number of the corresponding fluid layer, where $n=1$ and $n=i$ are the fluid layers in contact with the inner and outer wall of the microchannel, respectively. Fluids are immiscible and electrically conductive, and the liquid-liquid interfaces between them can be considered polarizable and impermeable to charged particles. As a result, electrostatic effects appear at the liquid-liquid interfaces caused by the formation of an electric double layer deriving in a potential difference $\Delta {\psi }_n$. In addition, an electrostatic interaction with the surrounding fluids arises on the microchannel walls, forming another double electric layer at the solid-liquid interfaces whose surface charge is represented by the zeta potentials ${\zeta }_{\mathrm{I}}$ and ${\zeta }_{\mathrm{O}}$ for the inner and outer walls, respectively. With the above, the movement of fluids is produced by applying a uniform electric field $E_z$ which induces electro-osmotic effects on the electric double layers. The electric field is produced by electrodes located at the ends of the conduit and the beginning of each green module. Additionally to the electro-osmotic forces, each fluid layer can be subjected to an external pressure gradient $p_z$.

2.2. Governing and Constitutive Equations

The electro-osmotic flow field for the incompressible n-fluid layers is governed by the following Poisson, continuity and Cauchy momentum equations:

$${\nabla }^2{\mathrm{\Phi }}_n=-\frac{{\rho }_{\mathrm{e},n}}{{ϵ}_n},$$

(1)

$$\nabla ·{\mathbf{v}}_n=0$$

(2)

and

$${\rho }_n\frac{D{\mathbf{v}}_n}{Dt}=-\nabla p+\nabla ·{\mathbf{\tau }}_n+{\rho }_n\mathbf{g}+{\rho }_{\mathrm{e},n}\mathbf{E},$$

(3)

where $\mathrm{\Phi }$ is the electric potential, ${\rho }_{\mathrm{e}}$ is the volumetric free charge density, $ϵ$ is the dielectric permittivity, $\mathbf{v}$ is the velocity vector, $\rho$ is the fluid density, t is the time, p is the pressure, $\mathbf{\tau }$ is the stress tensor, $\mathbf{g}$ is the gravitational acceleration vector and $\mathbf{E}$ is the electric field vector. Moreover, the stress tensor is related to the Maxwell rheological model as [84,85]:

$${\mathbf{\tau }}_n+{\lambda }_n\frac{\partial {\mathbf{\tau }}_n}{\partial t}={\eta }_{0,n}{\dot{\mathbf{\gamma }}}_n,$$

(4)

where $\lambda$ is the relaxation time, ${\eta }_0$ is the zero-shear-rate viscosity and ${\dot{\mathbf{\gamma }}}_n=(\nabla {\mathbf{v}}_n)+{(\nabla {\mathbf{v}}_n)}^T$ is the rate-of-strain tensor. The general governing and constitutive equations mentioned above can be simplified with the following assumptions:

• Transient and fully-developed flow.
• The physical and electrical properties of fluids are considered constant [67].
• Planar interfaces between the fluids [69,73,81] which can be assumed by considering the following. (i) very low Reynolds numbers, i.e., $R{e}_n(={\rho }_{n}R{u}_{\mathrm{HS}}/{\eta }_{0,n})\ll 1$, which results in parallel flows [68,81]. (ii) Uniform zeta potentials along the study section of the annular microchannel. And (iii) the gravity does not affect the flow in the channel [69,81]. The characteristic flow velocity is defined by the Helmholtz-Smoluchowski velocity as $u_{\mathrm{HS}}=-{ϵ}_{\mathrm{ref}}{\zeta }_{\mathrm{ref}}{E}_z/{\eta }_{\mathrm{ref}}$, where the subscript “$\mathrm{ref}$” indicates reference quantities.
• The electrochemical interface structure between the immiscible fluids follows the Verwey-Niessen theory [72,86]. Therefore, the electric double layers at the liquid-liquid interfaces are two diffuse charge layers separated by a central inner compact layer. The latter is characterized by a potential drop between the two diffuse layers due to the orientation of the solvent molecules. Also, the continuity of electrical displacements on both sides of the central compact layer in the absence of ions in the inner layer is considered [86,87].
• The annular microchannel is considered very long, and the study region neglects any end effect. Thus, the electric potential ${\mathrm{\Phi }}_n$ can be determined as the superposition of the potential into the electric double layer, ${\psi }_n$, with the external potential, $\varphi$, yielding [6]:

  $${\mathrm{\Phi }}_n(r,z)={\psi }_n\left(r\right)+\varphi \left(z\right),$$

  (5)

  where $\varphi \left(z\right)={\varphi }_0-z{E}_z$, ${\varphi }_0$ is the electric potential at the inlet of the microchannel at $z=0$, and $E_z$ is the constant electric field in the z-axis.
• The free charge density in terms of the electrical potential, ${\psi }_n$, can be described by the Boltzmann distribution as [6]:

  $${\rho }_{\mathrm{e},n}=-2Z_{n}e{n}_{0,n}sinh\left(\frac{Z_{n}e{\psi }_n\left(r\right)}{k_{\mathrm{B}}{T}_n}\right),$$

  (6)

  where $Z_n$ is the valence of electrolyte, e is the electron charge, $n_{0,n}$ is the ionic number concentration in the bulk solution, $k_{\mathrm{B}}$ is the Boltzmann constant and $T_n$ is the fluid temperature.
• Low interfacial potentials to consider the Debye-Hückel approximation at solid-liquid [6,88] and liquid-liquid interfaces [74,76]. This assumption is valid until 50 mV [6,15,25].
• The electric double layers do not overlap.
• The pressure gradient $p_z$ applied in the z-axis of the microchannel remains constant and is produced by syringe pumps.

Once the assumptions above are considered, the result is a unidirectional flow with flat liquid-liquid interfaces between the fluids. Then, (1)–(4) are rewritten as follows, the Poisson-Boltzmann equation as:

$$\frac{d^2{\psi }_n\left(r\right)}{d{r}^2}+\frac{1}{r}\frac{d{\psi }_n\left(r\right)}{dr}={\kappa }_n^2{\psi }_n\left(r\right),$$

(7)

while the momentum equation reduces to:

$${\rho }_n\frac{\partial u_n(r,t)}{\partial t}=-p_z-\frac{1}{r}{\tau }_{rz,n}(r,t)-\frac{\partial {\tau }_{rz,n}(r,t)}{\partial r}-{ϵ}_n{\kappa }_n^2{E}_z{\psi }_n\left(r\right)$$

(8)

and the constitutive equation for Maxwell fluids can be expressed as:

$${\tau }_{rz,n}(r,t)+{\lambda }_n\frac{\partial {\tau }_{rz,n}(r,t)}{\partial t}=-{\eta }_{0,n}\frac{\partial u_n(r,t)}{dr},$$

(9)

where ${\kappa }_n^2=2Z_n^2{e}^2{n}_{0,n}/{ϵ}_n{k}_B{T}_n$ is the Debye-Hückel parameter related to the Debye length ${\kappa }_n^{-1}={\left({ϵ}_n{k}_{\mathrm{B}}{T}_n/2Z_n^2{e}^2{n}_{0,n}\right)}^{1/2}$ [6], $u_n$ is the fluid velocity on the $z–$axis and ${\tau }_{rz,n}$ is the shear stress.

To obtain the momentum Equation (8) in terms of velocity and be more easily solved, (9) has to be derived with respect to the r-axis, resulting in:

$$\frac{\partial {\tau }_{rz,n}(r,t)}{\partial r}=-{\lambda }_n\frac{{\partial }^2{\tau }_{rz,n}(r,t)}{\partial t\partial r}-{\eta }_{0,n}\frac{{\partial }^2{u}_n(r,t)}{\partial r^2},$$

(10)

after, substituting (10) into (8), yields:

$${\rho }_n\frac{\partial u_n(r,t)}{\partial t}=-p_z+{\lambda }_n\frac{{\partial }^2{\tau }_{rz,n}(r,t)}{\partial t\partial r}+{\eta }_{0,n}\frac{{\partial }^2{u}_n(r,t)}{\partial r^2}-\frac{1}{r}{\tau }_{rz,n}(r,t)-{ϵ}_n{\kappa }_n^2{E}_z{\psi }_n\left(r\right).$$

(11)

In this direction, (8) is derived with respect to time to have:

$$\frac{{\partial }^2{\tau }_{rz,n}(r,t)}{\partial t\partial r}=-{\rho }_n\frac{{\partial }^2{u}_n(r,t)}{\partial t^2}-\frac{1}{r}\frac{\partial {\tau }_{rz,n}(r,t)}{\partial t}.$$

(12)

Therefore, by replacing (12) into (11) the following relationship is obtained:

$${\rho }_n{\lambda }_n\frac{{\partial }^2{u}_n(r,t)}{\partial t^2}+{\rho }_n\frac{\partial u_n(r,t)}{\partial t}=-p_z-\frac{1}{r}{\tau }_{rz,n}(r,t)-\frac{{\lambda }_n}{r}\frac{\partial {\tau }_{rz,n}(r,t)}{\partial t}+{\eta }_{0,n}\frac{{\partial }^2{u}_n(r,t)}{\partial r^2}-{ϵ}_n{\kappa }_n^2{E}_z{\psi }_n\left(r\right).$$

(13)

Finally, to eliminate ${\tau }_{rz,n}$ from (13), (9) is divided by the radial coordinate r and replaced into (13), resulting in the following momentum equation for $n–$fluid layers as a function of velocity only:

$${\rho }_n{\lambda }_n\frac{{\partial }^2{u}_n(r,t)}{\partial t^2}+{\rho }_n\frac{\partial u_n(r,t)}{\partial t}=-p_z+\frac{{\eta }_{0,n}}{r}\frac{\partial u_n(r,t)}{dr}+{\eta }_{0,n}\frac{{\partial }^2{u}_n(r,t)}{\partial r^2}-{ϵ}_n{\kappa }_n^2{E}_z{\psi }_n\left(r\right).$$

(14)

The governing equations given by (7) and (14), referring to the electric potential and velocity, are solved by applying the corresponding boundary conditions at $t>0$ as follows. First, for fluid layer 1 in contact with the microchannel inner wall, the boundary conditions that are taken into account and evaluated at $r=aR$ are a specified zeta potential and the no-slip boundary conditions as:

$${\psi }_1(r=aR)={\zeta }_{\mathrm{I}},u_1(r=aR,t)=0.$$

(15)

Second, the boundary conditions at the liquid-liquid interfaces have to be evaluated at $r=r_{n=1,2,3,\dots ,i-1}$. For the electric potential are a potential difference and the Gauss’s law as:

$${\psi }_{n+1}\left(r\right)-{\psi }_n\left(r\right)=\Delta {\psi }_n,$$

(16)

$${ϵ}_{n+1}\frac{d{\psi }_{n+1}\left(r\right)}{dr}={ϵ}_n\frac{d{\psi }_n\left(r\right)}{dr},$$

(17)

and for the velocity, the continuity of velocity:

$$u_{n+1}(r,t)=u_n(r,t)$$

(18)

and the electroviscous stress balance:

$${\tau }_{rz,n+1}(r,t)+{\tau }_{e,n+1}\left(r\right)={\tau }_{rz,n}(r,t)+{\tau }_{e,n}\left(r\right),$$

(19)

where the Maxwell shear stress is:

$${\tau }_{e,n}\left(r\right)=-{ϵ}_n{E}_z\frac{d{\psi }_n\left(r\right)}{dr}.$$

(20)

And third, for fluid layer i in contact with the microchannel outer wall, the boundary conditions that are taken into account and evaluated at $r=R$ are a specified zeta potential and the no-slip boundary conditions as:

$${\psi }_i(r=R)={\zeta }_{\mathrm{O}},u_i(r=R,t)=0.$$

(21)

Complementary, the initial conditions applicable to the momentum Equation (14) for all fluids and in the entire geometric domain of the annular microchannel ($aR\le r\le R$) are as follows:

$$u_n(r,t=0)=0,{\tau }_{rz,n}(r,t=0)=0,{\left.\frac{\partial u_n}{\partial t}\right|}_{r,t=0}=0.$$

(22)

2.3. Dimensionless Mathematical Model

To normalize the mathematical model defined above, are introduced the following dimensionless variables:

$$\overline{t}=\frac{{\eta }_{\mathrm{ref}}t}{{\rho }_{\mathrm{ref}}{R}^2},\overline{r}=\frac{r}{R},{\overline{\psi }}_n=\frac{{\psi }_n}{{\zeta }_{\mathrm{ref}}},{\overline{u}}_n=\frac{u_n}{u_{\mathrm{HS}}},{\overline{\tau }}_{rz,n}=\frac{R{\tau }_{rz,n}}{{\eta }_{\mathrm{ref}}{u}_{\mathrm{HS}}}.$$

(23)

Therefore, by replacing (23) into (7), (9) and (14), the Poisson-Boltzmann, momentum and Maxwell equations are obtained in their dimensionless form as follows:

$$\frac{d^2{\overline{\psi }}_n\left(\overline{r}\right)}{d{\overline{r}}^2}+\frac{1}{\overline{r}}\frac{d{\overline{\psi }}_n\left(\overline{r}\right)}{d\overline{r}}={\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right),$$

(24)

$${\overline{\rho }}_n{\overline{\lambda }}_n\frac{{\partial }^2{\overline{u}}_n(\overline{r},\overline{t})}{\partial {\overline{t}}^2}+{\overline{\rho }}_n\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{t}}=-\Gamma +{\overline{\eta }}_n\frac{{\partial }^2{\overline{u}}_n(\overline{r},\overline{t})}{\partial {\overline{r}}^2}+\frac{{\overline{\eta }}_n}{\overline{r}}\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{r}}+{\overline{ϵ}}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right)$$

(25)

and

$${\overline{\tau }}_{rz,n}(\overline{r},\overline{t})+{\overline{\lambda }}_n\frac{\partial {\overline{\tau }}_{rz,n}(\overline{r},\overline{t})}{\partial \overline{t}}=-{\overline{\eta }}_n\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{r}}.$$

(26)

As well, all boundary conditions are normalized by use the dimensionless variables given in (23) for $\overline{t}>0$ as follows. Equation (15) for solid-liquid interface placed at the microchannel inner wall at $\overline{r}=a$, can be rewritten as:

$${\overline{\psi }}_1(\overline{r}=a)={\overline{\zeta }}_{\mathrm{I}},{\overline{u}}_1(\overline{r}=a,\overline{t})=0.$$

(27)

Respectively, from (16)–(20) for all liquid-liquid interfaces placed at $\overline{r}={\overline{r}}_{n=1,2,3,\dots ,i-1}$, the following dimensionless equations are obtained:

$${\overline{\psi }}_{n+1}\left(\overline{r}\right)-{\overline{\psi }}_n\left(\overline{r}\right)=\Delta {\overline{\psi }}_n,$$

(28)

$${\overline{ϵ}}_{n+1}\frac{d{\overline{\psi }}_{n+1}\left(\overline{r}\right)}{d\overline{r}}={\overline{ϵ}}_n\frac{d{\overline{\psi }}_n\left(\overline{r}\right)}{d\overline{r}},$$

(29)

$${\overline{u}}_{n+1}(\overline{r},\overline{t})={\overline{u}}_n(\overline{r},\overline{t})$$

(30)

and

$${\overline{\tau }}_{rz,n+1}(\overline{r},\overline{t})+{\overline{ϵ}}_{n+1}\frac{d{\overline{\psi }}_{n+1}\left(\overline{r}\right)}{d\overline{r}}={\overline{\tau }}_{rz,n}(\overline{r},\overline{t})+{\overline{ϵ}}_n\frac{d{\overline{\psi }}_n\left(\overline{r}\right)}{d\overline{r}}.$$

(31)

Also, (21) for the solid-liquid interface placed at the microchannel outer wall at $\overline{r}=1$, can be rewritten as:

$${\overline{\psi }}_i(\overline{r}=1)={\overline{\zeta }}_{\mathrm{O}},{\overline{u}}_i(\overline{r}=1,\overline{t})=0.$$

(32)

Regarding the initial conditions in (22), these also are normalized by using (23), yielding:

$${\overline{u}}_n(\overline{r},\overline{t}=0)=0,{\overline{\tau }}_{rz,n}(\overline{r},\overline{t}=0)=0,{\left.\frac{\partial {\overline{u}}_n}{\partial \overline{t}}\right|}_{\overline{r},\overline{t}=0}=0,$$

(33)

for $a\le \overline{r}\le 1$.

The dimensionless parameters arising in this section are described mathematically as:

$${\overline{\kappa }}_n=\frac{R}{{\kappa }^{-1}},{\overline{\rho }}_n=\frac{{\rho }_n}{{\rho }_{\mathrm{ref}}},{\overline{ϵ}}_n=\frac{{ϵ}_n}{{ϵ}_{\mathrm{ref}}},{\overline{\eta }}_n=\frac{{\eta }_{0,n}}{{\eta }_{\mathrm{ref}}},\Gamma =\frac{p_z{R}^2}{U_{\mathrm{HS}}{\eta }_{\mathrm{ref}}},{\overline{\lambda }}_n=\frac{{\eta }_{\mathrm{ref}}{\lambda }_n}{{\rho }_{\mathrm{ref}}{R}^2},{\overline{r}}_n=\frac{r_n-a}{R}+a,\Delta {\overline{\psi }}_n=\frac{\Delta {\psi }_n}{{\zeta }_{\mathrm{ref}}},$$

(34)

where ${\overline{\kappa }}_n$ are electrokinetic parameters representing the ratios between the microchannel height to the Debye lengths, ${\overline{\rho }}_n$ are the densities ratios, ${\overline{ϵ}}_n$ are the dielectric permittivities ratios, ${\overline{\eta }}_n$ are the viscosity ratios, $\Gamma$ is the ratio of external pressure forces to electro-osmotic forces, and ${\overline{\lambda }}_n$ are the dimensionless relaxation times. Besides, ${\overline{r}}_n$ and $\Delta {\overline{\psi }}_n$ are the interface positions and the potential differences, respectively, both dimensionless and ranging from $n=1$ to $n=i-1$.

3. Solution Methodology

3.1. Electric Potential Distribution

Because the Poisson-Boltzmann equation for the electric potential (24) has the form of a modified Bessel differential equation, it results in the following well-known solution for each n-fluid layer:

$${\overline{\psi }}_n\left(\overline{r}\right)=C_{2n-1}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+C_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right),$$

(35)

where $C_{2n-1}$ and $C_{2n}$ are integration constants whose subscripts are according to each fluid layer, $I_0(·)$ and $K_0(·)$ are the modified Bessel functions of the first and second kind of zero-order, respectively. To obtain the integration constants value, the corresponding boundary conditions given in (27)–(29) and (32) were applied with the aid of (35). As result, a system of linear equations represented by a general matrix notation as $\mathbf{A}\mathbf{x}=\mathbf{b}$ is shown in Appendix A. Thus, the matrix inverse method is used in (A1) to find $C_{2n-1}$ and $C_{2n}$.

3.2. Velocity Profiles

The velocity distribution for this multilayer flow is obtained by employing the Laplace transform method as follows:

$$U_n(\overline{r},s)=\mathcal{L}\left\{{\overline{u}}_n(\overline{r},\overline{t})\right\}={\int }_0^{\infty }{\overline{u}}_n(\overline{r},\overline{t})e^{-s\overline{t}}d\overline{t}$$

(36)

and for the shear stress as:

$${\tilde{\tau }}_{rz,n}(\overline{r},s)=\mathcal{L}\left\{{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})\right\}={\int }_0^{\infty }{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})e^{-s\overline{t}}d\overline{t}.$$

(37)

Applying the above pair of equations to the momentum and constitutive equations given in (25) and (26), results in:

$$\begin{array}{c}\hfill {\overline{\lambda }}_n{\overline{\rho }}_n\left[s^2{U}_n(\overline{r},s)-s{\overline{u}}_n(\overline{r},\overline{t}=0)-{\left.\frac{\partial {\overline{u}}_n}{\partial \overline{t}}\right|}_{\overline{r},\overline{t}=0}\right]+{\overline{\rho }}_n\left[s{U}_n(\overline{r},s)-{\overline{u}}_n(\overline{r},\overline{t}=0)\right]=\\ \hfill {\overline{\eta }}_n\frac{{\partial }^2{U}_n(\overline{r},s)}{\partial {\overline{r}}^2}+\frac{{\overline{\eta }}_n}{\overline{r}}\frac{\partial U_n(\overline{r},s)}{\partial \overline{r}}+\frac{{\overline{ϵ}}_n{\overline{\kappa }}_n^2}{s}{\overline{\psi }}_n\left(\overline{r}\right)-\frac{\Gamma }{s}\end{array}$$

(38)

and respectively

$$\begin{array}{c}\hfill {\tilde{\tau }}_{rz,n}(\overline{r},s)+{\overline{\lambda }}_n\left[s{\tilde{\tau }}_{rz,n}(\overline{r},s)-{\overline{\tau }}_{rz,n}(\overline{r},\overline{t}=0)\right]=-{\overline{\eta }}_n\frac{\partial U_n(\overline{r},s)}{\partial \overline{r}}.\end{array}$$

(39)

The initial conditions given in the (33) satisfy the (38) and (39), and these can be rewritten as:

$$\begin{array}{c}\hfill \frac{{\partial }^2{U}_n(\overline{r},s)}{\partial {\overline{r}}^2}+\frac{1}{\overline{r}}\frac{\partial U_n(\overline{r},s)}{\partial \overline{r}}-{\alpha }_n^2{U}_n(\overline{r},s)={\beta }_n{\overline{\psi }}_n\left(\overline{r}\right)+\frac{\Gamma }{s{\overline{\eta }}_n}\end{array}$$

(40)

and

$${\tilde{\tau }}_{rz,n}(\overline{r},s)=-{\gamma }_n\frac{\partial U_n(\overline{r},s)}{\partial \overline{r}},$$

(41)

where ${\alpha }_n^2=({\overline{\rho }}_{n}s/{\overline{\eta }}_n)({\overline{\lambda }}_{n}s+1)$, ${\beta }_n=-{\overline{ϵ}}_n{\overline{\kappa }}_n^2/{\overline{\eta }}_{n}s$, and ${\gamma }_n={\overline{\eta }}_n/(1+{\overline{\lambda }}_{n}s)$. To solve the momentum Equation (40), requires suitable boundary conditions which arise by apply (36) and (37) in (27) and in (30)–(32). As a result, for the inner wall of the annular microchannel at $\overline{r}=a$:

$$U_1(\overline{r}=a,s)=0,$$

(42)

while for all liquid-liquid interfaces at $\overline{r}={\overline{r}}_{n=1,2,3,\dots i-1}$ and with the aid of (41), are obtained the following expressions:

$$U_{n+1}(\overline{r}={\overline{r}}_n,s)=U_n(\overline{r}={\overline{r}}_n,s),$$

(43)

$$-{\gamma }_{n+1}\frac{\partial U_{n+1}}{\partial \overline{r}}+\frac{{\overline{ϵ}}_{n+1}}{s}\frac{d{\overline{\psi }}_{n+1}}{d\overline{r}}=-{\gamma }_n\frac{\partial U_n}{\partial \overline{r}}+\frac{{\overline{ϵ}}_n}{s}\frac{d{\overline{\psi }}_n}{d\overline{r}}$$

(44)

and finally, for the outer wall of the annular microchannel at $\overline{r}=1$, it results in:

$$U_i(\overline{r}=1,s)=0.$$

(45)

Hence, the mathematical model for this electro-osmotic flow in Laplace space consists of (40) and (42)–(45), which is solved as described below. In general, (40) represents an non-homogeneous ordinary differential equation and its solution can be constructed by the principle of superposition of a homogeneous solution $U_{h,n}(\overline{r},s)$ and a particular solution $U_{p,n}(\overline{r},s)$, as follows:

$$U_n(\overline{r},s)=U_{h,n}(\overline{r},s)+U_{p,n}(\overline{r},s),$$

(46)

where the homogeneous and particular solution are:

$$U_{h,n}(\overline{r},s)=A_n{I}_0\left({\alpha }_n\overline{r}\right)+B_n{K}_0\left({\alpha }_n\overline{r}\right)$$

(47)

and

$$U_{p,n}(\overline{r},s)=F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+H_0.$$

(48)

The particular solution given in (48) consists of two parts. The first part corresponds to $F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)$ and retains the form of a modified Bessel differential equation, as seen on the right-hand side of (35). The second part corresponds to $H_0$, which is treated as a constant. Both parts are also related to the terms on the right-hand side of (40). Therefore, $A_n$, $B_n$, $F_n$, $G_n$, and $H_0$ are constants to be determined. In this sense, substituting (48) in (40), and with the aid of (35), yields:

$$\begin{array}{c}\hfill \frac{{\partial }^2}{\partial {\overline{r}}^2}\left[F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]+\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\left[F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]-{\alpha }_n^2\left[F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+H_0\right]=\\ \hfill {\beta }_n\left[C_{2n-1}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+C_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]+\frac{\Gamma }{s{\overline{\eta }}_n}.\end{array}$$

(49)

Considering that the Bessel functions depend only on the radial coordinate and by factoring $F_n$ and $G_n$, from (49) is obtained that:

$$\begin{array}{c}\hfill F_n\left[\frac{d^2}{d{\overline{r}}^2}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{I}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]+G_n\left[\frac{d^2}{d{\overline{r}}^2}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]-{\alpha }_n^2{H}_0=\\ \hfill {\beta }_n\left[C_{2n-1}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+C_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]+\frac{\Gamma }{s{\overline{\eta }}_n}.\end{array}$$

(50)

Following the principle of superposition of solutions, the global solution will be the sum of all solutions, which can be manipulated independently. Therefore, from (50) the following three relationships are obtained:

$$F_n\left[\frac{d^2}{d{\overline{r}}^2}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{I}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]={\beta }_n{C}_{2n-1}{I}_0\left({\overline{\kappa }}_n\overline{r}\right),$$

(51)

$$G_n\left[\frac{d^2}{d{\overline{r}}^2}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]={\beta }_n{C}_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)$$

(52)

and

$$-{\alpha }_n^2{H}_0=\frac{\Gamma }{s{\overline{\eta }}_n}.$$

(53)

However, by replacing (35) into (24), yields:

$$\frac{d^2}{d{\overline{r}}^2}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)={\overline{\kappa }}_n^2{I}_0\left({\overline{\kappa }}_n\overline{r}\right)$$

(54)

and

$$\frac{d^2}{d{\overline{r}}^2}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+\frac{1}{\overline{r}}\frac{d}{d\overline{r}}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)={\overline{\kappa }}_n^2{K}_0\left({\overline{\kappa }}_n\overline{r}\right).$$

(55)

Therefore, replacing the equalities of (54) and (55) into (51) and (52), respectively, these can be rewritten as:

$$F_n\left[{\overline{\kappa }}_n^2{I}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{I}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]={\beta }_n\left[C_{2n-1}{I}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]$$

(56)

and

$$G_n\left[{\overline{\kappa }}_n^2{K}_0\left({\overline{\kappa }}_n\overline{r}\right)-{\alpha }_n^2{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right]={\beta }_n\left[C_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)\right].$$

(57)

The constants $F_n$ and $G_n$ are obtained from the above pair of equations, yielding:

$$F_n=\frac{{\beta }_n{C}_{2n-1}}{{\overline{\kappa }}_n^2-{\alpha }_n^2},G_n=\frac{{\beta }_n{C}_{2n}}{{\overline{\kappa }}_n^2-{\alpha }_n^2}.$$

(58)

Taking into account the constant values of (53) and (58), the dimensionless velocity distribution for each one $n–$fluid layer is established from (46), yielding:

$$U_n(\overline{r},s)=A_n{I}_0\left({\alpha }_n\overline{r}\right)+B_n{K}_0\left({\alpha }_n\overline{r}\right)+F_n{I}_0\left({\overline{\kappa }}_n\overline{r}\right)+G_n{K}_0\left({\overline{\kappa }}_n\overline{r}\right)-\frac{\Gamma }{s{\alpha }_n^2{\overline{\eta }}_n}.$$

(59)

The values of constants $A_n$ and $B_n$ were obtained by applying the boundary conditions given in (42)–(45) to (59), that with aid of (35), results in a system of linear equations. This system of equations can be represented by a general matrix notation as $\mathbf{A}\mathbf{x}=\mathbf{b}$ in Appendix B. Therefore, to find the values of the constants $A_n$ and $B_n$, the matrix inverse method is used in (A2).

Thus, once found the constants $F_n$, $G_n$, $A_n$ and $B_n$, these are replaced into (59). In this last equation, the inverse of the Laplace transform is applied using the numerical method based on concentrated matrix exponential (CME) distributions [89]. Here, a finite linear combination of the transform values approximates $\overline{u}$, as follows:

$${\overline{u}}_n(\overline{r},\overline{t})\approx {\overline{u}}_n(\overline{r},\overline{t},M)=\frac{1}{\overline{t}}\sum _{k=1}^M{\omega }_k{U}_n\left(\overline{r},\frac{{\theta }_k}{\overline{t}}\right),$$

(60)

where ${\omega }_1$ and ${\theta }_1$ are real coefficients, and from ${\omega }_2$ to ${\omega }_M$, and from ${\theta }_2$ to ${\theta }_M$ are $(M-1)/2$ complex conjugate pairs that the authors in Horváth et al. [89] provide. Here, $M=50$.

3.3. Steady-State Velocity

The steady-state solution for the velocity distribution is obtained simplifying (25) as follows:

$$\frac{1}{\overline{r}}\frac{d}{d\overline{r}}\left(\overline{r}\frac{d{\overline{u}}_n}{d\overline{r}}\right)=\frac{1}{{\overline{\eta }}_n}\left(\Gamma -{\overline{ϵ}}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\right),$$

(61)

which is integrated to yields

$$\begin{array}{c}\hfill {\overline{u}}_n=\frac{1}{{\overline{\eta }}_n}\left\{\frac{\Gamma {\overline{r}}^2}{4}+{\overline{ϵ}}_n{C}_{2n-1}\left[1-I_0\left({\overline{\kappa }}_n\overline{r}\right)\right]-{\overline{ϵ}}_n{C}_{2n}{K}_0\left({\overline{\kappa }}_n\overline{r}\right)+D_{n}ln\left(\overline{r}\right)+E_n\right\},\end{array}$$

(62)

where $D_n$ and $E_n$ are constants which are determined by applying the corresponding boundary conditions given in (27) and (30)–(32) into (62). The result of this application is a system of linear equations which can be represented by a general matrix notation as $\mathbf{A}\mathbf{x}=\mathbf{b}$ in Appendix C. Therefore, to find the constants $D_n$ and $E_n$, the matrix inverse method is used in (A4).

4. Results and Discussion

In this investigation, the following range of values are taken into account: $1\le R\le$ 50 $\mathsf{\mu }$m, $100\le {\kappa }_n^{-1}\le 300$ nm, ${\rho }_n=1000$ kg m${}^{-3}$, ${10}^{-4}\le {\eta }_n\le {10}^{-2}$ kg m${}^{-1}$ s${}^{-1}$, $E_z\le {10}^4$ V m${}^{-1}$, ${10}^{-11}\le {ϵ}_n\le {10}^{-9}$ C V${}^{-1}$ m${}^{-1}$, $Z_n=O\left(1\right)$, $0\le {\lambda }_n\le {10}^{-3}$ s, $-25\le {\zeta }_{\mathrm{I},\mathrm{O}}\le 25$ mV, $-6.25\le \Delta {\psi }_n\le 6.25$ mV, $k_B=$$1.381\times {10}^{23}$ J K${}^{-1}$ and $e=1.602\times {10}^{-19}$ C. Due to any pressure value can be provided by syringe pumps on the microfluidic devices, the pressure gradient, $p_z$, is considered of the same order as the electro-osmotic forces. The value of n will be kept in a moderate and realistic order of $O\left(1\right)$, although the solution of this problem can take values of $n=1,2,3,...,\infty$. Thus, with an appropriate combination of the aforementioned values, the following orders of magnitude of the dimensionless parameters are obtained: $\Gamma =\pm O\left(1\right),-0.25\le \Delta {\overline{\psi }}_n\le 0.25,10\le {\overline{\kappa }}_n\le 100,0.5\le {\overline{ϵ}}_n\le 2$, $0.1\le {\overline{\eta }}_n\le 12$ and $0\le {\overline{\lambda }}_n\le 10$. The density of fluids is assumed to be the same with ${\overline{\rho }}_n=1$.

4.1. Validation

This section aims to validate the semi-analytical solutions of the electro-osmotic flow presented in (60) and (62) for the transient and steady-state regime, respectively. Figure 2 compares the velocity profiles of a multilayer flow handling four layers with two investigations reported by the scientific community where a single fluid is transported. In all cases for the multilayer flow, the fluids are Newtonian with ${\overline{\lambda }}_n=0$ and have the same physical properties as ${\overline{\eta }}_n={\overline{ϵ}}_n={\overline{\rho }}_n=1$; also, the electrostatic interactions between fluid layers are absent in the liquid-liquid interfaces with $\Delta {\overline{\psi }}_n=0$, the value of the electrokinetic parameters is ${\overline{\kappa }}_n=10$, and any external pressure gradient is neglected with $\Gamma =0$. Furthermore, the specified zeta potential is ${\overline{\zeta }}_{\mathrm{I},\mathrm{O}}=1$ and because $\Delta {\overline{\psi }}_n=0$, there are no electric double layers at liquid-liquid interfaces. Under the conditions mentioned earlier for the multilayer flow and with a dimensionless internal radius of the channel of $a=0.5$, it is possible to recover and compare the solutions of the present work with the following two studies. The first work carried out by Chang and Wang [23] describes an analytical solution of the electro-osmotic flow in an annular microchannel. The transient regime presents the velocity profiles at the dimensionless times $\overline{t}=0.001,0.01,0.03,$ and 10. As a result, an excellent convergence can be observed between the transient solutions. On the other hand, the second work developed by Tsao [22], corresponds to an electro-osmotic flow through an annulus. In the steady-state, there is also a very good match between the solutions. This validation shows that when there are no electric double layers at the liquid-liquid interfaces and the physical properties of the fluids are equal, the velocity profile of a multilayer flow will behave like a single fluid layer.

4.2. Electric Potential and Velocity Profiles

Figure 3 shows the dimensionless electric potential and velocity profiles as a function of the radial coordinate. The multilayer electro-osmotic flow through the annular microchannel consists of three fluid layers equidistantly distributed with ${\overline{r}}_2=0.533$ and ${\overline{r}}_3=0.766$. The fluid layers 1 and 3 are in contact with the internal and external walls of the microchannel, respectively. Figure 3a–c also present the transient development of the velocity profiles as a function of dimensionless time $\overline{t}$ until reaching the steady-state regime in $\overline{t}\to \infty$, with a potential difference at the liquid-liquid interfaces of $\Delta {\overline{\psi }}_n=0.25,0$ and $-0.25$, respectively. The other dimensionless parameters are shown in the graph. Figure 3a,c show the electric potential distribution, where the values of $\Delta {\overline{\psi }}_n=0.25$ and $\Delta {\overline{\psi }}_n=-0.25$ produce a discontinuity of this distribution at the liquid-liquid interfaces. In addition, the values of $\Delta {\overline{\psi }}_n=0.25$ and $\Delta {\overline{\psi }}_n=-0.25$ also produce a change in polarity and distribution of the electric potential on each side of liquid-liquid interfaces. The aforementioned is due to the change in polarity of the electric charges in the electric double layers close to the liquid-liquid interfaces, which depends on the positive or negative values of $\Delta {\overline{\psi }}_n$. On the other hand, Figure 3b represents the classical electro-osmotic flow with a uniform electric potential distribution along the cross-section of the annular channel with $\Delta {\overline{\psi }}_n=0$. In this graph, it can be observed that for short times ($\overline{t}=0.05$), the movement of fluids due to electro-osmotic effects begins from the double electrical layers located on the channel walls. Then, as time passes, the momentum transfer out of the electrical double layers is transmitted through viscous effects towards the center of the microchannel. Also, since there are no discontinuities in the distribution of the electric potential at the liquid-liquid interfaces of Figure 3b, the shape of the velocity profiles is undisturbed. Conversely, for any discontinuity and change in polarity of the electric potential distribution at liquid-liquid interfaces as is observed in Figure 3a with $\Delta {\overline{\psi }}_n=0.25$ and in Figure 3c with $\Delta {\overline{\psi }}_n=-0.25$, the velocity profiles will suffer certain disturbances like velocity jumps around the region of contact between the fluids. Furthermore, the continuity of viscous stresses at the liquid-liquid interfaces is affected by the electrostatic effects from Maxwell stresses (see (31)), contributing to the disturbances in the shape of velocity profiles. Therefore, these disturbances in Figure 3a,c result from the presence of double electrical layers around the liquid-liquid interfaces. Here, the magnitude and sign of the potential difference $\Delta {\overline{\psi }}_n=0.25$ or $\Delta {\overline{\psi }}_n=-0.25$ can be controlled by the excess electrical charge in the interfacial region [72,87].

In Figure 4 it is presented the influence of the dimensionless relaxation times with ${\overline{\lambda }}_n=0.1,2.0$ and 10 on the velocity profiles. Figure 4a–c, present an oscillatory behavior of velocity due to the memory effects of fluids. These oscillations increase with the value of the dimensionless relaxation time. Thus, the velocity profile has a higher magnitude when the dimensionless relaxation time is ${\overline{\lambda }}_n=10$. Also, the multilayer flow of Figure 4c takes longer to reach the steady-state regime than the case shown in Figure 4a. Once the velocity profiles of Figure 4a–c reach the steady-state regime (see the blue lines for $\overline{t}\to \infty$), all of them acquire the same magnitude along the cross-section of the microchannel, which indicates that the elastic effects have ended. From this moment, the fluids behave as Newtonian fluids.

The velocity profiles in Figure 5 show the behavior of a purely electro-osmotic flow under the influence of the viscosity ratios. In Figure 5a–c, the viscosity ratio of each fluid layer is the same being ${\overline{\eta }}_n=0.7,2$ and 6, respectively. The viscosity value can measure the resistance that the multilayer flow opposes to flow. Therefore, Figure 5a denotes higher velocity profiles than Figure 5b,c for having the lowest viscosity value with ${\overline{\eta }}_n=0.7$ and, consequently, the lowest flow resistance.

4.3. Velocity Tracking

Once the behavior of the velocity profiles in the flow field has been described, the following results will analyze the velocity tracking at different positions of the annular microchannel cross-section as a time function. With this, the influence of dimensionless parameters on the flow dynamics in its transitory period will be elucidated.

To explore the damped oscillatory motion of Maxwell fluids, Figure 6 shows the velocity tracking as a time function of the three cases of multilayer flow presented in Figure 3 for $\Delta {\overline{\psi }}_n=0.25,0$ and $-0.25$ as case (a), (b) and (c), respectively. The velocity tracking position is located at the middle of each fluid layer, at $\overline{r}=0.4166,0.65$ and $0.8833$ for fluids 1, 2, and 3, respectively. As can be seen, fluids 1 and 3 in contact with the microchannel walls begin their movement before intermediate fluid 2. This effect happens because the electric potential in the walls has a greater magnitude than in any other position, which promotes a faster effect of the externally applied electric field on the double electric layers at the solid-liquid interfaces. After starting its movement, each fluid layer increases its velocity and experiences an oscillatory behavior that dampens with time until it reaches a steady-state regime. It can be seen that the crest and trough of waves in the oscillatory movement of fluids neighboring the microchannel walls present a velocity peak, while the central fluid does not. It is also observed that the value of the potential difference $\Delta {\overline{\psi }}_n$ influences the velocity magnitude in each fluid layer. Additionally, the potential difference modifies the time in which the multilayer flow reaches the steady-state (${\overline{t}}_{\mathrm{ss}}$), having for the case (a), (b) and (c) that ${\overline{t}}_{\mathrm{ss}}=22.35,24.48$ and $23.05$, respectively.

Figure 7 shows the velocity tracking of the multilayer flow presented in Figure 4. Tracking it is carried out in the middle of each fluid layer. It is noticeable that fluids with a higher dimensionless relaxation time take longer to come out of rest and reach the steady-state regime. Approximately, for a value of ${\overline{\lambda }}_n=10$ it takes a time of ${\overline{t}}_{\mathrm{ss}}=141.01$ to reach the steady-state regime, while for ${\overline{\lambda }}_n=0.1$ it needs ${\overline{t}}_{\mathrm{ss}}=2.09$. In this context, the fluids with a higher value of the dimensionless relaxation time have a more significant number of oscillations and a larger velocity magnitude than the other cases shown. The above demonstrates the elastic and memory effects of this type of Maxwell fluid and that they are intrinsic in the hyperbolic partial differential Equation (25) and its solution.

In other results, Figure 8 shows the velocity tracking of the analyzed cases in Figure 5 for values of ${\overline{\eta }}_n=0.7,2$ and 6, respectively, and in the middle of the annular microchannel at $\overline{r}=0.65$. As expected, the multilayer flow with a value of the viscosity ratios of ${\overline{\eta }}_n=6$ has a lower velocity magnitude than the other cases presented due to a higher resistance to flow. As observed, the high viscosity of the fluids will increase the number of oscillations in the transient period of the flow. In this figure, can noticed that the fluids with higher viscosity reach the steady-state regime faster; for example, with ${\overline{\eta }}_n=6$, the steady-state regime is achieved in ${\overline{t}}_{ss}=19.01$, while multilayer flow with ${\overline{\eta }}_n=0.7$ takes longer to reach that condition with ${\overline{t}}_{ss}=26.02$. The above is because by increasing the viscosity of the fluids, the flow resistance also increases, acting as a type of brake, damping the oscillatory regime faster. On the other hand, as previously mentioned, the presence of the elastic contribution of the fluid is manifested in the number of oscillations. That indicates stored energy coming from the microstructure immersed in the fluid, which gives rise to the viscoelastic behavior of Maxwell fluids. Then this type of fluids, when deformed, dissipate energy generating oscillations until reaching an equilibrium state or steady-state. Therefore, fluids with little resistance to flow dissipate stored energy more easily, showing an increase in velocity magnitude and few oscillations before reaching steady-state (case with ${\overline{\eta }}_n=0.7$). Conversely (case with ${\overline{\eta }}_n=6.0$), the high viscosity of the fluids will make it more difficult to dissipate the stored energy, increasing the number of oscillations and reducing the velocity magnitude.

Figure 9 represents the velocity tracking of a multilayer flow of three immiscible fluids as a function of dimensionless time. Here, four cases of combined values of the viscosity ratios are analyzed, case (a) with ${\overline{\eta }}_{1,2,3}=(1.0,1.5,2.1)$, case (b) with ${\overline{\eta }}_{1,2,3}=(3.5,1.9,0.9)$, case (c) with ${\overline{\eta }}_{1,2,3}=(0.9,4.0,2.1)$ and case (d) with ${\overline{\eta }}_{1,2,3}=(3.0,9.0,1.0)$. The other dimensionless parameters remain the same as Figure 8. It is noticeable that any combination of the viscosity ratios of the fluids presented in this figure, and which are in the same order of magnitude, that is, approximately $O\left(1\right)$, the time to reach the steady-state regime normalizes to an average value of ${\overline{t}}_{\mathrm{ss}}=16.75$. Under these conditions, the waves number along the oscillatory motion is normalized to a similar number in all the cases presented. Unlike Figure 8, the oscillations in this figure show more noise or peaks in their periods. In all cases, the movement of fluids in contact with the microchannel walls will begin before the intermediate fluid layer because the higher electro-osmotic effects will begin for early times in the walls than in the liquid-liquid interfaces. Once normalized and constrained the multilayer flow behavior of the three immiscible fluids, the total energy of each fluid is reflected in the velocity magnitude; in general, for the fluid with lower viscosity ratio it reaches a higher velocity.

Figure 10 shows the velocity tracking of a multilayer electro-osmotic flow of three fluid layers. In this figure, four cases (a), (b), (c), and (d) are studied with different combinations of the viscosity ratios with ${\overline{\eta }}_{1,2,3}=(10,0.1,5.0)$, ${\overline{\eta }}_{1,2,3}=(0.1,9.0,0.1)$, ${\overline{\eta }}_{1,2,3}=(0.1,1.0,12)$ and ${\overline{\eta }}_{1,2,3}=(12,1.2,0.1)$, respectively. The other dimensionless parameters are the same as Figure 8 and Figure 9. Unlike Figure 8 and Figure 9, in this figure, the combinations of the viscosity ratios handle values in three different orders of magnitude, that is, $O\left({10}^{-1}\right),O\left({10}^0\right)$, and $O\left({10}^1\right)$. It is clear that the assumed viscosity values have a relevant influence in each case and in the time in which it reaches the steady-state regime, being ${\overline{t}}_{\mathrm{ss}}=14.45,20.25,18.89$ and $20.35$ for the case (a), (b), (c) and (d), respectively. Case (a) of this figure shows that when more viscous fluids contact the microchannel walls, they will produce a more significant wave number during the transient period. On the contrary, in case (b), it is shown that when the less viscous fluids are in this position, they will produce a smaller number of oscillations to reach the steady-state regime. Intermediate cases of this behavior are (c) and (d) cases. It is important to note that assigning different orders of magnitude to the viscosity ratios in the immiscible fluids that compose a multilayer array, increases the nonlinearities (amplitude variations) in the oscillatory response of the flow.

Figure 11 analyzes the value of the dimensionless internal radius of the microchannel a on the velocity tracking. This figure shows the velocity behavior at the middle of each fluid for $a=0.1,0.4,0.6$ and $0.8$. The above means, from the case with a widest annular microchannel cross-section ($a=0.1$) to the most narrow ($a=0.8$), respectively. As can be seen, narrow annular channels develop the first wave of the oscillatory regime more quickly than the wider channels. On the other hand, the first wave with the highest velocity is developed in the widest channel and reduces its magnitude as the channel becomes narrower. In addition, as the channel is thinner, the number and uniformity of the waves increase with the value of the dimensionless internal radius (compare case $a=0.1$ with $a=0.8$). With the above mentioned, it is clear that when the value of the dimensionless internal radius changes, it will affect the velocity magnitude, number and uniformity of the waves during the transient period. In another result, it can be concluded that the influence of a is small on the time in which the multilayer flow reaches the steady-state regime because the values of ${\overline{t}}_{\mathrm{ss}}$ are very similar in the presented cases.

Figure 12 presents the influence of the remaining dimensionless parameters reported in this investigation on the velocity tracking in the multilayer electro-osmotic flow. In case (a) of this figure, it can be seen that the zeta potential values in the walls of the annular microchannel directly affect the velocity magnitude of the multilayer electro-osmotic flow. The highest velocity magnitude is reached during the transient and steady-state periods with symmetric zeta potentials on both walls. Whereas with asymmetric zeta potentials in the channel walls, i.e., with ${\overline{\zeta }}_{\mathrm{O}}<{\overline{\zeta }}_{\mathrm{I}}$ or ${\overline{\zeta }}_{\mathrm{O}}>{\overline{\zeta }}_{\mathrm{I}}$ the velocity magnitude will decrease. This behavior is because the value of the zeta potential in the walls defines the intensity of the electro-osmotic forces on the flow through the electric potential distribution into the electric double layers. In case (b), it is clear that by increasing the fluid’s electrical susceptibility via the ${\overline{ϵ}}_n$ parameter, the intensity of the applied external electric field forces and the flow velocity magnitude will increase. Case (c) shows the influence of the electrokinetic parameter. It is observed in this figure that steeper velocity tracking slopes are obtained by increasing the value of the electrokinetic parameter. The above is because in electro-osmotic flows with very thin electric double layers ($\overline{\kappa }\to \infty$), the elastic or restitution effects of fluids, also called memory effects, will cause much faster velocity changes during the oscillatory period of the flow. On the other hand, the imposition of a favorable ($\Gamma =2$) or adverse ($\Gamma =-2$) pressure gradient to the flow will increase or decrease the velocity magnitude of the multilayer flow concerning a purely electro-osmotic flow ($\Gamma =0$). In case (e), the number of fluid layers within the annular microchannel has little influence on the flow by slightly modifying its velocity magnitude. Any dimensionless parameters presented in this Figure 12 have minimal impact on the time the multilayer flow reaches a steady-state regime.

Finally, the dimensionless results presented in this work together with the suitable combination of the physical parameters given at the beginning of this section, comply with the restriction $\lambda \dot{\mathbf{\gamma }}<1$ [84] to use the linear viscoelastic model given in (4). Therefore, the above restriction in terms of the variables used in the present study is $\lambda (\partial u/\partial r)\sim \lambda (u_{\mathrm{HS}}/R)(\partial \overline{u}/\partial \overline{r})<1$.

5. Conclusions

The present work studies start-up of the multilayer electro-osmotic flow of Maxwell fluids through an annular microchannel. The obtained results show that:

• The Maxwell fluids exhibit an oscillatory behavior in the transient evolution due to the elastic and memory effects of this type of fluids.
• As the dimensionless relaxation time value increases, the number of oscillations, velocity magnitude, and the time to reach the steady-state regime will increase too.
• The dimensionless viscosity ratios of fluids dictate the degree of resistance they oppose to flow, which, combined with the dimensionless values of the relaxation times, strongly control the magnitude and oscillatory phenomenon of velocity field.
• While the dimensionless internal radius of the annular microchannel directly influences the number of oscillations presented by Maxwell fluids, the number of fluid layers does not affect the general behavior of the flow.
• Higher velocity magnitudes are obtained with symmetric zeta potentials at the microchannel walls. A similar effect is obtained on the velocity with higher dielectric permittivity values and the imposition of a favorable pressure gradient to the flow.