Start-Up Electroosmotic Flow of Multi-Layer Immiscible Maxwell Fluids in a Slit Microchannel

In this investigation, the transient electroosmotic flow of multi-layer immiscible viscoelastic fluids in a slit microchannel is studied. Through an appropriate combination of the momentum equation with the rheological model for Maxwell fluids, an hyperbolic partial differential equation is obtained and semi-analytically solved by using the Laplace transform method to describe the velocity field. In the solution process, different electrostatic conditions and electro-viscous stresses have to be considered in the liquid-liquid interfaces due to the transported fluids content buffer solutions based on symmetrical electrolytes. By adopting a dimensionless mathematical model for the governing and constitutive equations, certain dimensionless parameters that control the start-up of electroosmotic flow appear, as the viscosity ratios, dielectric permittivity ratios, the density ratios, the relaxation times, the electrokinetic parameters and the potential differences. In the results, it is shown that the velocity exhibits an oscillatory behavior in the transient regime as a consequence of the competition between the viscous and elastic forces; also, the flow field is affected by the electrostatic conditions at the liquid-liquid interfaces, producing steep velocity gradients, and finally, the time to reach the steady-state is strongly dependent on the relaxation times, viscosity ratios and the number of fluid layers.


Introduction
Microfluidics is a term that is used in fields of science with miniaturized systems where fluids are used as key components of control and sensing [1]. The handling of small sample volumes, scalability, integration of multiple functions and fields, low operating costs, low energy consumption, and so forth are some of the already known advantages of systems miniaturization; however, an inherent problem in these small devices is in the manufacture of moving parts for the manipulation of fluids or samples. That is why techniques based on electrokinetic effects arise that do not need moving parts [2]. The microsystem known as a laboratory on a chip, with size from millimeters to centimeters, facilitates the implementation of many laboratory tasks, which include sample preparation, mixing, separation, manipulation, control, detection, and analysis [3]. Applications cover the areas of mechanics, biology, chemistry, and medicine, seeking to improve technologies to preserve human health and improve the quality of life [4].
In this context, and to cover the different applications above mentioned, the electroosmotic flow has emerged as an electrokinetic phenomenon to transport fluids in microsystems. The electroosmosis represents the movement, due to an applied electric field, of an electrolyte solution relative to a stationary charged surface [5]. This transport method has been theoretically studied since many years ago in small channels with Newtonian fluids, as the early works carried out by Burgreen and Nakache [6] and Rice and Whitehead [7]. Since then and until now, the scientific community has continued studies on the electroosmotic flow behavior using Newtonian [8][9][10][11][12] and non-Newtonian fluids [13][14][15][16][17][18][19], addressing emerging issues regarding the rheology of the fluids and ionic concentrations, channel geometry, wall zeta potentials effects, boundary slip effects, among other topics, and their implications on the flow characteristics.
Research progress on electroosmotic flows with single-phase fluids is broad, as shown through the references cited in the previous paragraph and the references contained within them in turn. However, several microdevices applications use parallel multiphase flows for continuous chemical processing in analyses and synthesis [20,21], realizing specific operations as mixing and reaction [22], phase confluence and separation [23], solvent extraction [24,25], liquid-liquid extraction [26], purification [27] and synthesis of polymer membranes [28]. These parallel flows are a kind of flow pattern generated by using flow-focusing techniques in microsystems [29,30].
Following this direction, the scientific community exploits the benefit of electrokinetic phenomena by not including moving parts in microdevices, developing electrokinetic flow-focusing techniques. Hence, Pan et al. [31] perform an experimental investigation in microfluidic chips for mixing enhancement through the electrokinetic flow focusing and valveless switching of multiple samples flows. Meanwhile, Jiang et al. [32] propose a new microfluidic method to transport samples between sheath streams; here, the sheath streams for flow-focusing are generated by electroosmotic effects. In another work, Li et al. [33] present a theoretical and experimental investigation on flow-focusing with valveless switching, using the coupled effect of hydrodynamics and electroosmosis; in this work, the applied technique for switching a nonconducting sample stream or droplets, use two sheath streams of conducting fluids in a microchannel under electroosmotic effects. In this sense, Jia et al. [34] carried out a continuous-flow focusing study for collecting microparticles using induced-charge electroosmosis in a microfluidic device.
Therefore, for about two decades, the research to understand the physical mechanisms for moving parallel flows of immiscible fluids under electrokinetic effects has considered the transport of two layers [35][36][37][38][39][40][41][42] and three layers [43]. In these investigations, the arrangement of fluids considers that one fluid layer is non-conducting and the other(s) fluid(s) layer(s) if, being the electrolytic fluid(s) which is(are) under electroosmotic effects. Consequently, interfacial phenomena include the formation of an electrical double layer both in solid-liquid interfaces and in liquid-liquid interfaces. Regarding the liquid-liquid interfaces, the hydrodynamic and electrostatic boundary conditions are established in a relatively simple form with a specified zeta potential and the partially or null employment of Maxwell electric stress. However, other studies about the flow of two immiscible parallel fluids, consider that both fluids are conductive (i.e., fluids based in electrolytic solutions), increasing the complexity of the electrostatic boundary conditions in the liquid-liquid interface through a potential difference and the Gauss's law for the electrical displacement, together with the hydrodynamic boundary conditions via the combination of viscous and electric Maxwell stresses [44][45][46][47]. In addition, to cover the different flow-focusing applications in microdevices, the study of parallel flows under electrokinetic effects also has been extended to multi-layer systems [48][49][50].
To complement the research of Liu et al. [36], Li et al. [48], Afonso et al. [38], Huang et al. [39], Jian et al. [46], Matías et al. [42], Escandón et al. [49] about parallel flows with non-Newtonian fluids under electroosmotic effects, and considering that many of handled fluids in microsystems have complex rheological behavior, the present investigation aims to realize a parametric study on the start-up of the electroosmotic flow of multi-layer immiscible Maxwell fluids in a microchannel. The semi-analytical solution for the velocity field that is based on the Laplace transform method, obtains the description of new liquid-liquid interfacial phenomena as well as combined effects from the physical, rheological and electrical properties of fluids, over the velocity profiles and on the tracking of the velocity magnitude in the time. Hence, this model will be able to analyze many cause-effect relationships and it should serve as a guide for starting an experimental set-up in microdevices which require high spatio-temporal precision to carry out different clinical, chemical and biological analyses.

Physical Model Description
In this work, the transient electroosmotic flow of multi-layer immiscible fluids in a slit microchannel is analyzed. The physical model of the flow phenomenon studied here is on a Cartesian coordinate system (x, y), as is shown in Figure 1. The conduit is integrated by two parallel flat plates separated by a distance of 2H and filled by fluid layers, in which each one is composed of a mixture of symmetrical electrolytes with solutes that exhibit viscoelastic behavior. In the sketch, each liquid-liquid interface is placed in a y n position; where, the subscript n = 1, 2, 3, ..., i represents the number of the fluid layer, and i is the fluid layer in contact with the upper microchannel wall. Due to the fluids are immiscible and electrically conductive, in addition to the fact that the interfaces between them are polarizable and impermeable to charged particles, in this region, a double electrical layer appears presenting electrostatic properties through a potential difference ∆ψ n . The microchannel walls are also polarizable and acquire a surface electric charge represented by the zeta potential ζ w , also promoting the formation of an electric double layer in the solid-liquid interfaces. The fluids movement is due to the ends of the conduit are subject to an electric potential generated by a pair of electrodes that gives rise to a uniform electric field E x inducing electroosmotic effects on the electric double layers described before.

Governing and Constitutive Equations
The flow field of multi-layer immiscible fluids is governed by the Poisson equation for the electric potential distribution on the microchannel where Φ n is the electric potential, ρ e,n is the volumetric free charge density and n is the dielectric permittivity. Also, with the continuity equation for incompressible fluids as where v is the velocity vector. And the Cauchy momentum governing equation ρ n Dv n Dt = −∇p + ∇ · τ n + ρ n g + ρ e,n E, where ρ n is the fluid density, t is the time, p is the pressure, τ n is the stress tensor, g is the gravitational acceleration vector and E is the electric field vector. The shear stress tensor is related with the Maxwell rheological model as [51,52] τ n + λ n ∂τ n ∂t = −η 0,nγn , where λ n is the relaxation time, η 0,n is the zero-shear-rate viscosity andγ n = (∇v n ) + (∇v n ) T is the rate-of-strain-tensor.

Simplified Mathematical Model
The general governing equations given in the previous section can be simplified taking into account the following assumptions: • Constant physical properties and independent of the local electric field, ion concentration, and temperature [36,53].

•
Fully-developed flow [35]. • Impermeable interfaces between the fluids and ideally polarizable to electric charges [5,[54][55][56]. The electric double layers at the liquid-liquid interfaces are composed of two diffuse charge layers separated by a central compact layer, characterized by a potential difference drop due to the orientation of the solvent molecules; also, the continuity of electrical displacements on both sides of the central inner layer, in absence of ions in the inner layer is considered [5,[55][56][57][58].

•
Long microchannel neglecting any end effects; hence, the electric potential Φ n , is the algebraic sum of the potential due to the electric double layer, ψ n , and the potential due to the imposed electric field, φ, as [5]: where φ(x) = φ 0 − xE x ; being φ 0 the electric potential at the inlet of the microchannel at x = 0 and E x is the external electric field independent of the position and constant along the axial direction.

•
The local distribution of the free charges, that is, ions, is governed by the electrical potential into the electric double layer, ψ n , through the Boltzmann distribution as [5] ρ e,n = −2z n en 0,n sinh z n eψ n k B T n , 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 B is the Boltzmann constant and T n is the fluid temperature.

•
There is no imposed pressure gradient on microchannel.

•
The electric double layers do not overlap.

•
Any physical or chemical modification on the wall surfaces to cause hydrophobic interactions at the solid-liquid interfaces [47,64] is negligible.
As result, Equations (1)-(4), can be rewritten for a unidirectional flow as follows, leaving the Poisson-Boltzmann equation the momentum equation and respectively the Maxwell's constitutive rheological equation where κ 2 n = 2z 2 n e 2 n 0,n / n k B T n is the Debye-Hückel parameter related to the Debye length κ −1 n = n k B T n /2z 2 n e 2 n 0,n 1/2 [5], u n is the velocity on the x−direction and τ xy,n is the shear stress.
However, to get a momentum equation only in terms of velocity, Equation (9) is derived with respect to the transverse coordinate, producing which is replaced into Equation (8) obtaining On the other hand, Equation (8) is derived with respect to time, yielding The previous result is replaced into Equation (11), obtaining a momentum equation of hyperbolic type for the n fluids in terms of the axial velocity as follows To solve the governing equations given in Equations (7) and (13), the following boundary conditions in t > 0 for the electric potential and velocity are considered. At the bottom wall of microchannel for the fluid layer n = 1, the boundary conditions at y = 0 are a specified zeta potential and the no-slip boundary condition respectively as In the case of each liquid-liquid interface at y = y n=1,2,3,...,i−1 , the boundary conditions that are considered are a potential difference, the Gauss's law for the electrical displacement, a velocity continuity, and a stresses balance that include the Maxwell stresses and viscous shear stresses (also called electro-viscous stresses balance), respectively as follows and τ xy,n (y, t) + τ e,n (y) = τ xy,n+1 (y, t) + τ e,n+1 (y), (18) where the Maxwell shear stress is Additionally, at the top wall of microchannel for the fluid layer n = i, the boundary conditions at y = 2H, are a specified zeta potential and the no-slip boundary condition respectively as Finally, the initial conditions to solve the momentum equation, Equation (13), for the entire geometric domain 0 ≤ y ≤ 2H, that is, for all fluid layers are
The dimensionless parameters in this section are defined as whereκ n are the ratios between the microchannel height to the Debye lengths or also known as electrokinetic parameters,ρ n are the densities ratios,¯ n are the dielectric permittivities ratios,η n are the viscosity ratios andλ n are the dimensionless relaxation times. On the other hand,ȳ n are the interface positions and ∆ψ n are the potential differences; these two dimensionless parameters ranging from n = 1 to n = i − 1.

Electric Potential Distribution
The Poisson-Boltzmann equation for the electric potential, Equation (23), has a well-known solution that, in terms of n-layers of fluid is given bȳ where C 2n−1 and C 2n are integration constants that are determined by applying the corresponding boundary conditions at solid-liquid and liquid-liquid interfaces given in Equations (26)- (28) and (31)- (34). As result, the following equation system is obtained . . .
The above system of linear algebraic equations contains the same number of variables as the equations. Hence, the integration constants C 2n−1 and C 2n are solved by the matrix inverse method [65].
The homogeneous solution and the particular one, are the following equations, respectively U h,n (ȳ, s) = A n e α nȳ + B n e −α nȳ (47) and U p,n (ȳ, s) = D n eκ nȳ + E n e −κ nȳ (48) where A n , B n , D n , and E n are constant to be determined. The constants D n and E n are obtained by the substitution of Equation (48) into the Equation (40), yielding Considering the constant values of the previous equation and Equation (46), the dimensionless velocity distribution of each fluid layer of electroosmotic flow is U n (ȳ, s) = A n e α nȳ + B n e −α nȳ + D n eκ nȳ + E n e −κ nȳ .
To find the constants A n and B n , it is necessary to apply the boundary conditions given from Equations (42)- (45) to Equation (50), and with aid of Equation (34), the following system of linear algebraic equations is obtained . . .
which has been solved using the inverse matrix method in a process analogous to that of the electric potential distribution. Therefore, the constants D n and E n in Equation (49), and the constants A n and B n found through Equation (51), are replaced into Equation (50), where the inverse Laplace transform is numerically applied to solve the velocity distribution in this electroosmotic flow. To this, the method based on concentrated matrix exponential (CME) distributions is used [66]; in this framework, a finite linear combination of the transform values approximatesū, viā where ω 1 and θ 1 are real coefficients, and from ω 2 to ω M , and from θ 2 to θ M are (M − 1)/2 complex conjugate pairs that the authors in Horváth et al. [66] provide. Here, M = 50.

Solution Validation
To validate the performance of the semi-analytical solution found in this work for the transient velocity distribution, a comparison was made with two investigations reported by the scientific community, considering the transport of Newtonian and Maxwell fluids, respectively. In the first case, in the research carried out by Yang et al. [10], they model an electroosmotic flow of an aqueous 1:1 electrolyte (NaCl) in a slit microchannel with the following physical properties: a density of ρ = 998 kg m −3 , a viscosity of η=0.90×10 −3 kg m −1 s −1 , a relative electrical permittivity of r = 80, and a concentration of 10 −4 M, at a temperature of T = 298 K; additionally, the microchannel size and the wall zeta potential were set at 2H = 10 µm and 50 mV, respectively. With that set of values, the following electrokinetic parameter is obtained,κ n =164.5, and the dimensionless times to evaluate the velocity profiles aret = 0.0036, 0.036, 0.36 and 3.6 (=0.1, 1, 10 and 100 µs). Therefore, by comparing the work of Yang et al. [10] with the present investigation for three immiscible fluid layers, in Figure 2, an excellent convergence between their results is shown.  [10] with n = 1, against the results of the present investigation with three fluid layers, n = 3 (ȳ 1 =2/3 and y 2 =4/3). The other parameters areρ n =η n =¯ n =1, and ∆ψ n =0.
In the second case, the analytical solution for the dimensionless velocity profiles obtained by Escandón et al. [16] on the transient electroosmotic flow with Maxwell fluids, are compared with the present study, as shown in Figure 3. Here, the electrokinetic parameter takes a value ofκ n = 20 and the viscoelastic behavior of fluids is presented trough the two dimensionless relaxations time values of λ n = 0.12 andλ n = 2.5 in Figure 3a,b, respectively, finding a very good match between the results.   [16] with n = 1, against the results of the present investigation with four fluid layers, n = 4 (ȳ 1 =1/2, y 2 =1.0, andȳ 3 =3/2). The other parameters areρ n =η n =¯ n =1, and ∆ψ n =0, for: (a)λ n = 0.12 and (b)λ n = 2.5. Figure 4 shows the dimensionless electric potential profile,ψ n , and the start-up of the electroosmotic flow velocity,ū n , of three layers (n = 3) of immiscible Maxwell fluids in a slit microchannel, both variables as a function of the transverse coordinateȳ and for three different values of the potential difference ∆ψ n (= 0.5, 0, −0.5). The interfaces between fluids have been placed inȳ 1 = 2/3 andȳ 2 = 4/3, respectively, and the other dimensionless parameters are specified in the caption of the figure. Because of the high ionic concentration in the diffuse layers within the electric double layers formed in the solid-liquid interfaces in the system, there is a higher magnitude of the electric potential in these zones. On the other hand, in Figure 4a with ∆ψ n = 0.5 and Figure 4c with ∆ψ n = −0.5, can be appreciated an electric slip at liquid-liquid interfaces due to counterions concentrations in each side of interfaces, while for Figure 4b with ∆ψ n = 0, the classical null distribution of electrical potential is recovered outside of the electrical double layers on the walls. Here, the potential difference or electric slip between immiscible layers is proportional to the difference in the magnitude of ∆ψ n given by Equation (27) at each interface, and the sign gives the orientation of the counterions and electric potential distribution. Regarding the velocity, in each Figure 4a-c, are shown the evolution of the velocity profiles since an early time oft = 0.05 to the steady-state whent → ∞. As can be seen, for the early times, the fluids movement beginning from the Debye length in the solid-liquid and liquid-liquid interfaces due to electroosmotic effects, transmitting the movement by viscous forces to the rest of the fluid layers as time progresses. The influence of the potential difference on velocity development is clear when comparing Figure 4a,c with Figure 4b, producing great disturbances and steep velocity gradients in the flow velocity. Figure 5 shows the elastic behavior of the Maxwell fluids via the dimensionless relaxation time on the flow dynamics. In this figure are presented three cases for the velocity profiles evolution, in Figure 5a-c, the selected dimensionless relaxation times values areλ n = 0.1,λ n = 2 and λ n = 10, respectively. In these figures, it is noticeable that the start-up of fluids is more slowly as the dimensionless relaxation time increases, due to the memory effects of the viscoelastic fluids also increase, delaying the start of movement of fluids. Hence, in the case of Figure 5a the time to reach the steady-state whent → ∞ is much shorter than Figure 5c; in this context, the velocity magnitude and oscillatory behavior are greater for the case withλ n = 10 than forλ n = 0.1, because there is a more severe competition between the viscous and the elastic forces in the first case, that is, withλ n = 10. Figure 5a-c, reverse flow is produced in early times due to the potential difference at liquid-liquid interfaces.

For all cases in
In Figure 6 is represented the evolution of the velocity profiles of three fluid layers as a function of the transverse coordinate, and three combinations of the viscosities ratiosη n (= 0.7, 2, 6). It is evident that in Figure 6a, being the case with the lower viscosity fluid layers withη n = 0.7, the flow has the highest magnitude of velocity profiles; this can be corroborated if it is compared with the case of Figure 6c with a viscosity ofη n = 6 where the fluids have a greater resistance to flow and a lower magnitude of velocity. In all Figure 6a-c, the oscillatory behavior from elastic effects of Maxwell fluids is maintained. Figure 7 shows a wide combination of all dimensionless parameters studied in the present work. Here is observed the response of the dimensionless electric potential and velocity profiles for the electroosmotic flow of four layers (n = 4) of immiscible Maxwell fluids, under different values of electrokinetic parameters (κ n=1,2,3,4 = 10, 20, 30, 40) and relaxation times (λ n=1,2,3,4 = 0.1, 2, 1, 10). The velocity evolution goes from the timet = 0.1 to the steady-state whent → ∞. It can be seen a constant and gradual velocity evolution for the fluid 1, due to the small value of the relaxation time while for a higher relaxation time value the velocity profile oscillates continuously until reaching the steady-state resulting in stronger memory effects from the viscoelastic fluid such is the case of the fluid 4. Regarding the electrokinetic parameter effect, a value ofκ 1 = 10 in fluid 1 yields a parabolic shape of the velocity profiles due to the low ionic concentration of the buffer solution, producing a thick electric double layer; contrary the aforementioned, as the electrokinetic parameter grows in fluid 4 to take the value ofκ 4 = 40, the electric double layer becomes thinner and results in more slanted and straighter velocity profiles.

Tracking of the Velocity
In Figure 8 the tracking of the dimensionless velocity on the transverse coordinateȳ = 1 as a function of the dimensionless time is presented. In all cases of the sub-figures contained in Figure 8, the start-up of movement of fluid(s) in the center of the microchannel is delayed by memory effects due to its viscoelastic properties. In general, after overcoming these memory effects, a sudden and severe increase in velocity occurs, beginning a continuous oscillatory movement of increasing and decreasing velocity until steady-state is reached. The results presented in Figure 8a-c are taken from Figures 4-6, respectively. It is clear from Figure 8a,d,f, that the time it will take for the fluids to reach the steady-state is independent of the dimensionless parameters ∆ψ n ,ρ n andκ n , respectively. However, from Figure 8b, the time to reach the steady-state in the multi-layer electroosmotic flow is strongly dependent of the relaxation time, where forλ n = 0.1,λ n = 2 andλ n = 10 the dimensionless times to reach the steady-state aret ss ≈ 2.27,t ss ≈ 36.36 andt ss ≈ 181.82, respectively, these time results are due to the increases of fluid elasticity via the parameterλ n . In this context, from Figure 8c the time to reach the steady-state in the flow is also dependent of the viscosity ratiosη n , being the less viscous fluids those that take longer to reach that regime withη n = 0.7 in a time oft ss ≈ 20, while on the contrary case, withη n = 6, the steady-state is reached in a shorter time witht ss ≈ 12.73, due to the increase of the viscous forces and the corresponding faster braking of the flow. Furthermore, it can be seen in Figure 8e that multi-layer flows with three or more fluids take less time to reach the permanent regime due to the combined effects at the liquid-liquid interfaces. Finally, the time to reach the steady-state regime is established for the present work as the time in which the absolute value of the velocity difference between two immediate times at the same position is less than 10 −3 .

Conclusions
In this investigation, a semi-analytical solution of the start-up of the electroosmotic flow of multi-layer immiscible Maxwell fluids in a slit microchannel was obtained. In the parametric study, different fluid properties, geometric characteristics of the number and thickness of fluid layers, and electrostatic boundary conditions at liquid-liquid interfaces were considered. The electrostatic conditions from the electric double layers between the fluids via the potential differences and electro-viscous shear stresses, break the continuity of the electric potential distribution and produce significant changes in the velocity profiles in these zones. Regarding the dimensionless relaxation time effects on the velocity profiles, as this parameter increases a longer oscillatory behavior is caused by the memory effects of Maxwell fluids. Likewise, the magnitude of the flow velocity will significantly reduce with layers of more viscous fluids due to greater resistance to flow. In other results on the fluid dynamics of the multilayer electroosmotic flow, the time to reach the steady-state regime is strongly controlled by some dimensionless parameters reported here, like the relaxation times, the viscosity ratios, and the number of fluid layers.