Start-Up Multilayer Electro-Osmotic Flow of Maxwell Fluids through an Annular Microchannel under Hydrodynamic Slip Conditions

Abstract

The present investigation analyzes the transient multilayer electro-osmotic flow through an annular microchannel with hydrophobic walls. The fluids are considered immiscible and viscoelastic, following the Maxwell rheological model. In the problem examined, the linearized Poisson–Boltzmann and Cauchy momentum equations are used to determine the electric potential distribution and the flow field, respectively. Here, different interfacial phenomena are studied through the imposed boundary conditions, such as the hydrodynamic slip and specified zeta potentials at solid–liquid interfaces, the velocity continuity, the electroviscous stresses balance, the potential difference, and the continuity of electrical displacements at the interfaces between fluids. The semi-analytic solution uses the Laplace transform theory. In the results, the velocity profiles and velocity tracking show the oscillatory behavior of flow, which strongly depends on the dimensionless relaxation time. Furthermore, the hydrodynamic slip on the channel walls contributes to the release of energy stored in the fluids due to elastic effects at the start-up of the flow. Similarly, other dimensionless parameters are also investigated. This research aims to predict the parallel flow behavior in microfluidic devices under electro-osmotic effects.

1. Introduction

Microfluidic devices and their applications are an active area of research, and pumping liquids within a fluidic network is essential to these devices [1]. In this context, electro-osmotic flow is the bulk liquid motion that results when an externally applied electric field interacts with the net surplus of charged ions in the diffuse part of an electric double layer [2]. This flow is commonly used as a pumping method in micro total analysis systems and lab-on-a-chip devices. Compared with other micropumps, electro-osmotic pumps offer several advantages, such as creating constant pulse-free flows and eliminating moving parts [3]. Furthermore, different applications have benefited from the electro-osmotic pumping method, such as high-flow-velocity applications [4], high-performance liquid chromatographic separations [5], electronics cooling [6], biological [7] and chemical [8] analysis systems, drug delivery [9], and cell culture [10].

Several investigations that involve different microchannel geometries have been reported in the literature to explain the electro-osmotic phenomenon. Concerning parallel-plates microchannels, the research by Yang et al. [11] showed the transient electro-osmotic flow of Newtonian fluids under the influence of the electric double-layer thickness and the geometric size of the microchannel. Zhao et al. [12] and Mahanta et al. [13] investigated the electro-osmotic flow characteristics of power-law and viscoelastic fluids through straight planar microchannels. In the case of rectangular geometries, Marcos et al. [14] undertook a theoretical study of the electro-osmotic flow under sinusoidally alternating and constant electric fields, and Miller et al. [15] analytically solved the start-up of the electro-osmotic flow, validating this experimentally via the velocity tracking located in the central region of the rectangular channels.

Another widely discussed geometry in the electro-osmotic flow is the cylindrical geometry, as represented by studies carried out by Keh and Tseng [16] and Kang et al. [17], investigating whether the flow field depends on the electrokinetic radius or the low and high zeta potentials. The electro-osmotic flow has also been studied in annular channels. Tsao [18] investigated the electro-osmotic flow in the steady-state, considering thin and thick electric double layers. Chang and Wang [19] studied the starting electro-osmotic flow in which the transient time depends on the inner-to-outer ratio for the annulus. More recently, Li et al. [20] reported the electro-osmotic flow of Jeffreys fluids in a microannulus, showing that, in particular, the wall gap reduction will restrict the fluid’s elastic effect. Also, the scientific community has extended its investigations with more complex geometries, such as elliptical [21], triangular [22], semicircular [23], and trapezoidal [24].

The interest in microscale fluid management has increased due to the development of new applications associated with chemical, biological, and medical areas that require the adequate control of fluid samples. Consequently, studying the rheology of fluids in electro-osmotic flows is fundamental in transporting biofluids, polymer solutions, or colloidal suspensions. In particular, fluids such as blood, DNA solutions [25], saliva, synovial fluid [26], mucus [27], and polymer solutions of polyethylene oxide [28], cetylpyridinium chloride, and sodium salicylate [29] are considered viscoelastic fluids studied by the scientific community under the rubric of microfluidics. Some of these fluids can be modeled as Maxwell fluids in the present investigation. In this context, Zhao et al. [12] analyzed the velocity profiles, flow rate, and shear stress distribution in an electro-osmotic flow, considering the viscosity’s shear-thinning and shear-thickening effects in power-law fluids. In other investigations, Wang et al. [30] and Escandón et al. [31] examined the viscoelastic effects of Maxwell fluids on the start-up of electro-osmotic flows, where, compared to Newtonian fluids, this kind of fluid produces a velocity overshoot and oscillations. In this context, other non-Newtonian fluids in electro-osmotic flows have been explored to highlight their rheological effect on the flow field, such as Burgers [32], Eyring [33], and Phan–Thien–Tanner fluids [34].

In the context of fluid manipulation in microfluidics and the hydrodynamic resistance that these offer when flowing in small conduits, the strategy of exploiting hydrophobic slippage arises [35,36]. Currently, the methods to prepare hydrophobic surfaces are the template method [37], the layer-by-layer self-assembly method [38], the electrospinning method [39], the chemical etching method [40], and the sol-gel method [41], among others [42,43]. Therefore, the benefits of electro-osmotic flow and hydrodynamic slip due to hydrophobic walls have been extensively reviewed. Ngoma and Erchiqui [44], Gaikwad and Mondal [45], and Chang et al. [46], among others, investigated the influence of slippage from hydrophobic microchannels on electrokinetic flows with Newtonian fluids. Additionally, other works have addressed the electro-osmotic flow with non-Newtonian fluids and hydrodynamic slip at the channel walls, such as those developed by Shit et al. [47] with power-law fluids, Awan et al. [48] with Oldroy-B fluids, Zhang et al. [33] with Eyring fluids, and Ferrás et al. [49] and Vaista et al. [50] with Phan–Thien–Tanner fluids. Also, references [51,52,53] present other approaches to address the slip problem for transporting non-Newtonian fluids in flow models.

On the other hand, in microfluidics, it is not only necessary to transport single fluids, but it is also necessary to pump several parallel layers of fluids without mixing to perform different analyses of samples, such as solvent extraction [54], chemical interactions at interfaces [55,56], nuclear chemistry [57], or to make polymer [58] and hydrogel [59] fibers. Previous works exploit the benefit of laminar flow at low Reynolds numbers. In this context, the scientific community has studied electro-osmotic pumps by viscous drag of non-conducting fluids through neighboring conducting fluids with two [60,61,62] and three layers [63] of immiscible fluids in microchannels. In other cases, all transported fluids are considered electrically conductive, providing additional electrostatic interactions at liquid–liquid interfaces for two [64,65,66] or three immiscible fluid layers [67]. Furthermore, the movement of several immiscible fluids under electro-osmotic effects has been extended to studies of multilayer flows to elucidate applications of the flow focusing and chemical processing of samples [68,69,70].

Therefore, the present investigation aims to increase the knowledge of transient multilayer electro-osmotic flows in an annular-type geometry. Also, it will highlight the interactions of the rheological behavior of Maxwell viscoelastic fluids with the hydrodynamic slip effects imposed on the hydrophobic walls of the microchannel, which have not been previously reported. For the above, a semi-analytical solution method based on the Laplace transform and a numerical method based on concentrated matrix exponential (CME) distributions will be used for the inversion of the Laplace transform.

2. Problem Formulation

2.1. Physical Model

The study investigates the start-up electro-osmotic flow exhibited by immiscible and incompressible fluids through an annular microchannel. In a multilayer pattern, fluids are considered viscoelastic, following the Maxwell rheological model. Regarding the geometry, two concentric tubes constitute the annular microchannel and, with the help of the microfluidic modules shown in Figure 1, allow the inclusion of several fluid layers in the range of $n=1,2,3,...,i$, where n is the fluid layer number. The cylindrical coordinates ($r,z$) shown by the enlarged view with a red square in Figure 1 represent the system’s origin in the center of the microfluidic device. For this configuration, $n=1$ and $n=i$ indicate the fluid layer in contact with the inner (at $r=aR$) and outer (at $r=R$) walls of the microchannel, respectively, where $0<a<1$ is a geometrical constant and R is the external radius. Each liquid–liquid interface in the annular cavity is located at $r_n$; the subscript range is $n=1,2,3,...,i-1$. Also, each fluid layer is composed of a mixture of a symmetrical electrolyte with a solute that provides the viscoelastic characteristics. Moreover, the immiscible fluids are considered electrically conductive, which causes electrostatic effects at the liquid–liquid interfaces via electric double layers with potential differences, $∆{\psi }_n$. Furthermore, the polarizable solid–liquid interfaces and zeta potentials ${\zeta }_{\mathrm{I}}$ and ${\zeta }_{\mathrm{O}}$ at the walls generate electric double layers in these regions. The movement of the fluids depends on the electro-osmotic forces produced by a uniform electric field $E_z$ that originates because a pair of electrodes induces an external electric potential at the ends of the conduits formed by the microfluidic modules. However, a pressure gradient $p_z$ in the z-direction can also influence the fluid layers. For this problem, it is considered that the channel walls are hydrophobic, generating a hydrodynamic slip velocity with a slip length $l_s$.

2.2. Governing and Constitutive Equations

The Poisson [61], continuity, and Cauchy momentum equations [65,71] govern this electro-osmotic flow for incompressible fluids with constant properties, given by

$${\nabla }^2{\psi }_n=-\frac{{\rho }_{e,n}}{{\epsilon }_n},$$

(1)

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

(2)

and

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

(3)

where $\psi$ is the electric potential within the EDL, ${\rho }_e$ is the volumetric free charge density, $\epsilon$ is the dielectric permittivity, $\mathbf{v}$ is the velocity vector, $\rho$ is the fluid density, t is the time, p is the pressure, $\mathit{\tau }$ is the stress tensor, and $\mathbf{E}$ is the electric field vector. The constitutive equations are the Boltzmann distribution for the free electric charge in the diffuse ionic layer adjacent to the charged surface [61]

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

(4)

and the linear rheological model for Maxwell fluids as [65,71,72]

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

(5)

In the two previous equations, it is found that, for each fluid layer n, Z is the valence of the electrolyte, e is the electron charge, $n_0$ is the ionic number concentration in the bulk solution, $k_{\mathrm{B}}$ is the Boltzmann constant, T is the fluid temperature, $\lambda$ is the relaxation time, ${\eta }_0$ is the zero-shear-rate viscosity, and $\dot{\gamma }=(\nabla \mathbf{v})+{(\nabla \mathbf{v})}^T$ is the rate-of-strain tensor.

For the mathematical modeling, transient, fully developed, and stratified laminar flow with very small Reynolds numbers is considered ($R{e}_n\ll 1$); here, $R{e}_n={\rho }_{n}R{u}_{\mathrm{c}}/{\eta }_{0,n}$, with $u_{\mathrm{c}}$ being the characteristic velocity. Also, it is established that the zeta potential is uniform throughout the study region of the channel walls. Taking into account the above, this work assumes that all liquid–liquid interfaces will be parallel to the microchannel walls. Additionally, the microchannel length given by L is larger than the external radius R ($R\ll L$) and the electric double layers do not overlap. Therefore, with the assumptions above, the governing and constitutive equations in cylindrical coordinates are simplified for a unidirectional flow [61,65,70] as follows: Considering the dimensionless variables given below

$$\overline{t}=\frac{{\eta }_{\mathrm{ref}}t}{{\rho }_{\mathrm{ref}}{R}^2},\overline{r}=\frac{r}{R},{\overline{\psi }}_n\left(\overline{r}\right)=\frac{{\psi }_n\left(r\right)}{{\zeta }_{\mathrm{ref}}},{\overline{u}}_n(\overline{r},\overline{t})=\frac{u_n(r,t)}{u_{\mathrm{c}}},{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})=\frac{R{\tau }_{rz,n}(r,t)}{{\eta }_{\mathrm{ref}}{u}_{\mathrm{c}}},$$

(6)

the following steps are given. First, for the electric potential from Equations (1), (4) and (6), the Poisson–Boltzmann equation is given as

$$\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).$$

(7)

In Equation (7), the Debye–Hückel approximation for low potentials in the range of $Z_{n}e{\psi }_n/k_{\mathrm{B}}{T}_n\le 1$ is considered, resulting in $sinh\left(Z_{n}e{\psi }_n/k_{\mathrm{B}}{T}_n\right)\approx Z_{n}e{\psi }_n/k_{\mathrm{B}}{T}_n$. Second, from Equations (3) and (6), the Cauchy momentum equation is rewritten for unidirectional flow as follows:

$${\overline{\rho }}_n\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{t}}=-\Gamma -\frac{1}{\overline{r}}{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})-\frac{\partial {\overline{\tau }}_{rz,n}(\overline{r},\overline{t})}{\partial \overline{r}}+{\overline{\epsilon }}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right).$$

(8)

And third, from Equations (5) and (6), the Maxwell rheological model is

$${\overline{\tau }}_{rz}(\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}}.$$

(9)

In Equations (7)–(9) the following dimensionless parameters appear:

$${\overline{\kappa }}_n=\frac{R}{{\kappa }_n^{-1}},{\overline{\rho }}_n=\frac{{\rho }_n}{{\rho }_{\mathrm{ref}}},\Gamma =\frac{p_z{R}^2}{u_{\mathrm{c}}{\eta }_{\mathrm{ref}}},{\overline{\epsilon }}_n=\frac{{\epsilon }_n}{{\epsilon }_{\mathrm{ref}}},{\overline{\eta }}_n=\frac{{\eta }_{0,n}}{{\eta }_{\mathrm{ref}}},{\overline{\lambda }}_n=\frac{{\eta }_{\mathrm{ref}}{\lambda }_n}{{\rho }_{\mathrm{ref}}{R}^2},$$

(10)

where, ${\overline{\kappa }}_n$ are defined as the ratios between the microchannel height to the Debye lengths ${\kappa }_n^{-1}={\left({\epsilon }_n{k}_{\mathrm{B}}{T}_n/2Z_n^2{e}^2{n}_{0,n}\right)}^{1/2}$, ${\overline{\rho }}_n$ are the densities ratios, ${\overline{\eta }}_n$ are the viscosity ratios, $\Gamma$ is the ratio of the external pressure forces to the electro-osmotic forces, and ${\overline{\lambda }}_n$ are the dimensionless relaxation times. Furthermore, the characteristic velocity $u_c=-{\epsilon }_{\mathrm{ref}}{\zeta }_{\mathrm{ref}}{E}_z/{\eta }_{\mathrm{ref}}$ is the Helmholtz–Smoluchowski velocity based in the reference properties.

However, for the momentum equation to remain only as a function of velocity, Equation (9) has to be derived with respect to the radial coordinate $\overline{r}$, yielding

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

(11)

Now, replacing Equation (11) into Equation (8) results in

$${\overline{\rho }}_n\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{t}}=-\Gamma +{\overline{\lambda }}_n\frac{{\partial }^2{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})}{\partial \overline{t}\partial \overline{r}}+{\overline{\eta }}_n\frac{{\partial }^2{\overline{u}}_n(\overline{r},\overline{t})}{\partial {\overline{r}}^2}-\frac{1}{\overline{r}}{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})+{\overline{\epsilon }}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right).$$

(12)

In addition to this, Equation (8) has to be derived in the time, resulting in

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

(13)

With the above, by replacing Equation (13) into Equation (12), the following relationship is obtained

$${\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 -\frac{1}{\overline{r}}{\overline{\tau }}_{rz,n}(\overline{r},\overline{t})-\frac{{\overline{\lambda }}_n}{\overline{r}}\frac{\partial {\overline{\tau }}_{rz,n}(\overline{r},\overline{t})}{\partial \overline{t}}+{\overline{\eta }}_n\frac{{\partial }^2{\overline{u}}_n(\overline{r},\overline{t})}{\partial {\overline{r}}^2}+{\overline{\epsilon }}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right).$$

(14)

Finally, dividing Equation (9) by the radial coordinate $\overline{r}$, and substituting the result in Equation (14), it is possible to remove the stress shear ${\overline{\tau }}_{rz,n}$, obtaining the following result:

$${\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 +\frac{{\overline{\eta }}_n}{\overline{r}}\frac{\partial {\overline{u}}_n(\overline{r},\overline{t})}{\partial \overline{r}}+{\overline{\eta }}_n\frac{{\partial }^2{\overline{u}}_n(\overline{r},\overline{t})}{\partial {\overline{r}}^2}+{\overline{\epsilon }}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right).$$

(15)

2.3. Boundary and Initial Conditions

Equations (7) and (15) to determine the electric potential and the velocity profiles require boundary conditions in $\overline{t}>0$ to be solved. At the solid–liquid interface placed at $\overline{r}=a$, the boundary conditions correspond to a specified zeta potential related to the EDL and the hydrodynamic slip given by the Navier slip condition [73,74,75], which must be applied to the first fluid layer $n=1$, as

$${\overline{\psi }}_1(\overline{r}=a)={\overline{\zeta }}_{\mathrm{I}},{\overline{u}}_1(\overline{r}=a,\overline{t})=-\frac{{\overline{l}}_{s,1}}{{\overline{\eta }}_1}{\overline{\tau }}_{rz,1}(\overline{r}=a,\overline{t}).$$

(16)

Then, at the liquid–liquid interfaces placed at $\overline{r}={\overline{r}}_{n=1,2,3,...,i-1}$, the corresponding boundary conditions to the electric potential related to EDLs are a potential difference and the continuity of electrical displacements in the absence of ions in the inner layer, respectively, as [65,76,77]

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

(17)

and

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

(18)

together with the boundary conditions for velocity, being the velocity continuity as follows:

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

(19)

and the electroviscous stress balance defined by [65]

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

(20)

Lastly, at the solid–liquid interface placed at $\overline{r}=1$, the boundary conditions correspond to a specified zeta potential related to the EDL and the Navier slip condition [73,74,75], which must be applied to the last fluid layer $n=i$, as

$${\overline{\psi }}_i(\overline{r}=1)={\overline{\zeta }}_{\mathrm{O}},{\overline{u}}_i(\overline{r}=1,\overline{t})=\frac{{\overline{l}}_{s,i}}{{\overline{\eta }}_i}{\overline{\tau }}_{rz,i}(\overline{r}=1,\overline{t}).$$

(21)

Complementarily, Equation (15) for all fluid layers requires initial conditions in the domain $a\le \overline{r}\le 1$ to be solved, with the velocity, shear stress, and acceleration being zero everywhere, that is, in the rest state, as follows:

$${\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.$$

(22)

In this section, the following dimensionless parameters arise

$${\overline{r}}_n=\frac{n}{i}(1-a)+a,∆{\overline{\psi }}_n=\frac{∆{\psi }_n}{{\zeta }_{\mathrm{ref}}},{\overline{l}}_{s;1,i}=\frac{l_{s;1,i}}{R},$$

(23)

where ${\overline{r}}_n$ are the liquid–liquid interface positions distributed equidistantly; however, for any other position ${\overline{r}}_n$ is designated arbitrarily according to the desired number of layers. $∆{\overline{\psi }}_n$ are the potential differences. ${\overline{r}}_n$ and $∆{\overline{\psi }}_n$ are in the range from $n=1$ to $n=i-1$. On the other hand, ${\overline{l}}_{s;1,i}$ are the ratios between the hydrodynamic slip lengths with the radius of the microchannel outer wall.

3. Solution Methodology

3.1. Electric Potential

Equation (7) recovers the modified Bessel differential equation whose solution is given by

$${\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),$$

(24)

where $C_{2n-1}$ and $C_{2n}$ are integration constants, while $I_0$ and $K_0$ are the modified Bessel functions of the first and second kind of zero-order, respectively. Equation (24) establishes the dimensionless electric potential, and its mathematical expression applies to any fluid layer. In particular, two constants C arise in this equation for each fluid layer. Therefore, the boundary conditions for each liquid–liquid and solid–liquid interface must be applied to determine the values of $C_{2n-1}$ and $C_{2n}$, which are given in Equations (16)–(18) and (21). Then, first the boundary condition in Equation (16) is applied into (24) for the fluid $n=1$ at the interface placed at $\overline{r}=a$, yielding

$$C_1{I}_0\left({\overline{\kappa }}_{1}a\right)+C_2{K}_0\left({\overline{\kappa }}_{1}a\right)={\overline{\zeta }}_{\mathrm{I}}.$$

(25)

Second, the boundary conditions in Equations (17) and (18) are applied to Equation (24) for any liquid–liquid interface placed at $\overline{r}={\overline{r}}_n$, yielding

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

(26)

and, respectively,

$${\overline{\epsilon }}_{n+1}{\overline{\kappa }}_{n+1}\left[C_{2(n+1)-1}{I}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)-C_{2(n+1)}{K}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)\right]={\overline{\epsilon }}_n{\overline{\kappa }}_n\left[C_{2n-1}{I}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)-C_{2n}{K}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)\right].$$

(27)

And third, the boundary condition in Equation (21) is applied into (24) for the fluid $n=i$ at the interface placed at $\overline{r}=1$, yielding

$$C_{2i-1}{I}_0\left({\overline{\kappa }}_i\right)+C_{2i}{K}_0\left({\overline{\kappa }}_i\right)={\overline{\zeta }}_{\mathrm{O}}.$$

(28)

The system of equations given by (25)–(28) can be rewritten in matrix form $\mathbf{Ax}=\mathbf{b}$, as presented in detail in Appendix A. Therefore, the inverse matrix method is used to find the constants C given by the vector $\mathbf{x}$.

3.2. Transient-State Velocity Profiles

The Laplace transform is applied to the velocity profile and shear stress to determine the velocity field 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}$$

(29)

and

$${\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}.$$

(30)

Therefore, by applying the above definitions to the constitutive and Cauchy momentum equations, that is, to Equations (9) and (15), the following relationships are obtained

$$\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{d^2{U}_n(\overline{r},s)}{d{\overline{r}}^2}+\frac{{\overline{\eta }}_n}{\overline{r}}\frac{d{U}_n(\overline{r},s)}{d\overline{r}}+\frac{{\overline{\epsilon }}_n{\overline{\kappa }}_n^2}{s}{\overline{\psi }}_n\left(\overline{r}\right)-\frac{\Gamma }{s}\end{array}$$

(31)

and

$${\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{d{U}_n(\overline{r},s)}{d\overline{r}}.$$

(32)

Once the initial conditions of Equation (22) are applied, Equations (31) and (32) are reduced as follows

$$\frac{d^2{U}_n(\overline{r},s)}{d{\overline{r}}^2}+\frac{1}{\overline{r}}\frac{d{U}_n(\overline{r},s)}{d\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}$$

(33)

and

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

(34)

where ${\alpha }_n^2=({\overline{\rho }}_{n}s/{\overline{\eta }}_n)({\overline{\lambda }}_{n}s+1)$, ${\beta }_n=-{\overline{\epsilon }}_n{\overline{\kappa }}_n^2/{\overline{\eta }}_{n}s$, and ${\gamma }_n={\overline{\eta }}_n/(1+{\overline{\lambda }}_{n}s)$. The Laplace transforms from Equations (29) and (30) are applied into Equations (16), (19)–(21) to establish the boundary conditions that allow solving Equation (33). In this direction, Equation (34) contributes to the previous solution process, and, as a result, for the interface placed at $\overline{r}=a$, the following expression is obtained

$$U_1(\overline{r}=a,s)=\frac{{\overline{l}}_{s,1}{\gamma }_1}{{\overline{\eta }}_1}\frac{d{U}_1(\overline{r}=a,s)}{d\overline{r}}.$$

(35)

Also, for each liquid–liquid interface placed at $\overline{r}={\overline{r}}_{n=1,2,3,...,i-1}$, the boundary conditions result in

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

(36)

and

$$-{\gamma }_{n+1}\frac{d{U}_{n+1}(\overline{r}={\overline{r}}_n,s)}{d\overline{r}}+\frac{{\overline{\epsilon }}_{n+1}}{s}\frac{d{\overline{\psi }}_{n+1}(\overline{r}={\overline{r}}_n)}{d\overline{r}}=-{\gamma }_n\frac{d{U}_n(\overline{r}={\overline{r}}_n,s)}{d\overline{r}}+\frac{{\overline{\epsilon }}_n}{s}\frac{d{\overline{\psi }}_n(\overline{r}={\overline{r}}_n)}{d\overline{r}},$$

(37)

and, finally, for the outer wall of the annular microchannel placed at $\overline{r}=1$, the following expression is obtained:

$$U_i(\overline{r}=1,s)=-\frac{{\overline{l}}_{s,i}{\gamma }_i}{{\overline{\eta }}_i}\frac{d{U}_i(\overline{r}=1,s)}{d\overline{r}}.$$

(38)

Consequently, Equation (33) and Equations (35)–(38) constitute the mathematical model for Maxwell immiscible fluids driven by the electro-osmotic effects, which is found in Laplace space. It should be noted that Equation (33) is an inhomogeneous ordinary differential equation (ODE) whose solution can be constructed by the sum of the homogeneous solution $U_{h,n}(\overline{r},s)$ with the particular solution $U_{p,n}(\overline{r},s)$, as shown below

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

(39)

Therefore, the corresponding homogeneous solution for Equation (33) is the modified Bessel differential equation, whose solution is

$$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),$$

(40)

while their particular solution takes into account the form of the electric potential solution and a constant for the pressure gradient term as follows:

$$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.$$

(41)

Here, $A_n$, $B_n$, $F_n$, $G_n$, and $H_0$ are integration constants that will be determined later. Therefore, by replacing Equation (41) into Equation (33), and with the aid of Equation (24), this results in

$$\begin{array}{c}\hfill \frac{d^2}{d{\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{d}{d\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}$$

(42)

Using the superposition principle of solutions, the following three relationships can be obtained from the above equation

$$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),$$

(43)

$$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)$$

(44)

and

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

(45)

On the other hand, when replacing Equation (24) in Equation (7), and following the principle of superposition of solutions, it 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)$$

(46)

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).$$

(47)

Therefore, $F_n$ and $G_n$ are calculated by substituting the equalities of Equations (46) and (47) into Equations (43) and (44), respectively, resulting in

$$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}.$$

(48)

With the above values of the constants given in Equations (45) and (48), the dimensionless velocity distribution for any fluid layer from Equation (39) takes the following form:

$$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}.$$

(49)

In order to obtain the values of the constants $A_n$ and $B_n$, the boundary conditions shown in Equations (35)–(38) have to be applied to Equation (49). First, by applying Equation (35) to Equation (49) corresponding to the fluid layer $n=1$ at the interface placed at $\overline{r}=a$, the following relationship is obtained:

$$\begin{array}{c}\hfill A_1\left[I_0\left({\alpha }_{1}a\right)-\frac{{\overline{l}}_{s,1}{\gamma }_1}{{\overline{\eta }}_1}{\alpha }_1{I}_1\left({\alpha }_{1}a\right)\right]+B_1\left[K_0\left({\alpha }_{1}a\right)+\frac{{\overline{l}}_{s,1}{\gamma }_1}{{\overline{\eta }}_1}{\alpha }_1{K}_1\left({\alpha }_{1}a\right)\right]=\\ \hfill -F_1\left[I_0\left({\overline{\kappa }}_{1}a\right)-\frac{{\overline{l}}_{s,1}{\gamma }_1}{{\overline{\eta }}_1}{\overline{\kappa }}_1{I}_1\left({\overline{\kappa }}_{1}a\right)\right]-G_1\left[K_0\left({\overline{\kappa }}_{1}a\right)+\frac{{\overline{l}}_{s,1}{\gamma }_1}{{\overline{\eta }}_1}{\overline{\kappa }}_1{K}_1\left({\overline{\kappa }}_{1}a\right)\right]+\frac{\Gamma }{s{\alpha }_1^2{\overline{\eta }}_1}.\end{array}$$

(50)

Second, by applying Equations (36) and (37) for the boundary conditions placed at $\overline{r}={\overline{r}}_{n=1,2,3,...,i-1}$ to Equation (49) yields

$$\begin{array}{c}\hfill A_n{I}_0\left({\alpha }_n{\overline{r}}_n\right)+B_n{K}_0\left({\alpha }_n{\overline{r}}_n\right)-A_{n+1}{I}_0\left({\alpha }_{n+1}{\overline{r}}_n\right)-B_{n+1}{K}_0\left({\alpha }_{n+1}{\overline{r}}_n\right)=\\ \hfill -F_n{I}_0\left({\overline{\kappa }}_n{\overline{r}}_n\right)-G_n{K}_0\left({\overline{\kappa }}_n{\overline{r}}_n\right)+F_{n+1}{I}_0\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)+G_{n+1}{K}_0\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)+\frac{\Gamma }{s}\left(\frac{1}{{\alpha }_n^2{\overline{\eta }}_n}-\frac{1}{{\alpha }_{n+1}^2{\overline{\eta }}_{n+1}}\right)\end{array}$$

(51)

and, respectively,

$$\begin{array}{c}\hfill A_n{\gamma }_n{\alpha }_n{I}_1\left({\alpha }_n{\overline{r}}_n\right)-B_n{\gamma }_n{\alpha }_n{K}_1\left({\alpha }_n{\overline{r}}_n\right)-A_{n+1}{\gamma }_{n+1}{\alpha }_{n+1}{I}_1\left({\alpha }_{n+1}{\overline{r}}_n\right)+B_{n+1}{\gamma }_{n+1}{\alpha }_{n+1}{K}_1\left({\alpha }_{n+1}{\overline{r}}_n\right)=\\ \hfill -{\gamma }_n{\overline{\kappa }}_n\left[F_n{I}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)-G_n{K}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)\right]+\frac{{\overline{\epsilon }}_n{\overline{\kappa }}_n}{s}\left[C_{2n-1}{I}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)-C_{2n}{K}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)\right]\\ \hfill +{\gamma }_{n+1}{\overline{\kappa }}_{n+1}\left[F_{n+1}{I}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)-G_{n+1}{K}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)\right]-\frac{{\overline{\epsilon }}_{n+1}{\overline{\kappa }}_{n+1}}{s}\left[C_{2(n+1)-1}{I}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)-C_{2(n+1)}{K}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)\right].\end{array}$$

(52)

And third, by applying Equation (38) to Equation (49) corresponding to the fluid layer $n=i$ at the interface placed at $\overline{r}=1$, the following relationship is obtained:

$$\begin{array}{c}\hfill A_i\left[I_0\left({\alpha }_i\right)+\frac{{\overline{l}}_{s,i}{\gamma }_i}{{\overline{\eta }}_i}{\alpha }_i{I}_1\left({\alpha }_i\right)\right]+B_i\left[K_0\left({\alpha }_i\right)-\frac{{\overline{l}}_{s,i}{\gamma }_i}{{\overline{\eta }}_i}{\alpha }_i{K}_1\left({\alpha }_i\right)\right]=\\ \hfill -F_i\left[I_0\left({\overline{\kappa }}_i\right)+\frac{{\overline{l}}_{s,i}{\gamma }_i}{{\overline{\eta }}_i}{\overline{\kappa }}_i{I}_1\left({\overline{\kappa }}_i\right)\right]-G_i\left[K_0\left({\overline{\kappa }}_i\right)-\frac{{\overline{l}}_{s,i}{\gamma }_i}{{\overline{\eta }}_i}{\overline{\kappa }}_i{K}_1\left({\overline{\kappa }}_i\right)\right]+\frac{\Gamma }{s{\alpha }_i^2{\overline{\eta }}_i}.\end{array}$$

(53)

The system of equations given by (50)–(53) can be rewritten in matrix form $\mathbf{Ax}=\mathbf{b}$, as is presented in detail in Appendix B. Therefore, the inverse matrix method is used to find the constants $A_n$ and $B_n$ given by the vector $\mathbf{x}$.

Once the values of $F_n$, $G_n$, $A_n$, and $B_n$ have been calculated, they are substituted in Equation (49) for evaluation. Then, the inverse of the Laplace transform is numerically applied to the results derived from the previous process by using the concentrated matrix exponential (CME) distributions [78], obtaining the approximate values of $\overline{u}$, as indicated below:

$${\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),$$

(54)

where, ${\omega }_1$ and ${\theta }_1$ are real coefficients; while, from ${\omega }_2$ to ${\omega }_M$, and from ${\theta }_2$ to ${\theta }_M$ are $(M-1)/2$ complex conjugate pairs given in [78]. Here, $M=75$. The above procedure is also applicable to the shear stress by changing ${\overline{u}}_n$ to ${\overline{\tau }}_{rz,n}$, and $U_n$ to ${\tilde{\tau }}_{rz,n}$ in Equation (54).

3.3. Steady-State Velocity Profiles

The simplified Equation (15) is used to determine the steady-state velocity profile as follows:

$$\frac{1}{\overline{r}}\frac{d}{d\overline{r}}\left(\overline{r}\frac{d{\overline{u}}_n\left(\overline{r}\right)}{d\overline{r}}\right)=\frac{1}{{\overline{\eta }}_n}\left[\Gamma -{\overline{\epsilon }}_n{\overline{\kappa }}_n^2{\overline{\psi }}_n\left(\overline{r}\right)\right].$$

(55)

Then, by integrating twice Equation (55) with respect to $\overline{r}$, yields

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

(56)

where, $D_n$ and $E_n$ are integration constants, which will be found by implementing the boundary conditions from Equations (16) and (19)–(21) to Equation (56). First, by applying Equation (16) to Equation (56) corresponding to the fluid layer $n=1$ at the interface placed at $\overline{r}=a$, the following relationship is obtained:

$$\begin{array}{c}\hfill D_1\left[\frac{a\ln \left(a\right)-{\overline{l}}_{s,1}}{a}\right]+E_1=\frac{\Gamma a}{4}\left(2{\overline{l}}_{s,1}-a\right)-{\overline{\epsilon }}_1\left\{C_1\left[-I_0\left({\overline{\kappa }}_{1}a\right)+{\overline{l}}_{s,1}{\overline{\kappa }}_1{I}_1\left({\overline{\kappa }}_{1}a\right)+1\right]-C_2\left[K_0\left({\overline{\kappa }}_{1}a\right)+{\overline{l}}_{s,1}{\overline{\kappa }}_1{K}_1\left({\overline{\kappa }}_{1}a\right)\right]\right\}.\end{array}$$

(57)

Second, Equations (19) and (20) for the boundary conditions placed at $\overline{r}={\overline{r}}_{n=1,2,3,...,i-1}$ are applied to Equation (56). To this, from Equation (9), it should be considered that, for the steady-state ${\overline{\tau }}_{rz,n}\left(\overline{r}\right)=-{\overline{\eta }}_n\left[d{\overline{u}}_n\left(\overline{r}\right)/d\overline{r}\right]$; therefore, yielding

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

(58)

and, respectively,

$$\begin{array}{c}\hfill D_n\left(\frac{1}{{\overline{r}}_n}\right)-D_{n+1}\left(\frac{1}{{\overline{r}}_n}\right)=2\left\{{\overline{\epsilon }}_n{\overline{\kappa }}_n\left[C_{2n-1}{I}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)-C_{2n}{K}_1\left({\overline{\kappa }}_n{\overline{r}}_n\right)\right]\right.\\ \hfill \left.-{\overline{\epsilon }}_{n+1}{\overline{\kappa }}_{n+1}\left[C_{2(n+1)-1}{I}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)-C_{2(n+1)}{K}_1\left({\overline{\kappa }}_{n+1}{\overline{r}}_n\right)\right]\right\}.\end{array}$$

(59)

And third, by applying Equation (21) to Equation (56) corresponding to the fluid layer $n=i$ at the interface placed at $\overline{r}=1$, results in

$$D_i{\overline{l}}_{s,i}+E_i=-\frac{\Gamma }{4}\left(2{\overline{l}}_{s,i}+1\right)+{\overline{\epsilon }}_i\left\{C_{2i-1}\left[I_0\left({\overline{\kappa }}_1\right)+{\overline{l}}_{s,i}{\overline{\kappa }}_i{I}_1\left({\overline{\kappa }}_i\right)-1\right]+C_{2i}\left[K_0\left({\overline{\kappa }}_i\right)-{\overline{l}}_{s,i}{\overline{\kappa }}_i{K}_1\left({\overline{\kappa }}_i\right)\right]\right\}.$$

(60)

The system of equations given by (57)–(60) can be rewritten in matrix form $\mathbf{Ax}=\mathbf{b}$, as is shown in detail in Appendix C. Therefore, the inverse matrix method is used to find the constants $D_n$ and $E_n$ given by the vector $\mathbf{x}$.

4. Results and Discussion

The results presented in this section were obtained by taking into account the following ranges of the physical parameters: $1\le R\le$50 μm, $50\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}$, $-25\le {\zeta }_{\mathrm{I},\mathrm{O}}\le 25$ mV, $Z_n\sim O\left(1\right)$, $0\le {\lambda }_n\le {10}^{-4}$ s, ${\epsilon }_n\sim O\left({10}^{-10}\right)$ C V${}^{-1}$ m${}^{-1}$, $-12.5\le ∆{\psi }_n\le 12.5$ mV, and $l_{s;1,i}\lesssim O($100 nm). Also, the constant values of $k_{\mathrm{B}}=1.381\times {10}^{23}$ J K${}^{-1}$ and $e=1.602\times {10}^{-19}$ C were considered. The reference properties are taken from an aqueous $1:1$ electrolyte (e.g., NaCl) as ${\epsilon }_{\mathrm{ref}}=7.083\times {10}^{-10}$ C V${}^{-1}$ m${}^{-1}$, ${\eta }_{\mathrm{ref}}=0.9\times {10}^{-3}$ kg m${}^{-1}$ s${}^{-1}$, ${\rho }_{\mathrm{ref}}=998$ kg m${}^{-3}$, $Z_{\mathrm{ref}}=1$, and $T_{\mathrm{ref}}=298$ K. With the above, ${\zeta }_{\mathrm{ref}}=k_{\mathrm{B}}{T}_{\mathrm{ref}}/Z_{\mathrm{ref}}e$ in mV. Therefore, the adequate combination of the physical parameters allows use of the following values of the dimensionless parameters: $-0.5\le ∆{\overline{\psi }}_n\le 0.5$, $10\le {\overline{\kappa }}_n\le 100$, $0.75\le {\overline{\epsilon }}_n\le 1.25$, $0.6\le {\overline{\eta }}_n\le 6$, $0\le {\overline{\lambda }}_n\le 10$, and $0\le {\overline{l}}_{s;1,i}\le 0.1$. Because any externally applied pressure can be reached by syringe pumps, in this work, the parameter $\Gamma$ will be maintained with a value comparable to the electro-osmotic effect, that is $\Gamma =O\left(1\right)$. Also, moderately and realistically, a finite number of fluid layers n with a value of $O\left(1\right)$ will be considered; however, the present solution is valid for any value of n. The fluids are stratified with similar densities, i.e., ${\overline{\rho }}_n=1$.

4.1. Validation

A comparison with published works was performed to verify the proper functioning of the analytical (steady-state regime) and semi-analytical (flow start-up) solutions that describe the electro-osmotic flow in an annular channel. In all the comparative cases, the percentage of the root mean square error (RMSE) was obtained between the solutions with the following:

$$\mathrm{RMSE}={\left[\sum _{p=1}^{100}{\left({\overline{u}}_{a,p}-{\overline{u}}_{b,p}\right)}^2/\sum _{p=1}^{100}{\overline{u}}_{a,p}^2\right]}^{1/2}\times 100,$$

(61)

where, ${\overline{u}}_{a,p}$ and ${\overline{u}}_{b,p}$ represent the velocity profile value arrays of the present work (“a”) and of the work with which it is compared (“b”), respectively. $p=1,2,...,m$ is the nodal counter; for the velocity profiles and shear stress distributions $m=100$, and for the velocity tracking $m=6600$. In the case of the shear stress distribution, the variable $\overline{u}$ is changed by ${\overline{\tau }}_{r,z}$ in Equation (61), following the same procedure.

4.1.1. Steady-State Regime

Figure 2 compares the results obtained in the steady-state regime of the flow field with the research performed by Tsao [18] on Newtonian fluids (${\overline{\lambda }}_n=0$). In this figure ${\overline{\kappa }}_n=10$; where, $n=1$ (one fluid layer) or $n=3$ (three fluid layers). The rest of the dimensionless parameters appear in the figure caption. The results in Figure 2a reflect velocity profiles that depend on the radial coordinate without hydrodynamic slip; that is, ${\overline{l}}_{s;1,i}=0$, corresponding to the results of Tsao [18] with one fluid layer (marks) against the present solution with three fluid layers (continuous line), where very good agreement between the solutions was found, with an $\mathrm{RMSE}=0.69\%$. Furthermore, Figure 2a shows the velocity profiles for a multilayer flow subject to hydrodynamic slip effects given by the parameter ${\overline{l}}_{s;1,i}$ (=0.01, 0.05, 0.1). Here, the hydrodynamic slip in the microchannel wall is evident, while the velocity magnitude increases with ${\overline{l}}_{s;1,i}$. Figure 2b develops the shear stress distribution through the annular microchannel for the values of Figure 2a. As can be seen, the stress distributions are independent of the variation in the parameter ${\overline{l}}_{s;1,i}$ since, in all cases, they overlap. This is acceptable since, in Figure 2a, velocity profiles with the same shape are shown, differing only in magnitude. The results in Figure 2b for the shear stress distribution with ${\overline{l}}_{s;1,i}=0$ from the present work compared with the one presented by Tsao [18], have an $\mathrm{RMSE}=2\%$. Additionally, to find the shear stress from Tsao [18], the following dimensionless velocity profile ${\overline{u}}_z=A{I}_0\left(\overline{\kappa }\overline{r}\right)+B{K}_0\left(\overline{\kappa }\overline{r}\right)+C\ln \left(\overline{\kappa }\overline{r}\right)+D$ was taken from this reference and replaced in Equation (9) for a Newtonian fluid; the values of $A,B,C$ and D were extracted from Tsao [18].

4.1.2. Flow Start-Up

In Figure 3, the transient flow of three Newtonian fluids driven by electro-osmotic effects in the annular microchannel was compared against that obtained by Chang and Wang [19] for one fluid layer without slippage (${\overline{l}}_{s;1,i}=0$). The comparison was performed in two different ways: firstly, by recovering the velocity profiles for the times ${\overline{t}}_{1,2,3}=0.001,0.01,0.03$ and ${\overline{t}}_4=10$ reported by Chang and Wang [19] in the bottom corner square of Figure 3, and, secondly, by performing the velocity tracking in the microchannel center (in the middle of the fluid layer 2), where the times reported in the velocity profiles are shown with white circles. Regarding the velocity profiles presented that are gradually increasing, there are values of ${\mathrm{RMSE}}_{1,2,3}=0.0104\%,0.0125\%,0.0237\%$, and ${\mathrm{RMSE}}_4=0.0001\%$, for the times ${\overline{t}}_{1,2,3}$ and ${\overline{t}}_4$, respectively. On the other hand, regarding the velocity tracking, which begins at $\overline{t}=0.0002$ and grows gradually until the steady-state is reached, an $\mathrm{RMSE}=0.0081\%$ is determined between solutions.

Figure 4 compares the velocity profiles and velocity tracking from this research with the work carried out by Escandón et al. [70] for a multilayer electro-osmotic flow of Maxwell fluids (${\overline{\lambda }}_n=1$) in an annular microchannel without slip on the walls. Here, a perfect convergence with $\mathrm{RSME}=0.0\%$ was found between the solutions. Additionally, in Figure 4, the results for the velocity profiles and velocity tracking with hydrodynamic slip derived from this research are shown (${\overline{l}}_{s;1,i}=0.05$), where significant changes in the velocity magnitude are observed with respect to the cases without slippage on the walls. At the start-up of the electro-osmotic flow, in all cases (${\overline{l}}_{s;1,i}=0$ and ${\overline{l}}_{s;1,i}=0.05$), the indicated times are ${\overline{t}}_{1,2,3}=0.01,0.6,1.2$, and ${\overline{t}}_4=80$ (already within steady-state), respectively. In particular, for the early time of ${\overline{t}}_1=0.1$, all the velocity profiles overlap—the one reported by Escandón et al. [70] without slippage, and the two in this work with ${\overline{l}}_{s;1,i}=0$ and ${\overline{l}}_{s;1,i}=0.05$. The above indicates that the hydrodynamic slip does not significantly affect the early stage of the start-up of the flow, as it also happens with the velocity tracking. The results obtained from Figure 4 also indicate that the presence of electrostatic effects ($∆{\overline{\psi }}_n=0.25$) at the interfaces between fluids causes protuberances in the velocity distribution in these regions, compared to the results achieved in Figure 2 and Figure 3, whose velocity distribution turns out to be more uniform. Also, it should be noted that the temporal evolution of the velocity profiles in Figure 4 (see for ${\overline{t}}_{1,2,3}=0.01,0.6,1.2$) is carried out with an oscillatory phenomenon around the velocity profile in steady-state (see for ${\overline{t}}_4=80$) due to the elastic characteristics of Maxwell fluids. In contrast, Figure 3 reflects a gradual growth over time for the evolution of the velocity profiles of Newtonian fluids.

With the above, it is concluded that the present solutions adequately describe, both qualitatively and quantitatively, the dynamics of the multilayer electro-osmotic flow studied here.

4.2. Velocity Tracking

Section 4.1 demonstrated the general behavior of fluids. Therefore, this section will focus on velocity tracking as a function of time for various positions within the annular microchannel. With this, it will be elucidated how the movement of fluids in the transient-state regime is affected by the dimensionless parameters, mainly through the effect of hydrodynamic slip.

Figure 5 and Figure 6 show the transport of three layers of fluids driven by electro-osmotic effects in an annular microchannel. In addition, the velocity tracking is developed on the solid–liquid interfaces placed at $\overline{r}=a=0.3$ and $\overline{r}=1$, and at the center of the fluid layers at $\overline{r}=0.417$, $\overline{r}=0.592$, and $\overline{r}=0.708$, respectively. The figure captions indicate the rest of the dimensionless parameters. Also, in Figure 5 and Figure 6, two panels of results are observed: the one on the left represents the velocity tracking of the electro-osmotic flow without slip effects on the walls (${\overline{l}}_{s;1,i}=0$), and the right panel represents this with hydrodynamic slip effects on the microchannel walls (${\overline{l}}_{s;1,i}=0.05$).

Figure 5a,b develop the velocity tracking for three values of the dimensionless relaxation time ${\overline{\lambda }}_n$ (=0.25, 2.5, 8). In all cases, by increasing the value of ${\overline{\lambda }}_n$, the elastic effects also increase, causing the behavior of viscoelastic fluids to approach the dynamics of a solid; that is, due to the elastic contribution, the system stores energy. Subsequently, this energy is dissipated through damped oscillations until a steady-state is reached. This behavior is more evident in intermediate fluid 2, where the fluid offers the least resistance to flow. In the case of Figure 5a, the oscillatory motion of the fluids tends to be symmetric around the equilibrium velocity condition. However, this kind of behavior disappears when the microchannel walls exhibit hydrodynamic slip, and the flow oscillations tend to dampen and develop below the equilibrium velocity condition, as presented in Figure 5b. In general, the time for the fluids to reach the steady-state (${\overline{t}}_{ss}$) is very sensitive to the value of the parameter ${\overline{\lambda }}_n$ since this increases considerably with its value. For Figure 5a, the values of ${\overline{\lambda }}_n$ (=0.25, 2.5, 8) produce times of ${\overline{t}}_{ss}=4.69,40.21,110.16$, respectively. However, for Figure 5b, these times are reduced to ${\overline{t}}_{ss}=3.85,31.38,89.09$, respectively, due to the effect of the hydrophobic walls. In other results, when comparing Figure 5a,b, in the case without slip (${\overline{l}}_{s;1,i}=0$), the velocity of fluids 1 and 3 at the solid–liquid interfaces is zero over time (no-slip condition); however, when slip is present (${\overline{l}}_{s;1,i}=0.05$), the velocity at these locations gradually increases, and the oscillations typical of Maxwell’s fluids seem almost imperceptible. From this discussion, it follows that although hydrodynamic slip increases the velocity magnitude, it has a damping effect on the oscillatory motion that propagates from the microchannel walls to the rest of the fluids.

In Figure 6a,b, the velocity tracking depends on the potential difference $∆{\overline{\psi }}_n$ (=−0.5, 0, 0.5). In particular, within the structure of the electric double layers at liquid–liquid interfaces, the potential difference represents the equilibrium value for the ionic activity across the diffuse layer [77]; for this problem, this parameter produces a change in the velocity magnitudes depending on their sign. This change in velocity is magnified at positions from the outer radius to the inner radius of the microchannel, thus having a higher velocity change in fluid 1. The above can be explained by the slightly antisymmetric behavior of the velocity profiles in an annular channel with a small inner radius, obtaining higher velocity magnitudes in positions close to the inner wall [19,20]. As can be seen in Figure 6a,b, the interfacial hydrodynamic slip condition at the microchannel walls with ${\overline{l}}_{s;1,i}=0.05$ produces a significant change in the velocity magnitude during the transient period and until reaching the equilibrium condition. In this sense, for Figure 6a, the average steady-state velocity is approximately 0.9 at the center of fluid layers 1, 2, and 3. In Figure 6b, an average velocity close to 2 indicates that the magnitude increased by more than double. Therefore, the hydrodynamic slip ${\overline{l}}_{s;1,i}$ has a more significant impact on the velocity magnitude and the time required for the flow to reach a steady-state than the interfacial potential difference $∆{\overline{\psi }}_n$.

Figure 7 presents an electro-osmotic flow of three layers of immiscible Maxwell fluids in an annular microchannel. Here, four values of the internal radius a (=0.15,0.35,0.65,0.85) determine the behavior of the velocity tracking. For this arrangement, the velocity tracking is performed halfway through each fluid layer, as indicated in the figure. Furthermore, Figure 7a,b represent the flow with walls without hydrodynamic slip (${\overline{l}}_{s;1,i}=0$) and with slip (${\overline{l}}_{s;1,i}=0.05$), respectively. It is observed that the channels with a wider cross-section ($a=0.15,0.35$) delay the onset of fluid movement. The aforementioned is because the transfer of momentum from regions of the electric double layers to the central part of each layer requires a longer time than in narrow channels ($a=0.65,0.85$). As seen in Figure 7a,b, higher velocity magnitudes are reached in the transient-state in the wider channels than the narrow ones. With this, it is confirmed that the elastic effects are restricted by reducing the gap between the walls of the annular channel [20]. Another effect of the gap between the channel walls is that during the transient flow period within wide channels ($a=0.15,0.35$), non-linear wave oscillations (complex shape) occur, which are damped over time until reaching the equilibrium condition (steady-state). On the other hand, when considering flows in narrow channels ($a=0.65,0.85$), the waves are linear (sinusoidal shape) and are of a greater number than in the case of wide channels. This is because the dissipation of energy stored by viscoelastic fluids confined in the widest channels is delayed and occurs suddenly, compared to what happens with fluids in narrow channels, where the energy dissipation is immediate and gradual. Additionally, as shown in Figure 7a, the reduction in the microchannel size causes the time required for the flow to acquire a steady-state to also tend to decrease, from the value of $a=0.15$ with ${\overline{t}}_{ss}=17.61$ to $a=0.8$ with ${\overline{t}}_{ss}=13.62$. However, this time reduction is accentuated in Figure 7b with $a=0.15$ with ${\overline{t}}_{ss}=14.23$ up to $a=0.85$ with ${\overline{t}}_{ss}=7.93$. With those mentioned earlier, it is concluded that the narrow channels release the energy stored by the Maxwell fluids faster than the wide channels. On the other hand, in Figure 7, the role of hydrodynamic slip is to increase the overall flow velocity, dampen the magnitude and number of the oscillations, and decrease the time required for the flow to acquire a steady-state, a behavior that can be verified by examining Figure 8 where the effect of hydrodynamic slip on the electro-osmotic flow is analyzed for three values of ${\overline{l}}_{s;1,i}$ (=0.01, 0.025, 0.1).

Figure 9 analyzes a multilayer electro-osmotic flow affected by the rest of the dimensionless parameters. The inner radius of the annular microchannel is $a=0.3$, and the velocity tracking is performed at half the annular space in all cases, at $\overline{r}=0.65$. The first dimensionless parameter analyzed in this figure is the zeta potential ${\overline{\zeta }}_{\mathrm{I},\mathrm{O}}$, and, as observed in Figure 9a, any combination of values in this parameter modifies the velocity magnitude because this parameter represents the impact of the electro-osmotic forces applied to the flow, achieving the highest velocity when there are symmetrical potentials on the walls, ${\overline{\zeta }}_{(\mathrm{I},\mathrm{O})}=(1,1)$. A difference in velocity due to the values of ${\overline{\zeta }}_{\mathrm{I},\mathrm{O}}$ is slightly accentuated when hydrodynamic slip occurs with ${\overline{l}}_{s;1,i}=0.05$ in Figure 9b. The second parameter analyzed is the permittivities ratio ${\overline{\epsilon }}_n$. The permittivity is a physical parameter of the fluids that describes how much they are affected by an electric field; therefore, when its value increases to ${\overline{\epsilon }}_n=1.25$, the flow will reach its maximum velocity, as observed in Figure 9a. Moreover, the combination of the values displayed for the permittivity ratio ${\overline{\epsilon }}_n$ with a hydrodynamic slip, increases the velocity, as shown in Figure 9b.

In this context, the third parameter ${\overline{\kappa }}_n$ shown in Figure 9a relates the characteristic scale of the annular microchannel to the electric double-layer thickness for each fluid layer. For wide channels and/or thin electric double layers from highly concentrated electrolytes, the movement of the fluids is delayed more, and its sudden start-up produces non-linear waves (for ${\overline{\kappa }}_n=40,90$), increasing the memory effects of the viscoelastic fluids. However, values of (${\overline{\kappa }}_n=15$) for narrow channels and/or thick electric double layers tend to produce waves with more linear behavior, and the movement of the fluids begins more quickly and gradually. The above is because, as seen in Figure 7, the elastic effects diminish as the gap between the walls of the annular microchannel decreases, which also happens for low concentration electrolytes (i.e., for values of ${\overline{\kappa }}_n=15$). Also, the combination between the parameter ${\overline{\kappa }}_n$ and the hydrodynamic slip of ${\overline{l}}_{s;1,i}=0.05$ is very relevant, as shown in Figure 9b, since it significantly modifies the velocity magnitude, increasing the velocity magnitude around 700% with respect to the case of ${\overline{\kappa }}_n=90$ in Figure 9a. The fourth dimensionless parameter analyzed in Figure 9 corresponds to the parameter $\Gamma$ (=3, 0, −3), which indicates the ratio between the external pressure forces to the electro-osmotic forces. Positive values of $\Gamma$ indicate that the pressure forces act in an opposite direction to the electro-osmotic effects, decreasing the velocity magnitude of the multilayer flow; conversely, negative values of $\Gamma$ act in the same direction, increasing the velocity. With $\Gamma =0$, this is the case where the only driving force is the electric field. By comparing Figure 9a,b, the hydrodynamic slip slightly affects the influence of $\Gamma$ on the flow.

The fifth parameter analyzed in Figure 9 very slightly affects the velocity magnitude of the fluids, for the cases of the number of fluid layers shown n (=3, 4, 5). This happens whether the walls are hydrophobic or not (see Figure 9a,b). However, it should be noted that, as the number of layers in the flow increases with $n=5$, the behavior of the waves will tend to be a little smoother with respect to the cases with $n=3$ and $n=4$. The aforementioned occurs because the individual widths of the fluid layers have decreased, and, as a result, there is a type of rectification of the elastic effects on the waveform in the oscillatory period of the flow. For the sixth and last dimensionless parameter presented in Figure 9, corresponding to the viscosity ratios, ${\overline{\eta }}_n$, the results show that the resistance to flow is higher for values of ${\overline{\eta }}_n=6$, producing lower velocity magnitudes and a greater number of oscillations in flow than is the case with ${\overline{\eta }}_n=0.6$. Also, more viscous Maxwell fluids dissipate the energy stored in the deformation process faster, as well as their equilibrium condition, acting as a brake. Furthermore, the combination of the hydrodynamic slip ${\overline{l}}_{s;1,i}=0.05$ with any value of ${\overline{\eta }}_n$, damps the number of oscillations, as shown in Figure 9a,b. Finally, for Figure 9, the time required for the flow to acquire a steady-state is slightly affected by most of the parameters, excluding the viscosity ratios that affect this time. As previously shown, the hydrodynamic slip does influence this time ${\overline{t}}_{ss}$, as can be seen by comparing Figure 9a,b.

Figure 10 evaluates the time to reach the steady-state for different dimensionless parameters in the presence and absence of hydrodynamic slip, as indicated in the figure caption. The ranges of the variables represent the general behaviors of the dimensionless parameters of this work. The results in Figure 10a suggest that the time to reach the steady-state significantly increases as the relaxation time rises, mainly when hydrodynamic slip is absent. However, this behavior contrasts with the time response to reach the steady-state to viscosity changes, which exhibits the opposite trend. The number of fluid layers and the external pressure slightly impact the time to reach the steady-state, as depicted in the figure. On the other hand, Figure 10b shows that increasing the annular gap of the microchannel reduces the time to reach the steady-state. Similarly, the time to reach the steady-state increases as the symmetrical wall zeta potentials and the electric permittivity ratios increase, and it is not significantly affected by variation in the electrokinetic parameter. In particular, Figure 10c indicates that the time to reach the steady-state gradually decreases due to the hydrodynamic slip condition. From these results and the previous figures, it can be concluded that the transient evolution provides a better understanding of the physical characteristics at the start-up of the flow; specifically, the time to reach the steady-state strongly depends on the relaxation times and partially depends on the rest of the dimensionless parameters. In this sense, analyzing the transient flow promotes the optimal and appropriate design of microfluidic platforms.

For the results of velocity tracking, a dimensionless time-step of 0.015 was taken into account. In Figure 3 and Figure 4, an extra evaluation of the velocity tracking was carried out in which only the times reported by Chang and Wang [19] and by Escandón et al. [70], respectively, were considered; these times are reported in these figures with marks in the form of white circles on the lines. The estimate of the time it takes for the flow to reach steady-state is obtained from the velocity tracking data by comparing the difference between two velocity points. This process starts from the final results of the velocity tracking in the steady-state and towards the flow oscillation condition, continuously comparing the velocity magnitudes over time until the velocity difference reaches the value of ≥10${}^{-5}$. Therefore, determining the time ${\overline{t}}_{ss}$, in which the flow acquires the steady-state, is referential since it depends on the time-step and the precision of the velocity difference between the two points considered. Thus, the general condition of the steady-state regime is when ${\overline{t}}_{ss}\to \infty$.

4.3. Shear Stress and Rate-of-Strain Distributions

Figure 11 shows the shear stress and rate-of-strain distributions from Figure 3, Figure 4 and Figure 5, respectively; all graphs are with respect to the radial coordinate and time. Figure 11a,b present the distributions of Figure 3, where electrostatic interactions at liquid–liquid interfaces ($∆{\overline{\psi }}_n=0$), fluid memory effects (${\overline{\lambda }}_n=0$), and the wall slip (${\overline{l}}_{s;1,i}=0$) are absent. The results are regular shape distributions of shear stress and rate-of-strain. Conversely, in Figure 11c–f, dealing with Maxwell fluids (with ${\overline{\lambda }}_n=1$ and ${\overline{\lambda }}_n=8$), as well as with electrostatic interactions at liquid–liquid interfaces ($∆{\overline{\psi }}_n=0.25$), distributions with irregular shapes are developed, showing a series of protuberances. As can be seen in Figure 11c–f, the profiles for shear stress and rate-of-strain overlap when ${\overline{l}}_{s;1,i}=0$ and ${\overline{l}}_{s;1,i}=0.05$, indicating that hydrodynamic slip will only affect the velocity magnitude (see also Figure 2 and Figure 4). Thus, the results shown in Figure 11 represent the general behavior of the shear stress and rate-of-strain distributions of the multilayer electro-osmotic flow analyzed here. Further, these results help to demonstrate that the linear viscoelastic Maxwell model given by Equation (5) is valid for this flow phenomenon since it satisfies the restriction that $\lambda {\dot{\gamma }}_n<1$ [72]. The previous restriction can be represented with the variables of the present study as $\lambda (du/dr)\sim \lambda (u_c/R)(d\overline{u}/d\overline{r})<1$.

4.4. Influence of the Type of Slip Boundary Condition

This section discusses two widely used models to represent the hydrodynamic slip boundary condition and its influence on this multilayer electro-osmotic flow field. The first model, also used in the present investigation, takes into account the general linear Navier slip condition for unidirectional flow as $u=(l_s/\eta )\left|{\tau }_{xy,wall}\right|$, where Newtonian and non-Newtonian fluid behavior characteristics are incorporated through the shear stress [73,75,79,80,81,82,83]; here, u is the axial velocity, $l_s$ is the slip length, $\eta$ denotes the viscosity, and ${\tau }_{xy,wall}$ is the shear stress at the wall. The second model arises when the shear stress represented by Newton’s law as ${\tau }_{xy,wall}=\eta \left|\frac{du}{dn}\right|$ is replaced in the above model to easily derive the equivalent model for the Navier slip condition for Newtonian fluids as $u=l_s\left|\frac{du}{dn}\right|$ (also known as the standard Navier slip condition), and which was used in studies of electro-osmotic flows in [44,45,46], or was employed to transport non-Newtonian fluids such as those reported by [33,47,48]. However, focusing on the context of non-Newtonian fluids in electro-osmotic flows, as the viscoelastic fluids treated in this work, as well as those reported by [49,50], the slip boundary condition takes the general form with the stress tensor $u=(l_s/\eta )\left|{\tau }_{xy,wall}\right|$. That is why it is essential to elucidate the impact of the fluid rheology through these two models of the slip condition on the behavior of the electro-osmotic flow.

In relation to the above discussion, Figure 12 presents the pumping of four layers of Maxwell fluids driven by purely electro-osmotic forces. The velocity tracking is shown with ${\overline{\lambda }}_n=1$ and ${\overline{\kappa }}_n=20$; other dimensionless values are placed in the figure caption. The velocity tracking is carried out in three positions: at the solid–liquid interfaces given by $\overline{r}=a$, $\overline{r}=1$, and $\overline{r}={\overline{r}}_2$ (central position in the annular microchannel). The electro-osmotic flow in Figure 12a details the velocity tracking considering the general Navier slip boundary conditions at the walls given by $U_{1,i}=\pm \left[{\overline{l}}_{s:1,i}/(1+{\overline{\lambda }}_{1,i}s)\right](d{U}_{1,i}/d\overline{r})$ in the Laplace transform space (from Equations (35) and (38)). This last equation represents the effect of the applied viscous stress when combining Equation (34) with Equations (16) and (21) in the Laplace transform space. While, the velocity tracking for Figure 12b uses the well-known standard or equivalent Navier slip condition, which in this work is recovered with the form $U_{1,i}=\pm {\overline{l}}_{s;1,i}(d{U}_{1,i}/d\overline{r})$ in the Laplace transform space. The oscillations in Figure 12a,b are a response to the elastic contribution of the fluids, indicating the stored energy coming from the microstructure immersed in the fluids, which gives rise to the viscoelastic behavior of Maxwell fluids [70]. However, it is clear when examining Figure 12a that the hydrodynamic slip condition in the walls, in combination with the fluid rheology, considerably dampens (decreases) the magnitude of the velocity in the oscillatory behavior, as well as the time required for the flow to acquire a steady-state regime ${\overline{t}}_{ss}$, with respect to the results shown in Figure 12b. In this comparison, the oscillations generated by the evolution of the fluid velocity at positions $\overline{r}=a=0.3$ and $\overline{r}=1$ in Figure 12b tend to disappear in Figure 12a. Therefore, the previous results regarding the electro-osmotic flow of Maxwell fluids indicate that the fluid velocity during the transient state and the time required for the flow to acquire a steady-state will be overestimated when using the standard model of the Navier equation ($U={\overline{l}}_s\left|(dU/d\overline{r})\right|$) compared to those obtained with the general model ($U=\left[{\overline{l}}_s/(1+{\overline{\lambda }}_{n}s)\right]\left|(dU/d\overline{r})\right|$).

5. Conclusions

The transient electro-osmotic flow of Maxwell immiscible fluids through an annular microchannel with hydrophobic walls was investigated in this work. The main findings are the following: The presence of hydrophobic walls in the microchannel increases the velocity magnitude of the fluids; however, it dampens the oscillatory phenomenon of the flow, partially decreasing the amplitude of the oscillations below the equilibrium condition. The gap of the annular channel and the electrokinetic parameter determine whether the waves of the oscillatory phenomenon in the flow will be linear or non-linear. In particular, the hydrodynamic slip effect drastically modifies the velocity of the fluid layers when the electrokinetic parameter is high. The dimensionless relaxation times play the dominant role in determining the time it takes for the flow to reach a steady-state, and they are followed in importance by the annular gap, the hydrodynamic slip lengths, the wall zeta potential, and the viscosity ratios. The above is because the relaxation times influence the elastic effects of the fluids, modifying the dynamic behavior through the energy stored in the system. Also, this work demonstrated that the appropriate selection of the slip condition is important for analyzing the transport of non-Newtonian fluids. On the other hand, the numerical solution of this study was successfully integrated with the previously solved analytical part, yielding validated and accurate results at a reasonable computational cost. Thus, the analytical-numerical model presented here is versatile and powerful, serving as a valuable tool for future research in resolving complex transient-state problems, particularly for very short evaluation times. Therefore, this investigation serves as a basis for designing microfluidic platforms with high spatio-temporal precision, facilitating various tasks related to fluid samples in medicine, biology, and chemistry. Finally, other implications not addressed in the present work, such as high zeta potentials, other rheological constitutive equations, variable physical properties of fluids, Joule heating due to high external electric fields, and deformation at the interfaces, can be considered in future work.