*Article* **Electroosmotic Flow of Viscoelastic Fluid through a Constriction Microchannel**

**Jianyu Ji <sup>1</sup> , Shizhi Qian 1,\* and Zhaohui Liu <sup>2</sup>**


Huazhong University of Science and Technology, Wuhan 430074, China; zliu@mail.hust.edu.cn

**\*** Correspondence: sqian@odu.edu; Tel.: +1-757-683-3304

**Abstract:** Electroosmotic flow (EOF) has been widely used in various biochemical microfluidic applications, many of which use viscoelastic non-Newtonian fluid. This study numerically investigates the EOF of viscoelastic fluid through a 10:1 constriction microfluidic channel connecting two reservoirs on either side. The flow is modelled by the Oldroyd-B (OB) model coupled with the Poisson–Boltzmann model. EOF of polyacrylamide (PAA) solution is studied as a function of the PAA concentration and the applied electric field. In contrast to steady EOF of Newtonian fluid, the EOF of PAA solution becomes unstable when the applied electric field (PAA concentration) exceeds a critical value for a fixed PAA concentration (electric field), and vortices form at the upstream of the constriction. EOF velocity of viscoelastic fluid becomes spatially and temporally dependent, and the velocity at the exit of the constriction microchannel is much higher than that at its entrance, which is in qualitative agreement with experimental observation from the literature. Under the same apparent viscosity, the time-averaged velocity of the viscoelastic fluid is lower than that of the Newtonian fluid.

**Keywords:** electroosmosis; microfluidics; elastic instability; non-Newtonian fluid; Oldroyd-B model

**Citation:** Ji, J.; Qian, S.; Liu, Z. Electroosmotic Flow of Viscoelastic Fluid through a Constriction Microchannel. *Micromachines* **2021**, *12*, 417. https://doi.org/10.3390/mi 12040417

Academic Editor: Stéphane Colin

Received: 6 February 2021 Accepted: 7 April 2021 Published: 9 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

## **1. Introduction**

Electroosmotic flow (EOF) uses electric field to control fluid motion, and has been widely used in various microfluidic and nanofluidic applications such as fluid pump [1], mixing [2], and polymer translocation in biosensing [3]. The existing studies of EOF have been mainly focusing on Newtonian fluids [4,5]. However, in reality, EOF has been widely used to control and manipulate biological fluids (e.g., blood, saliva, lymph, protein, and DNA solutions) [6–8] and polymeric solutions [9], which exhibit non-Newtonian characteristics. Therefore, investigating EOF of viscoelastic fluids is of practical importance.

Bello et al. [10] conducted the pioneering study on EOF of non-Newtonian fluid, and measured EOF velocity of methyl cellulose solution in a capillary. Their results show that EOF velocity of such polymer solutions is much higher than that predicted with the classic Helmholtz-Smoluchowski velocity. Chang and Tsao [11] conducted similar experiments and found the effective viscosity decreased because of the sheared polymeric molecules inside the electrical double layer (EDL). Theoretically, non-Newtonian effects can be characterized by proper constitutive models relating the dynamic viscosity and the rate of shear. Such constitutive models include power-law model [12], Carreau model [13], WhiteMetzner model [14], Bingham model [15], Oldroyd-B (OB) model [16], PTT model [17], Moldflow second-order model [18], Giesekus model [19], etc. Das [20] developed an approximate solution for EOF velocity of power-law fluid between two parallel plates. Zhao et al. [21,22] derived a generalized Helmholtz–Smoluchowski velocity for EOF of power-law fluid in a slit microchannel. Later, Zhao and Yang [23,24] extended the study to a cylindrical microcapillary. Olivares et al. [25] experimentally investigated

EOF of a non-Newtonian polymeric solution and verified the generalized Helmholtz– Smoluchowski velocity. Tang et al. [26] numerically investigated EOF of power-law fluid using Lattice–Boltzmann method. Zimmerman et al. [13] carried out numerical simulation of EOF of Carreau fluid in a T-junction microchannel, and found that the flow field significantly depended on the non-Newtonian characteristics of the fluid. The aforementioned studies on EOF of non-Newtonian fluid are limited to inelastic constitutive models (i.e., power-law and Carreau models). However, some fluids show both viscous and elastic behaviors, which can be presented by viscoelastic constitutive models. There are existing literatures investigating the characteristics of EOF of viscoelastic fluids [27–31], showing that the viscoelasticity of the fluid affects the flow pattern and flow rate. Note that in the aforementioned studies, the EOF of non-Newtonian fluid was assumed in a steady state.

Recently, EOFs of non-Newtonian fluids have been reported to be time-dependent and show instabilities even at low Reynolds number (*Re*). Such EOFs are time-dependent because of the nonlinear viscosity and elasticity of non-Newtonian fluids. Bryce and Freeman [32] first reported the electro-elastic instability in EOF of PAA solutions through a 2:1:2 micro-scale contraction/expansion when the applied electric field exceeded a threshold value. Later, Bryce and Freeman [33] reported that such instabilities insignificantly enhanced the mixing in micro flows. Pimenta and Alves [34,35] later experimentally and numerically studied the electro-elastic instabilities of PAA solutions in both cross-slot and flow-focusing micro devices, and found that mixing efficiency was not improved significantly. Song et al. [36] experimentally and numerically studied the elastic instability in EOF of viscoelastic polyethylene oxide (PEO) solutions through T-shaped microchannels, and results demonstrated that the threshold electric field for onset of instability highly depended on the PEO concentration. Song et al. [37] later extended the work by experimentally investigating the fluid rheological effects on the elastic instability in EOF of six types of phosphate buffer-based aqueous solutions through T-shaped microchannels. They found that shear thinning effect of the fluid might account for the electro-elastic instabilities, however, the fluid with high elasticity alone did not have instability, which is inconsistent with the results of Pimenta [35]. The authors attribute the inconsistency to the neglect of microstructural effects (e.g., polymer-wall interaction and electric effect on molecular structure of polymer, etc.) of shear-thinning polymer solutions. However, this experimental result shows similarity to the work of Ko et al. [38], in which weakly shear-thinning, viscoelastic polyvinylpyrrolidone (PVP), and PEO solutions exhibited Newtonian-like EOF patterns, while shear-thinning and weakly elastic xanthan gum (XG) solution exhibited disturbance and vortices, suggesting that fluid elasticity alone has an insignificant impact on the steady-state EOF pattern. More recently, Sadek [39] experimentally investigated EOF of viscoelastic fluids through different microchannel configurations, including hyperbolicshaped contractions followed by an abrupt expansion, and abrupt contractions followed by a hyperbolic-shaped expansion, and EOF showed instabilities of elastic origin at very low Weissenberg numbers (*Wi*) (i.e., *Wi* < 0.01).

There is only limited literature on numerical studies of electro-elastic instabilities. Afonso et al. [40] numerically investigated the elastic instability of EOF through a crossslot geometry using the upper-converted Maxwell and the simplified Phan-Thien-Tanner models, and a direct flow transition from steady symmetric state to unsteady flow without crossing the steady asymmetric state at a critical *Wi* was observed. Pimenta and Alves [35] numerically investigated the electro-elastic instabilities in cross-slot and flow-focusing micro devices using OB model and Poisson-Boltzmann (PB) model. They found that strong shear-dominated flow within the EDL at the corners had a more significant contribution to the elastic instabilities than the extensionally dominated bulk flow. Song et al. [36] numerically investigated EOF of PEO solution through a T-shaped microchannel. Their model considered only the influence of PEO solution on the fluid viscosity, conductivity, and zeta potential. Due to the neglect of fluid elasticity effect in the mathematical model, only the electrokinetic flow phenomena of dilute PEO solution (i.e., ≤750 ppm) were captured.

Both experimental and numerical investigations in the EOF instabilities of viscoelastic fluid are limited, and the conditions proposed by various researchers for triggering the instabilities in the EOF of viscoelastic fluids show inconsistency and remain unclear. Inspired by the existing literature, in this work we numerically study EOF of viscoelastic fluids through a 10:1:10 contraction microchannel. The geometry, which consists of a constriction microchannel connecting two relatively big reservoirs on either end, is close to actual microfluidic device. The time-dependent OB model and PB model are adopted to describe the constitutive characteristics and the electrokinetic phenomenon, respectively. EOFs of PAA solutions with various weight concentrations under different applied electric fields are investigated. The effects of polymer concentration and applied electric field on the elastic instability are studied, and a map in polymer concentration-electric field space for predicting the onset of upstream vortices is formed. Both experimental and numerical investigations in the EOF instabilities of viscoelastic fluid are limited, and the conditions proposed by various researchers for triggering the instabilities in the EOF of viscoelastic fluids show inconsistency and remain unclear. Inspired by the existing literature, in this work we numerically study EOF of viscoelastic fluids through a 10:1:10 contraction microchannel. The geometry, which consists of a constriction microchannel connecting two relatively big reservoirs on either end, is close to actual microfluidic device. The time-dependent OB model and PB model are adopted to describe the constitutive characteristics and the electrokinetic phenomenon, respectively. EOFs of PAA solutions with various weight concentrations under different applied electric fields are investigated. The effects of polymer concentration and applied electric field on the elastic instability are studied, and a map in polymer concentration-electric field space for predicting the onset of upstream vortices is formed.

only the electrokinetic flow phenomena of dilute PEO solution (i.e., ≤750 ppm) were cap-

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 3 of 29

#### **2. Mathematical Model 2. Mathematical Model**

tured.

We consider incompressible monovalent binary electrolyte solution such as KCl with bulk concentration *c*<sup>0</sup> mixed with PAA polymer solution of concentration *c*p, which fills a microchannel of height *H*c, length *L*c, and width *W* connecting two identical reservoirs of height *H*<sup>r</sup> and length *L*<sup>r</sup> on either side. The solid walls of the constriction microchannel and the reservoirs are assumed to carry a constant negative zeta potential, *ξ*0. When dealing with non-Newtonian fluids, a constant zeta potential has been widely accepted [41]. Huang et al. [41] compared theoretical and experimental results of PEO solutions, and a constant zeta potential was proven for various PEO concentrations. Therefore, in the current study, we neglect the effect of the polymer concentration on the wall zeta potential. Two electrodes are placed at both ends of the reservoirs, and an external potential bias *U*<sup>0</sup> is applied between the inlet (Anode) and outlet (Cathode). Through the interaction between the externally applied electric field and net charges accumulated within the EDL in the vicinity of the charged walls, EOF flowing from the anode reservoir through the constriction microchannel towards the cathode reservoir is generated. The apparent electric field between the inlet and outlet is defined as *E*app = *U*0/(2*L*<sup>r</sup> + *L*c). In some applications, there are slit microchannels with width much larger than height [42,43]. For example, two-phase flow patterns were studied in a microchannel with 10-mm width and 50-µm height [42]. For microchannels with such geometries, the flow can be simplified to a 2D problem [44]. Therefore, in the current study we assume that the channel width is much larger than the channel height, and the flow can be simplified to a 2D problem as schematically shown in Figure 1. A Cartesian coordinate system with origin fixed at the center of the microchannel is adopted with the *x*-axis along the length direction and the *y*-axis in the height direction. We consider incompressible monovalent binary electrolyte solution such as KCl with bulk concentration mixed with PAA polymer solution of concentration , which fills a microchannel of height , length , and width *W* connecting two identical reservoirs of height and length on either side. The solid walls of the constriction microchannel and the reservoirs are assumed to carry a constant negative zeta potential, . When dealing with non-Newtonian fluids, a constant zeta potential has been widely accepted [41]. Huang et al. [41] compared theoretical and experimental results of PEO solutions, and a constant zeta potential was proven for various PEO concentrations. Therefore, in the current study, we neglect the effect of the polymer concentration on the wall zeta potential. Two electrodes are placed at both ends of the reservoirs, and an external potential bias is applied between the inlet (Anode) and outlet (Cathode). Through the interaction between the externally applied electric field and net charges accumulated within the EDL in the vicinity of the charged walls, EOF flowing from the anode reservoir through the constriction microchannel towards the cathode reservoir is generated. The apparent electric field between the inlet and outlet is defined as ୟ୮୮ = /(2୰ + ୡ). In some applications, there are slit microchannels with width much larger than height [42,43]. For example, two-phase flow patterns were studied in a microchannel with 10 mm width and 50-μm height [42]. For microchannels with such geometries, the flow can be simplified to a 2D problem [44]. Therefore, in the current study we assume that the channel width is much larger than the channel height, and the flow can be simplified to a 2D problem as schematically shown in Figure 1. A Cartesian coordinate system with origin fixed at the center of the microchannel is adopted with the *x*-axis along the length direction and the *y*-axis in the height direction.

**Figure 1.** Schematic diagram of a constriction microchannel connecting two reservoirs at both ends. The solid walls of reservoirs and the constriction channel are negatively charged, and an **Figure 1.** Schematic diagram of a constriction microchannel connecting two reservoirs at both ends. The solid walls of reservoirs and the constriction channel are negatively charged, and an electric field is imposed by applying a potential difference between anode and cathode positioned in two fluid reservoirs.

The induced EOF of viscoelastic fluid is governed by the continuity and Navier–Stokes equations:

$$
\nabla \cdot \mathfrak{u} = 0,\tag{1}
$$

$$
\rho \left( \frac{\partial \mathfrak{u}}{\partial t} + \mathfrak{u} \cdot \nabla \mathfrak{u} \right) = -\nabla p + \eta\_{\mathrm{s}} \nabla^{2} \mathfrak{u} + \nabla \cdot \mathfrak{r} - \rho\_{\mathrm{E}} \nabla \phi\_{\mathrm{Ext}\prime} \tag{2}
$$

where *u* is the velocity; *p* is the pressure; *ρ* is the fluid density; *η<sup>s</sup>* denotes the solvent dynamic viscosity; τ is the polymeric stress tensor accounting for the memory of the viscoelastic fluid; *t* represents time; *ρ*<sup>E</sup> and *φ*Ext represent, respectively, the volume charge density within the electrolyte solution and the externally applied electric potential. For different types of viscoelastic fluids, various constitutive models have been developed to relate the polymeric stress tensor *τ* and the deformation rate of the fluid, including WhiteMetzner model [14], which is commonly used for shear-thinning fluid; PTT model [17], which has good performance for prediction of viscosity at low shear rates; Giesekus model [19], which is suitable for concentrated polymer solutions; and OB model [16], which is suitable for dilute polymer solutions. Since OB model can properly fit the rheological behavior of aqueous PAA solutions [45], OB model is adopted in this study. In the OB model, the polymeric stress tensor, *τ*, is described as [16],

$$
\boldsymbol{\pi} = \frac{\eta\_{\rm P}}{\lambda} (\mathbf{c} - \mathbf{I}) \,\tag{3}
$$

where *η*<sup>p</sup> is the polymer dynamic viscosity; *λ* is the relaxation time of the polymer, which refers to the time it takes for polymer chains to return to equilibrium after being disturbed; *c* is the symmetric conformation tensor of the polymer molecules; and *I* is the identity matrix.

For the OB model, the conformation tensor *c* is governed by [16],

$$\frac{\partial \mathbf{c}}{\partial t} + \mathfrak{u} \cdot \nabla \mathbf{c} = \mathbf{c} \cdot \nabla \mathfrak{u} + (\nabla \mathfrak{u})^T \cdot \mathbf{c} - \frac{1}{\lambda} (\mathbf{c} - \mathbf{I}).\tag{4}$$

Typically, numerical simulation of viscoelastic flow is difficult to converge for high Weissenberg number problem [46,47]. Computations were found to break down at frustratingly low values of Weissenberg number (usually around *Wi* = 1; precise critical value also depends on the flow geometry) [48]. Therefore, the log-conformation tensor approach [47] is adopted. In the log-conformation tensor method, a new tensor (**Θ**) is defined as the natural logarithm of the conformation tensor,

$$\Theta = \ln(\mathcal{c}) = \mathbf{R}\ln(\Lambda)\mathbf{R},\tag{5}$$

where **Λ** is a diagonal matrix whose diagonal elements are the eigenvalues of *c*; and **R** is an orthogonal matrix with its columns being the eigenvalues of *c*. Equation (4) for the conformation tensor written in terms of **Θ** then becomes [46],

$$\frac{\partial \Theta}{\partial t} + \boldsymbol{\mu} \cdot \nabla \Theta = \boldsymbol{\Omega} \Theta - \boldsymbol{\Theta} \boldsymbol{\Omega} + 2 \mathbf{B} + \frac{1}{\lambda} \left( e^{\boldsymbol{\Theta}} - \boldsymbol{I} \right). \tag{6}$$

in the above, **Ω** and **B** are, respectively, the anti-symmetric matrix and the symmetric traceless matrix of the decomposition of the velocity gradient tensor ∇*u* [46].

Then, the conformation tensor *c* is recovered from **Θ**,

$$\mathfrak{c} = \exp(\mathfrak{G}).\tag{7}$$

The total electric potential, Ψ, is decomposed in two variables, Ψ = *φ*Ext + *ψ* [35], with *φ*Ext representing the potential originated from the externally applied electric potential while *ψ* being the potential arising from the charge of channel walls. In this study, the EDL thickness is on the order of nanometers (the calculation of EDL thickness will be shown

in the next paragraph), while the microchannel height is on the order of micrometers. Therefore, the Poisson–Boltzmann equation [49] is used to describe the potential, *ψ*:

$$\nabla \cdot (\varepsilon \nabla \psi) = \rho\_E = Fc\_0 \left( \exp\left(\frac{e\psi}{kT}\right) - \exp\left(-\frac{e\psi}{kT}\right) \right). \tag{8}$$

In the above, *<sup>F</sup>* is the Faraday's constant (i.e., 96,485.33289 <sup>C</sup>·mol−<sup>1</sup> ); *e* is the elementary charge (i.e., 1.6021766341×10−<sup>19</sup> C); *<sup>k</sup>* is Boltzmann's constant (i.e., 1.380649×10−<sup>23</sup> <sup>J</sup>·<sup>K</sup> −1 ); *T* is the absolute temperature of the fluid (i.e., 295 K); and *ε* represents the permittivity of the solution (i.e., 6.906266 <sup>×</sup> <sup>10</sup>−<sup>10</sup> <sup>F</sup>·m−<sup>1</sup> ). In this study, the bulk concentration *c*<sup>0</sup> is 0.01 mM; *z*<sup>1</sup> = 1 and *z*<sup>2</sup> = −1. In a biocompatible solution with pH of 7.4, the concentration of H<sup>+</sup> is 10−7.4 mol/L, and the concentration of OH<sup>−</sup> is 10−6.6 mol/L. The concentration of weak electrolyte is relatively low comparing with the background salt. Therefore, the weak electrolyte is not considered in the current study. The EDL thickness can be calculated by *λ<sup>D</sup>* = r *εkT eF*(*z* 2 1 *c*0+*z* 2 2 *c*0) , which is 95 nm. The potential *φ*Ext is governed by the following Laplace equation [50],

$$
\nabla^2 \phi\_{Ext} = 0.\tag{9}
$$

The boundary conditions are given as follows (Figure A4):

(1) At the Anode (edge AG in Figure A4): **n**·∇*u* = 0; *p* = 0; τ = **0**; *φ*Ext = *U*0; **n**·∇*ψ* = 0; **Θ** = **0**; where **n** denotes the normal unit vector on the surface.

(2) At the Cathode (edge FL in Figure A4): **n**·∇*u* = 0; *p* = 0; **n**·∇τ = 0; *φ*Ext = 0; **n**·∇*ψ* = 0; **n**·∇**Θ** = 0.

(3) On the reservoir walls (edges ABC, DEF, GHI, and JKL in Figure A4) and the microchannel walls (edges CD and IJ in Figure A4): *u* = **0**; **n**·∇*φ*Ext = 0; *ψ* = *ξ*0; **n**·∇**Θ** = 0; **n**·∇*p* is obtained from the momentum equation; the components of τ are linearly extrapolated.

The following initial conditions are specified within the domain: *u* = **0**; *p* = 0; τ = **0**; *φ*Ext = 0; *ψ* = **0**; **Θ** = 0.

Note that the electric potentials and the flow are only one-way coupling. The electric potentials *φ*Ext and *ψ* are in a steady state, and they are independent on the flow. However, the electric potentials affect the flow through the electrostatic force, which is the last term in Equation (2). For Newtonian fluid, the third term, ∇·τ, in the right-hand-side of Equation (2) is dropped, and the model includes Equations (1), (2), (8), and (9).

#### **3. Numerical Method and Code Validation**

The governing equations are numerically solved using the finite volume method by RheoTool (version 4.1, https://github.com/fppimenta/rheoTool, accessed on 1 June 2020), an open-source viscoelastic EOF solver [35] implemented in the open-source OpenFOAM platform. The details of the solver can be found from the work of Pimenta and Alves [34,50]. To numerically solve the coupled Equations (1), (2), (6), (8), and (9) along with the boundary and initial conditions, CUBISTA scheme [51] is used to discretize the convective terms in Equations (2) and (6). Central differences are used for the discretization of Laplacian and gradient terms. The time derivatives are discretized with three-time level explicit difference scheme [52], which is of the second order of accuracy. The exponential source term in Equation (8) is linearized using Taylor expansion up to the second term [53]. All of the terms in the momentum equation (i.e., Equation (2)), except the pressure gradient and the electric contribution, are discretized implicitly. A small time-step, ∆*t* = *λ*/10<sup>5</sup> , is used to ensure the accuracy. The well-known SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations-Consistent) algorithm [54] is used to resolve the velocity-pressure coupling. An inner-iteration loop is used to reduce the explicitness of the method and increase its accuracy and stability. The pressure field is computed by PCG (Preconditioned Conjugate Gradient) solver, of which the tolerance and maximum iteration are set to be <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>8</sup> and 800, respectively. The velocity field is computed by PBiCG (Preconditioned Biconjugate

Gradient) solver, of which the tolerance and the maximum iteration are set to be <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> and 1000, respectively. The computational steps of the solver are as follows [34]: be 1 × 10ିଵ and 1000, respectively. The computational steps of the solver are as follows [34]: Step 1. Initialize the fields ሼ, , , ୶୲, , ሽ and time (=0).

Step 1. Initialize the fields {*u*, *p*, τ, *φ*Ext, *ψ*, **Θ**}<sup>0</sup> and time (*t* = 0). Step 1.1. Compute steady state ୶୲ from Equation (9) and from Equation

Step 1.1. Compute steady state *φ*Ext from Equation (9) and *ψ* from Equation (8). (8).

Step 2. Enter the time loop (*t* = ∆*t*). Step 2. Enter the time loop ( = ∆).

Step 2.1. Enter the inner iteration loop (*i* = 0). Step 2.1. Enter the inner iteration loop (=0).

Step 2.1.1. Compute **Θ***<sup>i</sup>* and τ*<sup>i</sup>* by log-conformation method. Step 2.1.1. Compute and by log-conformation method.

Step 2.1.2. Compute estimated velocity field *u* ∗ *i* by solving the momentum equation. Step 2.1.2. Compute estimated velocity field ∗ by solving the momentum equation.

Step 2.1.3. ompute pressure field *p<sup>i</sup>* by enforcing the continuity equation. Step 2.1.3. Compute pressure field by enforcing the continuity equation.

Step 2.1.4. Correct the previously estimated velocity field using the correct pressure field. Step 2.1.4. Correct the previously estimated velocity field using the correct pressure field.

Step 2.1.5. Increase the inner iteration index (*i* = *i* + 1) and repeat the computation from Step 2.1.1, until the inner iteration criteria (i.e., maximum tolerance) is satisfied. Step 2.1.5. Increase the inner iteration index (=+1) and repeat the computation from Step 2.1.1, until the inner iteration criteria (i.e., maximum

Step 2.1.6. Set {*u*, *p*, τ, *φ*Ext, *ψ*, **Θ**}*<sup>t</sup>* = {*u<sup>i</sup>* , *p<sup>i</sup>* , τ*<sup>i</sup>* , *φ*Ext*<sup>i</sup>* , *ψ<sup>i</sup>* , **Θ***i*}. tolerance) is satisfied. Step 2.1.6. Set ሼ, , , ୶୲, , ሽ௧ = ൛, , , ୶୲ , , ൟ.

Step 2.2. Increase time, *t* = *t* + ∆*t*, and return to Step 2.1 until the simulation time is reached. Step 2.2. Increase time, = + ∆, and return to Step 2.1 until the simulation

Step 3. Stop the simulation and exit. time is reached.

Structural mesh is adopted to discretize the computational domain. 90◦ corners of the contraction channel (points I, J, C, and D in Figure 1) are smoothed by a fillet of 1 µm in radius to avoid sharp turns. The 90◦ corners of the reservoirs (points H, K, B, and E in Figure 1) are smoothed by a fillet of 2 µm in radius. To capture the EDL in the vicinity of the charged walls, a finer mesh is distributed near the charged reservoir and channel walls as shown in Figure 2. To reduce the number of mesh, we use a relatively low bulk concentration *c*<sup>0</sup> = 0.01 mM, and the EDL thickness is 95 nm in this study. In order to capture the details in the EDL and to guarantee the accuracy, the mesh size near the charged wall is set to be 10 nm so that there are 10 meshes within the EDL. There are 77,192 meshes in the whole geometry. A mesh independence study, described in the Appendix A, is performed to ensure the accuracy of the simulation. Step 3. Stop the simulation and exit. Structural mesh is adopted to discretize the computational domain. 90° corners of the contraction channel (points I, J, C, and D in Figure 1) are smoothed by a fillet of 1 μm in radius to avoid sharp turns. The 90° corners of the reservoirs (points H, K, B, and E in Figure 1) are smoothed by a fillet of 2 μm in radius. To capture the EDL in the vicinity of the charged walls, a finer mesh is distributed near the charged reservoir and channel walls as shown in Figure 2. To reduce the number of mesh, we use a relatively low bulk concentration = 0.01 mM, and the EDL thickness is 95 nm in this study. In order to capture the details in the EDL and to guarantee the accuracy, the mesh size near the charged wall is set to be 10 nm so that there are 10 meshes within the EDL. There are 77,192 meshes in the whole geometry. A mesh independence study, described in the Appendix A, is per-

formed to ensure the accuracy of the simulation.

**Figure 2.** Computational mesh used in the numerical simulations. Mesh of the whole geometry (**a**) and detailed view of the mesh at channel corner (**b**), at reservoir corner (**c**), and in the constriction microchannel (**d**). **Figure 2.** Computational mesh used in the numerical simulations. Mesh of the whole geometry (**a**) and detailed view of the mesh at channel corner (**b**), at reservoir corner (**c**), and in the constriction microchannel (**d**).

tration, ୮.

In this work, *η*<sup>p</sup> and *λ* for 100 ppm, 250 ppm, and 1000 ppm PAA-water solutions [55] are adopted to accomplish curve fitting as shown in Figure 3. The values of *η*<sup>p</sup> and *λ* were experimentally measured [56], and the slow retraction method was used to measure the relaxation time. The polymer dynamic viscosity can be expressed as *<sup>η</sup>*<sup>p</sup> <sup>=</sup> 2.22 <sup>×</sup> <sup>10</sup>−<sup>5</sup> ·*c*p*,* and the relaxation time can be expressed as *<sup>λ</sup>* <sup>=</sup> 3.69 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>+</sup> 3.94222 <sup>×</sup> <sup>10</sup>−<sup>5</sup> ·*c*<sup>p</sup> <sup>+</sup> 9.68889 <sup>×</sup> <sup>10</sup>−<sup>5</sup> ·*c*p 2 , where *c*<sup>p</sup> represents the weight concentration of PAA solution with the unit of ppm. The *η*<sup>p</sup> and *λ* for other *c*<sup>p</sup> studied in this work are estimated by the curve-fitting expressions. [55] are adopted to accomplish curve fitting as shown in Figure 3. The values of ୮ and were experimentally measured [56], and the slow retraction method was used to measure the relaxation time. The polymer dynamic viscosity can be expressed as ୮ = 2.22 × 10ିହ ∙ ୮*,* and the relaxation time can be expressed as = 3.69 × 10ିଷ + 3.94222 × 10ିହ ∙ ୮ + 9.68889 × 10ିହ ∙ ୮ <sup>ଶ</sup>, where ୮ represents the weight concentration of PAA solution with the unit of ppm. The ୮ and for other ୮ studied in this work are estimated by the curve-fitting expressions.

In this work, ୮ and for 100 ppm, 250 ppm, and 1000 ppm PAA-water solutions

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 7 of 29

**Figure 3.** (**a**) Polymer dynamic viscosity ୮ and (**b**) relaxation time as a function of the polyacrylamide (PAA) concen-**Figure 3.** (**a**) Polymer dynamic viscosity *η*p and (**b**) relaxation time *λ* as a function of the polyacrylamide (PAA) concentration, *c*p.

chowski velocity formula [57],

In a microfluidic channel with EDL thickness much smaller than the channel height, the EOF velocity of a Newtonian fluid can be approximated by the Helmholtz–Smolu-In a microfluidic channel with EDL thickness much smaller than the channel height, the EOF velocity of a Newtonian fluid can be approximated by the Helmholtz–Smoluchowski velocity formula [57],

$$
\mu\_0 = -\frac{\epsilon \mathfrak{E}\_0 E\_x}{\eta\_0},
\tag{10}
$$

 = − ఌకబாೣ ఎబ , (10) where ௫ is the actual local electric field in the main stream direction, and is the total viscosity of the fluid. To check the accuracy of our code, we simulate EOFs of both Newtonian and viscoelastic fluids in the same geometry with ୡ = 40 μm, ୡ = 200 μm, ୰ = 400 μm, and ୰ = 400 μm. Other parameters are set as = 60 V, = −0.11 V [58], and = 6.906266 × 10ିଵ F ∙ mିଵ. For Newtonian fluid, the total viscosity is set as = ୱ = 0.00322 kg/(m ∙ s). When the concentration of PAA solution is less than 2 ppm, the relaxation time is less than 0.1 ms [56], and the fluid can be approximately treated as Newtonian fluid. Therefore, for the OB model, parameters are set as ୱ = 0.00317 kg/(m ∙ s), ୮ = 0.00005 kg/(m ∙ s), = ୱ + ୮ = 0.00322 kg/(m ∙ s), and = 0.1 ms. Figure 4a depicts electric potential ୶୲(, 0) along the *x*-axis when ୟ୮୮ = 600 V/cm. The electric field in the *x*-direction, <sup>−</sup> డథు౮౪ డ௫ , in the constriction microchannel is 1820 V/cm, which is about 10 times of the electric field in the reservoirs. This is because of the 10:1:10 contraction geometry and current conservation. With the same electric conductivity, the electric field is inversely proportional to the cross-sectional area of the geometry. Note that the actual electric field within the constriction microchannel is about three times of the apparent electric field, ୟ୮୮, which does not consider the cross-sectional variation of the geometry. EOFs of both Newtonian fluid and viscoelastic fluid reach a steady state. Figure 4b shows the *x*component velocity profiles, (0, ), of the Newtonian fluid (solid line) and the OB model (circles). The velocity first rises rapidly within the thickness of EDL, then reaches a plateau in the cross section of the channel. When ௫ = 1820 V/cm, the calculated Helmholtz– Smoluchowski velocity is 4.29 mm/s, and the velocity at the center of the channel is 4.27 mm/s for both Newtonian and OB models. The relative difference between the approximated velocity and the simulated velocity is less than 0.5%. In addition, the result for OB model matches that of Newtonian fluid. Such consistency between Newtonian model and where *E<sup>x</sup>* is the actual local electric field in the main stream direction, and *η*<sup>0</sup> is the total viscosity of the fluid. To check the accuracy of our code, we simulate EOFs of both Newtonian and viscoelastic fluids in the same geometry with *H*<sup>c</sup> = 40 µm, *L*<sup>c</sup> = 200 µm, *H*<sup>r</sup> = 400 µm, and *L*<sup>r</sup> = 400 µm. Other parameters are set as *U*<sup>0</sup> = 60 V, *ξ*<sup>0</sup> = −0.11 V [58], and *<sup>ε</sup>* <sup>=</sup> 6.906266 <sup>×</sup> <sup>10</sup>−<sup>10</sup> <sup>F</sup>·m−<sup>1</sup> . For Newtonian fluid, the total viscosity is set as *η*<sup>0</sup> = *η*<sup>s</sup> = 0.00322 kg/(m·s). When the concentration of PAA solution is less than 2 ppm, the relaxation time is less than 0.1 ms [56], and the fluid can be approximately treated as Newtonian fluid. Therefore, for the OB model, parameters are set as *η*<sup>s</sup> = 0.00317 kg/(m·s), *η*<sup>p</sup> = 0.00005 kg/(m·s), *η*<sup>0</sup> = *η*<sup>s</sup> + *η*<sup>p</sup> = 0.00322 kg/(m·s), and *λ* = 0.1 ms. Figure 4a depicts electric potential *φ*Ext(*x*, 0) along the *x*-axis when *E*app = 600 V/cm. The electric field in the *x*-direction, − *∂φ*Ext *∂x* , in the constriction microchannel is 1820 V/cm, which is about 10 times of the electric field in the reservoirs. This is because of the 10:1:10 contraction geometry and current conservation. With the same electric conductivity, the electric field is inversely proportional to the cross-sectional area of the geometry. Note that the actual electric field within the constriction microchannel is about three times of the apparent electric field, *E*app, which does not consider the cross-sectional variation of the geometry. EOFs of both Newtonian fluid and viscoelastic fluid reach a steady state. Figure 4b shows the *x*-component velocity profiles, *u*(0, *y*), of the Newtonian fluid (solid line) and the OB model (circles). The velocity first rises rapidly within the thickness of EDL, then reaches a plateau in the cross section of the channel. When *<sup>E</sup><sup>x</sup>* <sup>=</sup> <sup>1820</sup> V/cm, the calculatedHelmholtz–Smoluchowski velocity is 4.29 mm/s, and the velocity at the center of the channel is 4.27 mm/s for both Newtonian and OB models. The relative difference between the approximated velocity and the simulated velocity is less than 0.5%. In addition, the result for OB model matches that of Newtonian fluid. Such consistency between Newtonian model and OB model is because the polymer dynamic viscosity *η*<sup>p</sup> is much smaller than the solvent dynamic viscosity *η*s, and the relaxation time of the polymer *λ* is also tiny. Under

find that the numerical result also agrees well with the analytical solution when the zeta potential is -110 mV. Therefore, the agreement of results attained from the OB model and Newtonian model, which are also validated by the Helmholtz–Smoluchowski ap-

OB model is because the polymer dynamic viscosity ୮ is much smaller than the solvent

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 8 of 29

the same as Newtonian fluid with the same total viscosity.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 8 of 29

the same as Newtonian fluid with the same total viscosity.

OB model is because the polymer dynamic viscosity ୮ is much smaller than the solvent dynamic viscosity ୱ, and the relaxation time of the polymer is also tiny. Under the considered condition, the elastic effect of the fluid is negligible and the OB fluid is almost

Afonso et al.[28] derived an analytical solution of viscoelastic EOF between two parallel plates based on the Debey–Hückel approximation, which is valid under the condition of low zeta potential (i.e., ൏ 25 mV). To further validate our code for OB model, EOF of viscoelastic fluid with ୱ = 0.001 kg/(m ∙ s), ୮ = 0.00222 kg/(m ∙ s), = ୱ + ୮ = 0.00322 kg/(m ∙ s), and = 8.6 ms in a straight 2D channel (with height of 40 μm) is studied. These rheology parameters are corresponding to those of 100 ppm PAA solution. is set as 10 V, while is chosen as −10 mV and −110 mV, respectively. Under the considered conditions, the flows are steady state due to relatively low electric field strength. Figure 5 depicts the *x*-component velocity profile at the center of the channel, and our numerical results (triangles) are in excellent agreement with the analytical result

the considered condition, the elastic effect of the fluid is negligible and the OB fluid is almost the same as Newtonian fluid with the same total viscosity. proximation, as well as the agreement of analytical solution of OB model and numerical results for EOF of viscoelastic fluid in a straight channel, validate our code. dynamic viscosity ୱ, and the relaxation time of the polymer is also tiny. Under the considered condition, the elastic effect of the fluid is negligible and the OB fluid is almost

**Figure 4.** (**a**) Electric potential distribution (blue dash line shows the relative position of the geometry) along the *x*-axis; (**b**) the *x*-component velocity at the center of the constriction microchannel, (0, ), for Newtonian model (solid line) and OB model (symbol). **Figure 4.** (**a**) Electric potential distribution (blue dash line shows the relative position of the geometry) along the *x*-axis; (**b**) the *x*-component velocity at the center of the constriction microchannel, *u*(0, *y*), for Newtonian model (solid line) and OB model (symbol). and Newtonian model, which are also validated by the Helmholtz–Smoluchowski approximation, as well as the agreement of analytical solution of OB model and numerical results for EOF of viscoelastic fluid in a straight channel, validate our code.

**Figure 5.** The *x*-component velocity profile of viscoelastic electroosmotic flow (EOF) between two parallel plates: (**a**) Zeta potential is −10 mV; (**b**) zeta potential is −110 mV. Analytical result of Afonso et al. [28] derived an analytical solution of viscoelastic EOF between two parallel plates based on the Debey–Hückel approximation, which is valid under the condition of low zeta potential (i.e., *ξ*<sup>0</sup> < 25 mV). To further validate our code for OB model, EOF of viscoelastic fluid with *η*<sup>s</sup> = 0.001 kg/(m·s), *η*<sup>p</sup> = 0.00222 kg/(m·s), *η*<sup>0</sup> = *η*<sup>s</sup> + *η*<sup>p</sup> = 0.00322 kg/(m·s), and *λ* = 8.6 ms in a straight 2D channel (with height of 40 µm) is studied. These rheology parameters are corresponding to those of 100 ppm PAA solution. *U*<sup>0</sup> is set as 10 V, while *ξ*<sup>0</sup> is chosen as −10 mV and −110 mV, respectively. Under the considered conditions, the flows are steady state due to relatively low electric field strength. Figure 5 depicts the *x*-component velocity profile at the center of the channel, and our numerical results (triangles) are in excellent agreement with the analytical result (line). Although the analytical solution is based on the Debey–Hückel approximation, we find that the numerical result also agrees well with the analytical solution when the zeta potential *ξ*<sup>0</sup> is −110 mV. Therefore, the agreement of results attained from the OB model and Newtonian model, which are also validated by the Helmholtz–Smoluchowski approximation, as well as the agreement of analytical solution of OB model and numerical results for EOF of viscoelastic fluid in a straight channel, validate our code. **Figure 4.** (**a**) Electric potential distribution (blue dash line shows the relative position of the geometry) along the *x*-axis; (**b**) the *x*-component velocity at the center of the constriction microchannel, (0, ), for Newtonian model (solid line) and OB model (symbol).

**Figure 5.** The *x*-component velocity profile of viscoelastic electroosmotic flow (EOF) between two parallel plates: (**a**) Zeta potential is −10 mV; (**b**) zeta potential is −110 mV. Analytical result of **Figure 5.** The *x*-component velocity profile of viscoelastic electroosmotic flow (EOF) between two parallel plates: (**a**) Zeta potential is −10 mV; (**b**) zeta potential is −110 mV. Analytical result of Afonso et al. [28] (solid line) and current numerical result (symbol). The analytical solution is described in Appendix A.

#### **4. Results and Discussion 4. Results and Discussion**  Newtonian fluid is investigated to provide the reference flow characteristic for the

scribed in Appendix A.

Newtonian fluid is investigated to provide the reference flow characteristic for the contraction geometry. For the PAA solution with different concentration *c*p, the applied apparent electric field *E*app is varied from low values to high (i.e., 100–600 V/cm). In this section, first, we describe the flow pattern of Newtonian fluid and the time-dependent flow patterns of PAA solutions. Then, the instabilities of PAA solutions with various *E*app and *c*<sup>p</sup> are discussed and a flow map is formed based on the investigated values of *E*app and *c*p. Finally, statistical results of cross-sectional average velocity are presented. contraction geometry. For the PAA solution with different concentration ୮, the applied apparent electric field ୟ୮୮ is varied from low values to high (i.e., 100–600 V/cm). In this section, first, we describe the flow pattern of Newtonian fluid and the time-dependent flow patterns of PAA solutions. Then, the instabilities of PAA solutions with various ୟ୮୮ and ୮ are discussed and a flow map is formed based on the investigated values of ୟ୮୮ and ୮. Finally, statistical results of cross-sectional average velocity are presented.

Afonso et al. [28] (solid line) and current numerical result (symbol). The analytical solution is de-

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 9 of 29

#### *4.1. Instability of PAA Solutions 4.1. Instability of PAA Solutions*

For Newtonian fluids with various total viscosities, the EOF reaches a steady state under all conditions of the applied electric field strengths. There is no vortex occurring in the reservoirs and the constriction microchannel. The streamlines of Newtonian fluid show excellent symmetry about the *x*-axis. Additionally, the magnitude of the velocity, U(*x*, *y*), is symmetric about the *y*-axis, U(*x*, *y*) = U(−*x*, *y*). For EOF of PAA solutions, when *E*app and *c*<sup>p</sup> are relatively low, the flow pattern is similar to that of Newtonian fluid, and the flow reaches a steady state without vortex. With increasing *E*app and *c*p, however, the viscoelastic flow becomes time dependent and significant instabilities are observed. Figure 6 depicts the streamlines at different times when *E*app = 100 V/cm and *c*<sup>p</sup> = 500 ppm. Figure 7 depicts the streamlines at different times when *E*app = 600 V/cm and *c*<sup>p</sup> = 150 ppm. Figure 8 shows the velocity magnitudes as a function of time at three different locations, namely, upstream of the constriction microchannel (−3*H*c, 0), center of the constriction microchannel (0, 0), and downstream of the constriction microchannel (3*H*c, 0). For the EOF of both *c*<sup>p</sup> = 150 ppm and *c*<sup>p</sup> = 500 ppm, we observe strong instabilities and upstream vortices. For Newtonian fluids with various total viscosities, the EOF reaches a steady state under all conditions of the applied electric field strengths. There is no vortex occurring in the reservoirs and the constriction microchannel. The streamlines of Newtonian fluid show excellent symmetry about the *x*-axis. Additionally, the magnitude of the velocity, U(, ), is symmetric about the *y*-axis, U(, ) = U(−, ). For EOF of PAA solutions, when ୟ୮୮ and ୮ are relatively low, the flow pattern is similar to that of Newtonian fluid, and the flow reaches a steady state without vortex. With increasing ୟ୮୮ and ୮, however, the viscoelastic flow becomes time dependent and significant instabilities are observed. Figure 6 depicts the streamlines at different times when ୟ୮୮ = 100 V/cm and ୮ = 500 ppm. Figure 7 depicts the streamlines at different times when ୟ୮୮ = 600 V/cm and ୮ = 150 ppm. Figure 8 shows the velocity magnitudes as a function of time at three different locations, namely, upstream of the constriction microchannel (−3ୡ, 0), center of the constriction microchannel (0,0), and downstream of the constriction microchannel (3ୡ, 0). For the EOF of both ୮ = 150 ppm and ୮ = 500 ppm, we observe strong instabilities and upstream vortices.

**Figure 6.** Instability of EOF with ୮ = 500 ppm and ୟ୮୮ = 100 V/cm. Streamlines at different times: (**a**) 1.71 s, (**b**) 1.75 s, (**c**) 1.79 s, (**d**) 1.83 s, (**e**) 1.87 s, and (**f**) 1.91 s. The color bar represents the elastic normal stress ௫௫. (Supplementary materials **Video S1** ). **Figure 6.** Instability of EOF with *<sup>c</sup>*<sup>p</sup> <sup>=</sup> <sup>500</sup> ppm and *<sup>E</sup>*app <sup>=</sup> <sup>100</sup> V/cm. Streamlines at different times: (**a**) 1.71 s, (**b**) 1.75 s,(**c**) 1.79 s, (**d**) 1.83 s, (**e**) 1.87 s, and (**f**) 1.91 s. The color bar represents the elastic normal stress *<sup>τ</sup>xx*. (Supplementary Materials Video S1).

**Figure 7.** Instability of EOF with ୮ = 150 ppm and ୟ୮୮ = 600 V/cm. Streamlines at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, and (**f**) 1.80 s. The color bar represents the elastic normal stress ௫௫. (Supplementary materials **Video S2**). **Figure 7.** Instability of EOF with *c*p = 150 ppm and *E*app = 600 V/cm. Streamlines at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, and (**f**) 1.80 s. The color bar represents the elastic normal stress *τxx*. (Supplementary Materials Video S2). **Figure 7.** Instability of EOF with ୮ = 150 ppm and ୟ୮୮ = 600 V/cm. Streamlines at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, and (**f**) 1.80 s. The color bar represents the elastic normal stress ௫௫. (Supplementary materials **Video S2**).

**Figure 8.** Velocity magnitudes at three different locations ((−3ୡ, 0), (0,0), (3ୡ, 0)). (**a**) ୮ = 150 ppm and ୟ୮୮ = 600 V/cm, (**b**) ୮ = 500 ppm and ୟ୮୮ = 100 V/cm. **Figure 8.** Velocity magnitudes at three different locations ((−3ୡ, 0), (0,0), (3ୡ, 0)). (**a**) ୮ = 150 ppm and ୟ୮୮ = 600 V/cm, (**b**) ୮ = 500 ppm and ୟ୮୮ = 100 V/cm. **Figure 8.** Velocity magnitudes at three different locations ((−3*H*c, 0), (0, 0), (3*H*c, 0)). (**a**) *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm, (**b**) *c*p = 500 ppm and *E*app = 100 V/cm.

Figures 6 and 7 show that the viscoelastic EOF is time-dependent. The streamlines in the left inlet reservoir far away from the solid walls (AB and GH in Figure 1) and near the entrance of the constriction microchannel show significant fluctuation and also become asymmetric about the *x*-axis. However, the streamlines near the solid walls of both reservoirs and in the outlet reservoir show insignificant change with time. Within 0.1 s, vortices continuously form and disappear within the inlet reservoir right before the entrance of the constriction microchannel. In Figure 6a, we observe significant curvature at the streamlines of the EOF upstream of the constriction microchannel. Then, the curvature of the streamlines further develops into a pair of vortices as shown in Figure 6b. Such vortices are also time-dependent. The size and shape of the vortices show notable differences at different times. After growing to the maximum size, the vortices start to shrink until the vortices break and disappear, as shown in Figure 6,d. Next, the vortices in the EOF Figures 6 and 7 show that the viscoelastic EOF is time-dependent. The streamlines in the left inlet reservoir far away from the solid walls (AB and GH in Figure 1) and near the entrance of the constriction microchannel show significant fluctuation and also become asymmetric about the *x*-axis. However, the streamlines near the solid walls of both reservoirs and in the outlet reservoir show insignificant change with time. Within 0.1 s, vortices continuously form and disappear within the inlet reservoir right before the entrance of the constriction microchannel. In Figure 6a, we observe significant curvature at the streamlines of the EOF upstream of the constriction microchannel. Then, the curvature of the streamlines further develops into a pair of vortices as shown in Figure 6b. Such vortices are also time-dependent. The size and shape of the vortices show notable differences at different times. After growing to the maximum size, the vortices start to shrink until the vortices break and disappear, as shown in Figure 6,d. Next, the vortices in the EOF Figures 6 and 7 show that the viscoelastic EOF is time-dependent. The streamlines in the left inlet reservoir far away from the solid walls (AB and GH in Figure 1) and near the entrance of the constriction microchannel show significant fluctuation and also become asymmetric about the *x*-axis. However, the streamlines near the solid walls of both reservoirs and in the outlet reservoir show insignificant change with time. Within 0.1 s, vortices continuously form and disappear within the inlet reservoir right before the entrance of the constriction microchannel. In Figure 6a, we observe significant curvature at the streamlines of the EOF upstream of the constriction microchannel. Then, the curvature of the streamlines further develops into a pair of vortices as shown in Figure 6b. Such vortices are also time-dependent. The size and shape of the vortices show notable differences at different times. After growing to the maximum size, the vortices start to shrink until the vortices break and disappear, as shown in Figure 6,d. Next, the vortices in the EOF keep forming and breaking repeatedly as shown in Figure 6e,f. Comparing Figure 6b,c, the central locations of the vortices are both spatially and temporally dependent. In Figure 6b,

keep forming and breaking repeatedly as shown in Figure 6e,f. Comparing Figure 6b,c,

keep forming and breaking repeatedly as shown in Figure 6e,f. Comparing Figure 6b,c,

6b, the direction of the circulation is marked by red curved lines. The pair of vortices are

the direction of the circulation is marked by red curved lines. The pair of vortices are in opposite directions and form a stagnant region right before the entrance of the constriction microchannel. Therefore, we call the induced vortices as entrance-centerline vortices.

For *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm, the width and length of the vortices are nearly the same as the height of the constriction microchannel (*H*c), while for *c*<sup>p</sup> = 500 ppm and *E*app = 100 V/cm, the width and the length of the vortices are about 2*H*c. In the similar geometry, however, Ko [38] did not observe vortices in their experiments with 200 ppm PAA solution under *E*app ranging from 75 V/cm to 200 V/cm. The elastic instability increases with increasing polymer concentration and the applied electric field. Our later results discussed in the following section show that for *c*<sup>p</sup> = 200 ppm, the vortices occur when the applied electric field exceeds the threshold value of 300 V/cm. Therefore, our numerical results qualitatively agree with the experimental observation of Ko [38] under their experimental condition. Table A1 in the Appendix A summarizes the EOF instabilities from the literature. In Sadek's [39] experimental study, small vortices at the entrance and large upstream circulation flows were observed. For highly concentrated polymer solution, downstream circulation flows were observed at a critical voltage. The upstream vortices found in this study are distinct from the small vortices and large circulation flows found in Sadek's study [39] in terms of location. Note that the geometry in our study differs significantly from the experimental study of Sadek [39], and we do not observe large circulating flows near the reservoir corners and channel lips.

For Newtonian fluid, the EOF reaches a steady state, and the velocity magnitudes are symmetric about the *y*-axis. Therefore, for Newtonian fluid, we have U(−3*Hc*, 0) = U(3*H*c, 0). However, as shown in Figure 8, the velocity magnitudes at three points (−3*H*c, 0), (0, 0), and (3*H*c, 0) fluctuate around certain values and velocity magnitudes do not show symmetry about the *y*-axis. For *c*<sup>p</sup> = 500 ppm and *E*app = 100 V/cm, as shown in Figure 8b, the time-averaged velocity at the channel center is 0.201 mm/s with a standard deviation of 0.013 mm/s. The time-averaged velocity at downstream of the constriction microchannel is 0.108 mm/s, which is about 2 times of that at upstream of the constriction microchannel (i.e., 0.047 mm/s). For *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm, the time-averaged velocities at the upstream, center, and downstream of the constriction microchannel are, respectively, 0.19 mm/s, 2.82 mm/s, and 2.05 mm/s. The ratio of the downstream velocity to upstream velocity is about 10 times. In contrast to Newtonian EOF, the flow velocity of viscoelastic fluid at the downstream is significantly higher than that at the upstream, which has also been experimentally observed in Ko's [38] experiments, where a fluid jet after the constriction microchannel was observed and the ratio of the velocity at the downstream centerline to that at upstream of the constriction microchannel varies between 1 and 2 under *E*app ranging from 75 V/cm to 200 V/cm and *c*<sup>p</sup> = 200 ppm.

Figure 8 also shows that EOF of *c*<sup>p</sup> = 500 ppm and *E*app = 100 V/cm presents stronger instabilities than that of *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm. Comparing Figures 6 and 7, the streamlines show stronger fluctuation and larger upstream vortices for solution with relatively high polymer concentration. Such trend suggests that although the increase of both *E*app and *c*<sup>p</sup> can enhance the instabilities of the viscoelastic EOF, the polymer concentration, *c*p, affects the instabilities of the EOF more significantly, which will be further discussed in next section.

Figures 6 and 7 also show the spatial distribution of elastic normal stress *τxx* with the color bar representing its magnitude. To clearly reveal it, Figure 9 depicts the spatial distribution of *τxx* in the whole geometry for *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm at *t* = 1.78 s. Within the two reservoirs, the elastic normal stress is nearly zero at location far away from the constriction microchannel. However, significant elastic normal stress is induced near the entrance of the constriction microchannel and near the downstream lips. Due to the contraction geometry, the electric field within the constriction microchannel is about 10 times of that within the inlet reservoir as shown in Figure 4a, and the flow velocity in the microchannel is significantly higher than that in the reservoir. For example, Figure 8 shows that the ratio of the time-averaged velocity within the microchannel to that in the

inlet reservoir, *U*(0, 0)/U(−3*H*c, 0), is 4.28 for *c*<sup>p</sup> = 500 ppm and *E*app = 100 V/cm and 10.79 for *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm. Near the entrance of the microchannel, the high velocity gradient results in a strong extension of polymer molecules, and consequently induces significant elastic normal stress. Therefore, *τxx* experiences a rapid increase near the entrance of the constriction microchannel. At the exit of the constriction microchannel, similarly, a significant increase of *τxx* is induced at the exit lips. Figure 10a,b depict the streamlines in the constriction microchannel and the color bar represents the velocity magnitude, U, for *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm and Newtonian fluid with the same total viscosity and *E*app at *t* = 1.78 s, respectively. For the viscoelastic fluid, at both the entrance and exit of the constriction microchannel, as shown by the dashed circles in Figure 10a, velocity becomes spatially dependent along the *y*-axis. Velocity near the walls of the constriction microchannel is significantly higher than that at the centerline of the microchannel, and a local maximum occurs near the inlet/outlet corners of the constriction microchannel. However, in the EOF of Newtonian fluid, as shown in Figure 10b, at both the entrance and exit of the constriction microchannel, the velocity magnitude is more evenly distributed in the cross section of the constriction microchannel. Figure 10c depicts the velocity magnitude profile at the entrance (2*x*/*H*<sup>c</sup> = −5) and exit (2*x*/*H*<sup>c</sup> = 5) of the constriction microchannel. For Newtonian fluid, the velocity magnitude profile is identical at 2*x*/*H*<sup>c</sup> = ±5 and is symmetric about the *x*-axis. The ratio of the maximum velocity magnitude near the channel walls to that at the centerline is 1.6. However, for PAA solution, due to the elastic instability, the velocity magnitude profile is asymmetric about the *x*-axis at 2*x*/*H*<sup>c</sup> = ±5. In addition, the ratios of the maximum velocity magnitude near the channel walls to that at the centerline are 9.7 and 4 at 2*x*/*H*<sup>c</sup> = −5 and 2*x*/*H*<sup>c</sup> = 5, respectively, which are much higher than that of the Newtonian fluid. For Newtonian fluid, the velocity profile is symmetric about the centerline of the microchannel (i.e., *y* = 0). However, Figure 10c shows that the local maximum velocity near the top channel wall differs from that near the bottom channel wall, and the velocity profile is asymmetric about *y* = 0. microchannel, similarly, a significant increase of ௫௫ is induced at the exit lips. Figure 10a,b depict the streamlines in the constriction microchannel and the color bar represents the velocity magnitude, U, for ୮ = 150 ppm and ୟ୮୮ = 600 V/cm and Newtonian fluid with the same total viscosity and ୟ୮୮ at = 1.78 s, respectively. For the viscoelastic fluid, at both the entrance and exit of the constriction microchannel, as shown by the dashed circles in Figure 10a, velocity becomes spatially dependent along the *y*-axis. Velocity near the walls of the constriction microchannel is significantly higher than that at the centerline of the microchannel, and a local maximum occurs near the inlet/outlet corners of the constriction microchannel. However, in the EOF of Newtonian fluid, as shown in Figure 10b, at both the entrance and exit of the constriction microchannel, the velocity magnitude is more evenly distributed in the cross section of the constriction microchannel. Figure 10c depicts the velocity magnitude profile at the entrance (2 ⁄ <sup>ୡ</sup> = −5) and exit (2 ⁄ <sup>ୡ</sup> = 5) of the constriction microchannel. For Newtonian fluid, the velocity magnitude profile is identical at 2 ⁄ <sup>ୡ</sup> = ±5 and is symmetric about the *x*-axis. The ratio of the maximum velocity magnitude near the channel walls to that at the centerline is 1.6. However, for PAA solution, due to the elastic instability, the velocity magnitude profile is asymmetric about the *x*-axis at 2 ⁄ <sup>ୡ</sup> = ±5. In addition, the ratios of the maximum velocity magnitude near the channel walls to that at the centerline are 9.7 and 4 at 2 ⁄ <sup>ୡ</sup> = −5 and 2 ⁄ <sup>ୡ</sup> = 5, respectively, which are much higher than that of the Newtonian fluid. For Newtonian fluid, the velocity profile is symmetric about the centerline of the microchannel (i.e., *y*=0). However, Figure 10c shows that the local maximum velocity near the top channel wall differs from that near the bottom channel wall, and the velocity profile is asymmetric about =0.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 12 of 29

ity in the microchannel is significantly higher than that in the reservoir. For example, Figure 8 shows that the ratio of the time-averaged velocity within the microchannel to that in the inlet reservoir, (0,0)/U(−3ୡ, 0), is 4.28 for ୮ = 500 ppm and ୟ୮୮ = 100 V/cm and 10.79 for ୮ = 150 ppm and ୟ୮୮ = 600 V/cm. Near the entrance of the microchannel, the high velocity gradient results in a strong extension of polymer molecules, and consequently induces significant elastic normal stress. Therefore, ௫௫ experiences a rapid increase near the entrance of the constriction microchannel. At the exit of the constriction

**Figure 9.** Spatial distribution of the elastic normal stress ௫௫ for ୮ = 150 ppm and ୟ୮୮ = 600 V/cm at = 1.78 s. **Figure 9.** Spatial distribution of the elastic normal stress *τxx* for *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm at *t* = 1.78 s.

As PAA solution flows from the microchannel into the outlet reservoir, fluid velocity first decreases when fluid exits the microchannel and then increases in the outlet reservoir, as shown by the region marked with a circle in Figure 10a and by the velocity magnitude as a function of *x* at =0 in Figure 10d. The two dashed lines in Figure 10d represent the As PAA solution flows from the microchannel into the outlet reservoir, fluid velocity first decreases when fluid exits the microchannel and then increases in the outlet reservoir, as shown by the region marked with a circle in Figure 10a and by the velocity magnitude as a function of *x* at *y* = 0 in Figure 10d. The two dashed lines in Figure 10d represent the entrance and exit of the constriction microchannel. EOF of Newtonian fluid within the constriction is a plateau, and its velocity magnitude within the constriction is

entrance and exit of the constriction microchannel. EOF of Newtonian fluid within the

much higher than those at both reservoirs, and this is because the electric field within the constriction microchannel is significantly higher than that in the reservoirs. However, the velocity of PAA solution becomes spatially dependent within the constriction, and a local maximum occurs before the exit and a local minimum occurs at the exit of the constriction microchannel. In addition, a local maximum occurs at the downstream outlet reservoir. Figure 10d also clearly shows that the velocity in the downstream outlet reservoir is significantly higher than that at the upstream inlet reservoir. For example, U(2*x*/*H*<sup>c</sup> = 10.0)/U(2*x*/*H*<sup>c</sup> = −10.0) = 3.37. The unexpected velocity decrease at the microchannel exit and velocity increase at the downstream outlet reservoir do not occur in Newtonian fluid as shown in Figure 10b,d. Such a phenomenon is probably because of the extrudate swell effect of polymers [59]. At the exit of the constriction microchannel, curved streamlines tilting toward the walls of the constriction microchannel are observed in viscoelastic fluid, suggesting that fluid tends to flow toward the charged walls of the microchannel. Such lateral velocity component results in the velocity's increase near the microchannel walls and velocity's decrease near the centerline at the exit of the constriction microchannel. In addition, the significant increase of *τxx* near the downstream lips observed in Figure 9 can also be attributed to the extrudate swell effect of polymers when polymer exits from the constriction microchannel to larger outlet reservoir. than those at both reservoirs, and this is because the electric field within the constriction microchannel is significantly higher than that in the reservoirs. However, the velocity of PAA solution becomes spatially dependent within the constriction, and a local maximum occurs before the exit and a local minimum occurs at the exit of the constriction microchannel. In addition, a local maximum occurs at the downstream outlet reservoir. Figure 10d also clearly shows that the velocity in the downstream outlet reservoir is significantly higher than that at the upstream inlet reservoir. For example, U(2 ⁄ <sup>ୡ</sup> = 10.0)⁄U(2 ⁄ <sup>ୡ</sup> = −10.0) = 3.37. The unexpected velocity decrease at the microchannel exit and velocity increase at the downstream outlet reservoir do not occur in Newtonian fluid as shown in Figure 10b,d. Such a phenomenon is probably because of the extrudate swell effect of polymers [59]. At the exit of the constriction microchannel, curved streamlines tilting toward the walls of the constriction microchannel are observed in viscoelastic fluid, suggesting that fluid tends to flow toward the charged walls of the microchannel. Such lateral velocity component results in the velocity's increase near the microchannel walls and velocity's decrease near the centerline at the exit of the constriction microchannel. In addition, the significant increase of ௫௫ near the downstream lips observed in Figure 9 can also be attributed to the extrudate swell effect of polymers when polymer exits from the constriction microchannel to larger outlet reservoir.

constriction is a plateau, and its velocity magnitude within the constriction is much higher

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 13 of 29

**Figure 10.** Streamlines and velocity magnitude for ୮ = 150 ppm and ୟ୮୮ = 600 V/cm and Newtonian fluid at = 1.78 s: (**a**) 150 ppm PAA solution, (**b**) Newtonian fluid with same total viscosity as 150 ppm PAA solution, (**c**) velocity magnitude profiles at 2 ⁄ <sup>ୡ</sup> = ±5, (**d**) velocity magnitudes profiles at =0 (The blue dash lines show the position of the contraction microchannel). The color bar represents the velocity magnitude U. **Figure 10.** Streamlines and velocity magnitude for *c*<sup>p</sup> = 150 ppm and *E*app = 600 V/cm and Newtonian fluid at *t* = 1.78 s: (**a**) 150 ppm PAA solution, (**b**) Newtonian fluid with same total viscosity as 150 ppm PAA solution, (**c**) velocity magnitude profiles at 2*x*/*H*<sup>c</sup> = ±5, (**d**) velocity magnitudes profiles at *y* = 0 (The blue dash lines show the position of the contraction microchannel). The color bar represents the velocity magnitude U.

The generated elastic normal stress ௫௫ was typically used to explain the formation of vortices in pressure-driven viscoelastic flows within curved geometries [60–62]. It has been reported that the development of the polymeric elastic stresses is caused by the flowinduced changes of the polymer conformation in the solution. Such changes of the polymer conformation are strain-dependent, anisotropic, and dependent on the flow. The extra elastic stresses are nonlinear under shear and can alter the flow behavior. At low Reynolds numbers where inertia is negligible, when the elastic normal stress exceeds by a certain amount the local shear stress, the flow transits from stable to unstable, and the vortices form at upstream of the constriction microchannel. Such elastic instabilities are often The generated elastic normal stress *τxx* was typically used to explain the formation of vortices in pressure-driven viscoelastic flows within curved geometries [60–62]. It has been reported that the development of the polymeric elastic stresses is caused by the flow-induced changes of the polymer conformation in the solution. Such changes of the polymer conformation are strain-dependent, anisotropic, and dependent on the flow. The extra elastic stresses are nonlinear under shear and can alter the flow behavior. At low Reynolds numbers where inertia is negligible, when the elastic normal stress exceeds by a certain amount the local shear stress, the flow transits from stable to unstable, and the vortices form at upstream of the constriction microchannel. Such elastic instabilities are often observed in flows with sufficient curvature [63–65], and some argue that curvature is necessary for infinitesimal perturbations to be amplified by the normal stress imbalances

in the viscoelastic flows [66]. However, other theoretical studies reported that viscoelastic flows also showed a nonlinear instability in parallel shear flows, such as in viscoelastic flows within straight pipes at low Reynolds numbers [67]. Although the formation of the upstream vortices in the viscoelastic EOF shares the same mechanism as the pressure-driven flow, the locations of the vortices found in this study are distinct from the typical lip and corner vortices occurring in pressure-driven viscoelastic flows. This is probably because of the different velocity profiles in the pressure-driven flow and the EOF. In pressure-driven flow, the velocity is zero at solid walls and increases to a maximum at the centerline of the geometry. However, the EOF velocity profile is nearly a plug flow as shown in Figure 4b. The velocity increases from zero to a plateau within the EDL thickness, which is only on the order of a few nanometers. For the pressure-driven flow, the highest velocity is at the centerline of the geometry and the velocity near the wall is relatively low, resulting in the stagnant region near the solid boundaries (lips and corners). However, EOF velocity in the vicinity of the charged wall is almost the same as that in the channel centerline. For the extensional flow of viscoelastic fluids, the stretched polymer molecules lead to large elastic stresses, which significantly depend on the geometry and velocity profile. The induced elastic stresses render the primary flow unstable and cause an irregular secondary flow. The flow subsequently acts back on the polymer molecules and stretches them further, causing a strong disturbance of the EOF and yielding a time-dependent EOF.

## *4.2. Elastic Instabilities under Various Eapp and c<sup>p</sup>*

In order to study the effects of *E*app and *c*<sup>p</sup> on the instabilities of viscoelastic EOF, *c*<sup>p</sup> is varied from 100 ppm to 500 ppm and *E*app is varied from 100 V/cm to 600 V/cm. Flow patterns under different conditions of *E*app and *c*<sup>p</sup> are shown in Figures 11–14. At certain *E*app (*c*p), EOF becomes more unstable with the increase of *c*<sup>p</sup> (*E*app). Figure 11 shows the streamlines for different PAA concentrations under *E*app = 600 V/cm, and Figure 12 shows the streamlines within the constriction microchannel with the color bar representing pressure for Newtonian fluid and *τxx* for PAA solutions. As shown in Figures 11 and 12, when *c*<sup>p</sup> increases, the polymeric stress *τxx* at the entrance of the constriction microchannel increases rapidly, resulting in the fluctuation of the streamlines at upstream of the microchannel. When *c*<sup>p</sup> is relatively low (100 ppm), EOF of viscoelastic fluid is similar to that of Newtonian fluid, and the flow is in a steady state. With an increase in the PAA concentration up to 150 ppm, significant curvature of the centerline streamlines is observed, and the streamlines become asymmetric about the *x*-axis and *y*-axis. As *c*<sup>p</sup> continuously increases up to 200 ppm, a pair of upstream vortices in opposite flow directions are induced at upstream of the constriction microchannel, forming a stagnant region as shown in Figure 11d. The width and length of the pair of vortices are about 1.6 time of the constriction microchannel height (i.e., 1.6*H*c). Such vortices are found to grow significantly in size with increasing *c*p, which is in qualitative agreement with the experimental observations of Ko [38]. Within the constriction microchannel, as shown in Figure 12d, nearly 1/4 of the microchannel length (i.e., <sup>1</sup> 4 *Lc*) from the entrance shows a significant increase of *τxx*. Near the downstream lips of the microchannel, a local maximum of the polymeric stress *τxx* is observed. When *c*<sup>p</sup> increases to 250 ppm, the fluctuation of the streamlines and the size of the vortices grow dramatically as shown in Figure 11e, in which the width and length of the vortices are about 2.9 times of the microchannel height (i.e., 2.9*H*c). The region with significant value of *τxx* is near 1/3 of the microchannel length (i.e., <sup>1</sup> 3 *Lc*), as shown in Figure 12e. When *c*<sup>p</sup> further increases to 500 ppm, as shown in Figure 11f, the vortices grow into 4.4 times of the constriction microchannel height (i.e., 4.4*H*c). More than half of the microchannel length shows a significant increase in *τxx*. Furthermore, a small vortex is induced near the downstream lip of the microchannel, which is also reported in experimental studies of Ko [38].

which is also reported in experimental studies of Ko [38].

which is also reported in experimental studies of Ko [38].

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 15 of 29

Furthermore, a small vortex is induced near the downstream lip of the microchannel,

Furthermore, a small vortex is induced near the downstream lip of the microchannel,

**Figure 11.** Streamlines of Newtonian fluid and PAA solutions with different concentrations under ୟ୮୮ = 600 V/cm at 1.70 s: (**a**) Newtonian fluid, (**b**) ୮ = 100 ppm, (**c**) ୮ = 150 ppm, (**d**) ୮ = 200 ppm, (**e**) ୮ = 250 ppm, and (**f**) ୮ = 500 ppm. The color bar represents the elastic normal stress ௫௫. **Figure 11.** Streamlines of Newtonian fluid and PAA solutions with different concentrations under *E*app = 600 V/cm at 1.70 s: (**a**) Newtonian fluid, (**b**) *c*p = 100 ppm, (**c**) *c*p = 150 ppm, (**d**) *c*p = 200 ppm, (**e**) *c*p = 250 ppm, and (**f**) *c*p = 500 ppm. The color bar represents the elastic normal stress *τxx*. **Figure 11.** Streamlines of Newtonian fluid and PAA solutions with different concentrations under ୟ୮୮ = 600 V/cm at 1.70 s: (**a**) Newtonian fluid, (**b**) ୮ = 100 ppm, (**c**) ୮ = 150 ppm, (**d**) ୮ = 200 ppm, (**e**) ୮ = 250 ppm, and (**f**) ୮ = 500 ppm. The color bar represents the elastic normal stress ௫௫.

**Figure 12.** Streamlines in microchannel of Newtonian fluid and PAA solutions with different concentrations under ୟ୮୮ = 600 V/cm at 1.7 s: (**a**) Newtonian fluid, (**b**) ୮ = 100 ppm, (**c**) ୮ = 150 ppm, (**d**) ୮ = 200 ppm, (**e**) ୮ = 250 ppm, and (**f**) ୮ = 500 ppm. The color bar represents the elastic normal stress ௫௫. **Figure 12.** Streamlines in microchannel of Newtonian fluid and PAA solutions with different concentrations under ୟ୮୮ = 600 V/cm at 1.7 s: (**a**) Newtonian fluid, (**b**) ୮ = 100 ppm, (**c**) ୮ = 150 ppm, (**d**) ୮ = 200 ppm, (**e**) ୮ = 250 ppm, and (**f**) ୮ = 500 ppm. The color bar represents the elastic normal stress ௫௫. **Figure 12.** Streamlines in microchannel of Newtonian fluid and PAA solutions with different concentrations under *<sup>E</sup>*app <sup>=</sup> <sup>600</sup> V/cm at 1.7 s: (**a**) Newtonian fluid, (**b**) *<sup>c</sup>*<sup>p</sup> <sup>=</sup> <sup>100</sup> ppm, (**c**) *<sup>c</sup>*<sup>p</sup> <sup>=</sup> <sup>150</sup> ppm, (**d**) *<sup>c</sup>*<sup>p</sup> <sup>=</sup> <sup>200</sup> ppm, (**e**) *<sup>c</sup>*<sup>p</sup> <sup>=</sup> <sup>250</sup> ppm,and (**f**) *<sup>c</sup>*<sup>p</sup> <sup>=</sup> 500 ppm. The color bar represents the elastic normal stress *<sup>τ</sup>xx*.

Figures 13 and 14 show the streamlines for ୮ = 150 ppm when ୟ୮୮ is varied from 100 V/cm to 600 V/cm with the color bar representing the magnitude of ௫௫. At a relatively low electric field such as ୟ୮୮ = 100 V/cm, EOF of PAA solution is similar to the Newtonian fluid, and the flow is in a steady state and symmetric about channel centerline. In addition, the induced polymeric stress ௫௫ in the constriction microchannel is relatively small. When ୟ୮୮ increases up to 400 V/cm, centerline streamlines start to show notable fluctuation and become asymmetric about the *x*-axis. At the entrance of the microchannel, a slight increase of ௫௫ is observed, however, significant increase of ௫௫ is observed near the microchannel walls and the downstream lips, as shown in Figure 14d. When ୟ୮୮ increases to 500 V/cm, a pair of upstream vortices are induced at upstream of Figures 13 and 14 show the streamlines for ୮ = 150 ppm when ୟ୮୮ is varied from 100 V/cm to 600 V/cm with the color bar representing the magnitude of ௫௫. At a relatively low electric field such as ୟ୮୮ = 100 V/cm, EOF of PAA solution is similar to the Newtonian fluid, and the flow is in a steady state and symmetric about channel centerline. In addition, the induced polymeric stress ௫௫ in the constriction microchannel is relatively small. When ୟ୮୮ increases up to 400 V/cm, centerline streamlines start to show notable fluctuation and become asymmetric about the *x*-axis. At the entrance of the microchannel, a slight increase of ௫௫ is observed, however, significant increase of ௫௫ is observed near the microchannel walls and the downstream lips, as shown in Figure 14d. When ୟ୮୮ increases to 500 V/cm, a pair of upstream vortices are induced at upstream of Figures 13 and 14 show the streamlines for *c*<sup>p</sup> = 150 ppm when *E*app is varied from 100 V/cm to 600 V/cm with the color bar representing the magnitude of *τxx*. At a relatively low electric field such as *E*app = 100 V/cm, EOF of PAA solution is similar to the Newtonian fluid, and the flow is in a steady state and symmetric about channel centerline. In addition, the induced polymeric stress *τxx* in the constriction microchannel is relatively small. When *E*app increases up to 400 V/cm, centerline streamlines start to show notable fluctuation and become asymmetric about the *x*-axis. At the entrance of the microchannel, a slight increase of *τxx* is observed, however, significant increase of *τxx* is observed near the microchannel walls and the downstream lips, as shown in Figure 14d. When *E*app increases to 500 V/cm, a pair of upstream vortices are induced at upstream of the microchannel, and the size of vortices is about the height of constriction microchannel. Additionally, as shown in Figure 14e, significant *τxx* is induced near the entrance of the constriction microchannel. However, when *E*app further increases to 600 V/cm, the size of the upstream vortices do

not show notable increase in comparison with that of *E*app = 500 V/cm. The results clearly show that increase of *c*<sup>p</sup> and/or *E*app can magnify the elastic instabilities of the viscoelastic EOF. However, the increase of *c*<sup>p</sup> has a more significant enhancing effect on the elastic instabilities of the viscoelastic EOF than the increase of *E*app. 500 V/cm. The results clearly show that increase of ୮ and/or ୟ୮୮ can magnify the elastic instabilities of the viscoelastic EOF. However, the increase of ୮ has a more significant enhancing effect on the elastic instabilities of the viscoelastic EOF than the increase of ୟ୮୮. the upstream vortices do not show notable increase in comparison with that of ୟ୮୮ = 500 V/cm. The results clearly show that increase of ୮ and/or ୟ୮୮ can magnify the elastic instabilities of the viscoelastic EOF. However, the increase of ୮ has a more significant enhancing effect on the elastic instabilities of the viscoelastic EOF than the increase of ୟ୮୮.

the microchannel, and the size of vortices is about the height of constriction microchannel. Additionally, as shown in Figure 14e, significant ௫௫ is induced near the entrance of the constriction microchannel. However, when ୟ୮୮ further increases to 600 V/cm, the size of the upstream vortices do not show notable increase in comparison with that of ୟ୮୮ =

the microchannel, and the size of vortices is about the height of constriction microchannel. Additionally, as shown in Figure 14e, significant ௫௫ is induced near the entrance of the constriction microchannel. However, when ୟ୮୮ further increases to 600 V/cm, the size of

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 16 of 29

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 16 of 29

**Figure 13.** Streamlines for 150 ppm PAA solution under different ୟ୮୮ at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm, (**d**) 400 V/cm, (**e**) 500 V/cm, (**f**) 600 V/cm. The color bar represents the elastic normal stress ௫௫. **Figure 13.** Streamlines for 150 ppm PAA solution under different *<sup>E</sup>*app at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm,(**d**) 400 V/cm, (**e**) 500 V/cm, (**f**) 600 V/cm. The color bar represents the elastic normal stress *<sup>τ</sup>xx*. **Figure 13.** Streamlines for 150 ppm PAA solution under different ୟ୮୮ at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm, (**d**) 400 V/cm, (**e**) 500 V/cm, (**f**) 600 V/cm. The color bar represents the elastic normal stress ௫௫.

**Figure 14.** Streamlines in microchannel of 150 ppm PAA solution under different ୟ୮୮ at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm, (**d**) 400 V/cm, (**e**) 500 V/cm, and (**f**) 600 V/cm. The color bar represents the elastic normal stress ௫௫. **Figure 14.** Streamlines in microchannel of 150 ppm PAA solution under different ୟ୮୮ at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm, (**d**) 400 V/cm, (**e**) 500 V/cm, and (**f**) 600 V/cm. The color bar represents the elastic normal stress ௫௫. **Figure 14.** Streamlines in microchannel of 150 ppm PAA solution under different *E*app at 1.78 s: (**a**) 100 V/cm, (**b**) 200 V/cm. (**c**) 300 V/cm, (**d**) 400 V/cm, (**e**) 500 V/cm, and (**f**) 600 V/cm. The color bar represents the elastic normal stress *τxx*.

Figure 15 depicts a flow map for the onset of vortices in unstable EOF as functions of ୮ and ୟ୮୮. At a fixed ୮ (ୟ୮୮), vortices and unstable EOF occur when ୟ୮୮ (୮) exceeds a certain threshold value. For example, for ୮ = 200 ppm, the flow becomes unstable with the occurrence of vortices when ୟ୮୮ is above 300 V/cm. At relatively low PAA concentration (i.e., ୮ = 100 ppm), it requires a very high electric field (up to 850 V/cm) to yield unstable EOF with upstream vortices. In contrast, at relatively high ୮ (i.e*.*, ୮ = Figure 15 depicts a flow map for the onset of vortices in unstable EOF as functions of ୮ and ୟ୮୮. At a fixed ୮ (ୟ୮୮), vortices and unstable EOF occur when ୟ୮୮ (୮) exceeds a certain threshold value. For example, for ୮ = 200 ppm, the flow becomes unstable with the occurrence of vortices when ୟ୮୮ is above 300 V/cm. At relatively low PAA concentration (i.e., ୮ = 100 ppm), it requires a very high electric field (up to 850 V/cm) to yield unstable EOF with upstream vortices. In contrast, at relatively high ୮ (i.e*.*, ୮ = Figure 15 depicts a flow map for the onset of vortices in unstable EOF as functions of *c*<sup>p</sup> and *E*app. At a fixed *c*<sup>p</sup> (*E*app), vortices and unstable EOF occur when *E*app (*c*p) exceeds a certain threshold value. For example, for *c*<sup>p</sup> = 200 ppm, the flow becomes unstable with the occurrence of vortices when *E*app is above 300 V/cm. At relatively low PAA concentration (i.e., *c*<sup>p</sup> = 100 ppm), it requires a very high electric field (up to 850 V/cm) to yield unstable EOF with upstream vortices. In contrast, at relatively high *c*<sup>p</sup> (i.e., *c*<sup>p</sup> = 500 ppm), the onset of vortices occurs at *E*app between 50 V/cm and 100 V/cm. An asymptotic curve fitting is implemented to illustrate the transition condition from no upstream vortices to the formation of upstream vortices, which is given as *<sup>E</sup>*app <sup>=</sup> 47.49 <sup>+</sup> 2892.25·0.987*c*<sup>p</sup> , where *c*<sup>p</sup> represents the polymer concentration in ppm, and *E*app is the apparent electric field in

V/cm. Above the curve in Figure 15, the EOF becomes time-dependent with upstream vortices, and no vortex forms under the conditions below the curve. Note that the flow map is only valid for the geometry considered in this study with the zeta potential of −110 mV. The instabilities of the viscoelastic EOF are dependent on the value of zeta potential. A comparison of the flow patterns of lower zeta potential (−70 mV) and higher zeta potential (−150 mV) for 150 ppm PAA solution under *E*app = 600 V/cm in the Appendix A clearly shows that higher zeta potential triggers stronger instabilities under the same condition. In addition, the dimensionless numbers (Reynolds number and Weissenberg number) of the studied conditions corresponding to Figure 15 are given in the Appendix A. with upstream vortices, and no vortex forms under the conditions below the curve. Note that the flow map is only valid for the geometry considered in this study with the zeta potential of −110 mV. The instabilities of the viscoelastic EOF are dependent on the value of zeta potential. A comparison of the flow patterns of lower zeta potential (−70 mV) and higher zeta potential (−150 mV) for 150 ppm PAA solution under ୟ୮୮ = 600 V/cm in the Appendix A clearly shows that higher zeta potential triggers stronger instabilities under the same condition. In addition, the dimensionless numbers (Reynolds number and Weissenberg number) of the studied conditions corresponding to Figure 15 are given in the Appendix A.

500 ppm), the onset of vortices occurs at ୟ୮୮ between 50 V/cm and 100 V/cm. An asymptotic curve fitting is implemented to illustrate the transition condition from no upstream vortices to the formation of upstream vortices, which is given as ୟ୮୮ = 47.49 + 2892.25 ∙ 0.987౦, where ୮ represents the polymer concentration in ppm, and ୟ୮୮ is the apparent electric field in V/cm. Above the curve in Figure 15, the EOF becomes time-dependent

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 17 of 29

**Figure 15.** Flow map in ୮-ୟ୮୮ space for EOF of PAA solutions flowing through a 10:1:10 constriction/expansion microchannel. Up-right of the fitting curve are the conditions that trigger the vortex in the EOF. **Figure 15.** Flow map in *c*p-*E*app space for EOF of PAA solutions flowing through a 10:1:10 constriction/expansion microchannel. Up-right of the fitting curve are the conditions that trigger the vortex in the EOF.

Since the flow velocity is time-dependent, we first calculate the cross-sectional average velocity over a period of ∆ = ଶ − ଵ, and then take the time-average to obtain the averaged velocity as, Since the flow velocity is time-dependent, we first calculate the cross-sectional average velocity over a period of ∆*t* = *t*<sup>2</sup> − *t*1, and then take the time-average to obtain the averaged velocity as,

$$\overline{\mathbf{U}} = \frac{\int\_{t\_1}^{t\_2} \int\_{-H\_\mathbb{C}/2}^{H\_\mathbb{C}/2} \mathbf{U}(\mathbf{0}, \mathbf{y}) dydt}{\Delta t \cdot H\_\mathbb{C}},\tag{11}$$

We choose ∆ = 1 s in the current study. Figure 16a shows the average velocity at the center of the constriction as a function of ୮ at different values of ୟ୮୮. In comparison, it also shows the result of Newtonian fluid whose viscosity is the same as the total viscosity of PAA solution under ୟ୮୮ = 600 V/cm. Under the same ୟ୮୮ = 600 V/cm, the average velocity in the constriction microchannel of the Newtonian flow is about 6%-12% higher than that of the PAA solution. The decrease of the average velocity in PAA solution is attributed to the induced polymeric stress at the entrance of the constriction microchannel. For Newtonian fluid, the average velocity decreases as ୮ increases, which is due to the increase of viscosity to make its viscosity be the same as that of PAA solution with the concentration of ୮. For PAA solutions, under the same ୟ୮୮ the average velocity exponentially decreases as the polymer concentration increases. One reason is attributed to the increase of total viscosity with the increase in ୮. In addition, the induced polymer stress within the constriction increases with the increase of polymer concentration, as shown in We choose ∆*t* = 1 s in the current study. Figure 16a shows the average velocity at the center of the constriction as a function of *c*<sup>p</sup> at different values of *E*app. In comparison, it also shows the result of Newtonian fluid whose viscosity is the same as the total viscosity of PAA solution under *E*app = 600 V/cm. Under the same *E*app = 600 V/cm, the average velocity in the constriction microchannel of the Newtonian flow is about 6–12% higher than that of the PAA solution. The decrease of the average velocity in PAA solution is attributed to the induced polymeric stress at the entrance of the constriction microchannel. For Newtonian fluid, the average velocity decreases as *c*<sup>p</sup> increases, which is due to the increase of viscosity to make its viscosity be the same as that of PAA solution with the concentration of *c*p. For PAA solutions, under the same *E*app the average velocity exponentially decreases as the polymer concentration increases. One reason is attributed to the increase of total viscosity with the increase in *c*p. In addition, the induced polymer stress within the constriction increases with the increase of polymer concentration, as shown in Figure 12, and the induced polymer stress slows down the flow.

Figure 12, and the induced polymer stress slows down the flow. Under the considered condition of *c*<sup>0</sup> = 0.01 mM, the EDL thickness is only 95 nm, which is much smaller than the height of the constriction. In Newtonian fluid, the EOF velocity can be approximated by the well-known Helmholtz–Smoluchowski velocity formula as described in Equation (10). We wonder if the average velocity of PAA solutions can be still approximated with the Helmholtz–Smoluchowski velocity formula. Since the actual local electric field within the constriction is much higher than the apparent

electric field *E*app, time-averaged electric field in the x-direction at the center of the constriction is used in the calculation of the Helmholtz–Smoluchowski velocity. Figure 16b shows the average velocity as a function of the apparent electric field under various PAA concentrations. The lines in Figure 16b represent the corresponding EOF velocity predicted by the Helmholtz–Smoluchowski velocity formula. At a fixed *c*p, as expected, the EOF velocity increases with an increase in the applied electric field. In general, the Helmholtz–Smoluchowski formula over predicts the velocity, and at a fixed *c*<sup>p</sup> the relative error increases with the increasing *E*app. For example, for *c*<sup>p</sup> = 100 ppm, the relative errors under *E*app =100 V/cm and 600 V/cm are, respectively, 1.2% and 8.6%. At a fixed *E*app, the absolute error, which is the difference between the Helmholtz–Smoluchowski approximated velocity and the average velocity obtained from the full numerical simulation, increases with the increasing PAA concentration. However, the relative error shows no notable change with the increasing *c*p. For example, at *E*app = 600 V/cm, the relative errors for *c*<sup>p</sup> = 100 ppm, 300 ppm, and 500 ppm are, respectively, 8.7%, 9.6%, and 9.0%. To evaluate the applicability of the Helmholtz–Smoluchowski formula to approximate the velocity of viscoelastic fluids for *c*<sup>p</sup> ranging from 100 ppm to 500 ppm, the minimum, average, and maximum relative errors at different *E*app are calculated as shown in Table 1. When *E*app ≤ 300 V/cm, the relative error is less than 5%. However, when *E*app ≥ 400 V/cm, the relative error is larger than 5%, and the Helmholtz–Smoluchowski formula failed to predict the velocity of the viscoelastic fluids accurately. In this study, the largest relative error is 9.6% when *c*<sup>p</sup> = 300 ppm and *E*app = 600 V/cm. field ୟ୮୮, time-averaged electric field in the x-direction at the center of the constriction is used in the calculation of the Helmholtz–Smoluchowski velocity. Figure 16b shows the average velocity as a function of the apparent electric field under various PAA concentrations. The lines in Figure 16b represent the corresponding EOF velocity predicted by the Helmholtz–Smoluchowski velocity formula. At a fixed ୮, as expected, the EOF velocity increases with an increase in the applied electric field. In general, the Helmholtz–Smoluchowski formula over predicts the velocity, and at a fixed ୮ the relative error increases with the increasing ୟ୮୮ . For example, for ୮ = 100 ppm , the relative errors under ୟ୮୮ =100 V/cm and 600 V/cm are, respectively, 1.2% and 8.6%. At a fixed ୟ୮୮, the absolute error, which is the difference between the Helmholtz–Smoluchowski approximated velocity and the average velocity obtained from the full numerical simulation, increases with the increasing PAA concentration. However, the relative error shows no notable change with the increasing ୮. For example, at ୟ୮୮ = 600 V/cm, the relative errors for ୮ = 100 ppm, 300 ppm, and 500 ppm are, respectively, 8.7%, 9.6%, and 9.0%. To evaluate the applicability of the Helmholtz–Smoluchowski formula to approximate the velocity of viscoelastic fluids for ୮ ranging from 100 ppm to 500 ppm, the minimum, average, and maximum relative errors at different ୟ୮୮ are calculated as shown in Table 1. When ୟ୮୮ 300 V/cm, the relative error is less than 5%. However, when ୟ୮୮ 400 V/cm, the relative error is larger than 5%, and the Helmholtz–Smoluchowski formula failed to predict the velocity of the viscoelastic fluids accurately. In this study, the largest relative error is 9.6% when ୮ = 300 ppm and ୟ୮୮ = 600 V/cm.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 18 of 29

Under the considered condition of = 0.01 mM, the EDL thickness is only 95 nm, which is much smaller than the height of the constriction. In Newtonian fluid, the EOF velocity can be approximated by the well-known Helmholtz–Smoluchowski velocity for-

can be still approximated with the Helmholtz–Smoluchowski velocity formula. Since the actual local electric field within the constriction is much higher than the apparent electric

**Figure 16.** Time averaged cross-sectional average velocity at the center of the constriction microchannel (*x* = 0): (**a**) Average velocity, (**b**) comparison of average velocity (with deviation) and **Figure 16.** Time averaged cross-sectional average velocity at the center of the constriction microchannel (*x* = 0): (**a**) Average velocity, (**b**) comparison of average velocity (with deviation) and Helmholtz–Smoluchowski velocity (lines).

**Table 1.** Relative error between the Helmholtz–Smoluchowski velocity and the average velocity **Table 1.** Relative error between the Helmholtz–Smoluchowski velocity and the average velocity from the full mathematical model.


#### **5. Conclusions**

Helmholtz–Smoluchowski velocity (lines).

Electroosmotic flow (EOF) of viscoelastic fluid through a 10:1:10 constriction microchannel is numerically investigated as functions of the applied electric field and the polymer concentration. In the current study, we neglect the effect of the polymer concentration on the zeta potential of the channel walls. Comparing to the EOF of Newtonian fluid, the following distinct results for viscoelastic EOF through a 10:1:10 constriction microchannel are obtained:


**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/mi12040417/s1. Video S1: The development of vortices for 500 ppm PAA solution under 100 V/cm apparent electric field. t<sup>0</sup> = 1.7 s and ∆t = 4.8 ms for each frame. Video S2: The development of vortices for 150 ppm PAA solution under 600 V/cm apparent electric field. t<sup>0</sup> = 1.7 s and ∆t = 2.0 ms for each frame.

**Author Contributions:** J.J. conducted simulations and drafted the manuscript; S.Q. supervised and revised the manuscript; and Z.L. analyzed and interpreted results. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded, in part, by the Foundation of State Key Laboratory of Coal Combustion (FSKLCCA1802).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Appendix A.1. Mesh Independence Study*

Three different meshes are used to conduct the mesh independent study with *E*app = 600 V/cm and *c*<sup>p</sup> = 100 ppm. As shown in Figure A1, there are 135,192, 95,252, and 77,192 cells in mesh 1, mesh 2, and mesh 3, respectively. For mesh 1, the meshes near the charged wall are 7 nm, so there are 14 meshes within the EDL thickness. In mesh 2 and mesh 3, the meshes near the charged wall are 10 nm, and there are 10 meshes within the EDL thickness. However, there are less meshes within the two reservoirs in mesh 3 than in mesh 2. Figure A2 shows spatial distribution of the normal polymer stress and streamlines of three different meshes at t = 1.78 s. The normal polymer stress and streamlines of three meshes show no notable difference.

of three meshes show no notable difference.

of three meshes show no notable difference.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 20 of 29

**Figure A1.** Three different meshes used for the mesh independence study. The meshes are symmetric with respect to the *x*-axis and *y*-axis, and only 1/4 of the total meshes are showed. (**a**) Mesh 1: 135.192 cells. (**b**) Mesh 2: 95.252 cells. (**c**) Mesh 3: 77.192 cells. **Figure A1.** Three different meshes used for the mesh independence study. The meshes are symmetric with respect to the *x*-axis and *y*-axis, and only 1/4 of the total meshes are showed. (**a**) Mesh 1: 135.192 cells. (**b**) Mesh 2: 95.252 cells. (**c**) Mesh 3: 77.192 cells. **Figure A1.** Three different meshes used for the mesh independence study. The meshes are symmetric with respect to the *x*-axis and *y*-axis, and only 1/4 of the total meshes are showed. (**a**) Mesh 1: 135.192 cells. (**b**) Mesh 2: 95.252 cells. (**c**) Mesh 3: 77.192 cells.

Three different meshes are used to conduct the mesh independent study with ୟ୮୮ = 600 V/cm and ୮ = 100 ppm. As shown in Figure A1, there are 135,192, 95,252, and 77,192 cells in mesh 1, mesh 2, and mesh 3, respectively. For mesh 1, the meshes near the charged wall are 7 nm, so there are 14 meshes within the EDL thickness. In mesh 2 and mesh 3, the meshes near the charged wall are 10 nm, and there are 10 meshes within the EDL thickness. However, there are less meshes within the two reservoirs in mesh 3 than in mesh 2. Figure A2 shows spatial distribution of the normal polymer stress and streamlines of three different meshes at t = 1.78 s. The normal polymer stress and streamlines

Three different meshes are used to conduct the mesh independent study with ୟ୮୮ = 600 V/cm and ୮ = 100 ppm. As shown in Figure A1, there are 135,192, 95,252, and 77,192 cells in mesh 1, mesh 2, and mesh 3, respectively. For mesh 1, the meshes near the charged wall are 7 nm, so there are 14 meshes within the EDL thickness. In mesh 2 and mesh 3, the meshes near the charged wall are 10 nm, and there are 10 meshes within the EDL thickness. However, there are less meshes within the two reservoirs in mesh 3 than in mesh 2. Figure A2 shows spatial distribution of the normal polymer stress and streamlines of three different meshes at t = 1.78 s. The normal polymer stress and streamlines

Figure A3 shows the spatial distribution of velocity magnitudes along *x* = 0 and *y* = 0 at t = 1.78 s. For velocity magnitude distribution along *x* = 0 (Figure A3a), the maximum relative error occurs at *y* = 0. The maximum relative error for mesh 2 is |Umesh2- Umesh1|/Umesh1 = 0.9%, and the maximum relative error for mesh 3 is |Umesh3-Umesh1|/Umesh1= 1.1%. The average relative errors for all sampling points are 0.45% for mesh 2 and 0.64% Figure A3 shows the spatial distribution of velocity magnitudes along *x* = 0 and *y* = 0 at t = 1.78 s. For velocity magnitude distribution along *x* = 0 (Figure A3a), the maximum relative error occurs at *y* = 0. The maximum relative error for mesh 2 is |Umesh2- Umesh1|/Umesh1 = 0.9%, and the maximum relative error for mesh 3 is |Umesh3-Umesh1|/Umesh1= 1.1%. The average relative errors for all sampling points are 0.45% for mesh 2 and 0.64% Figure A3 shows the spatial distribution of velocity magnitudes along *x* = 0 and *y* = 0 at t = 1.78 s. For velocity magnitude distribution along *x* = 0 (Figure A3a), the maximum relative error occurs at *y* = 0. The maximum relative error for mesh 2 is |Umesh2 − Umesh1|/Umesh1 = 0.9%, and the maximum relative error for mesh 3 is |Umesh3 − Umesh1|/Umesh1 = 1.1%. The average relative errors for all sampling points are 0.45% for mesh 2 and 0.64% for mesh 3. For velocity magnitude distribution along *y* = 0 (Figure A3b), comparing to mesh 1, the average relative errors for all sampling points are, respectively, 0.29% for mesh 2 and 0.72% for mesh 3.

Since the results from the above three meshes are in good agreement, we use mesh 3 to perform our other simulations.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 21 of 29

for mesh 3. For velocity magnitude distribution along *y* = 0 (Figure A3b), comparing to mesh 1, the average relative errors for all sampling points are, respectively, 0.29% for mesh

for mesh 3. For velocity magnitude distribution along *y* = 0 (Figure A3b), comparing to mesh 1, the average relative errors for all sampling points are, respectively, 0.29% for mesh

Since the results from the above three meshes are in good agreement, we use mesh 3

Since the results from the above three meshes are in good agreement, we use mesh 3

**Figure A3.** Spatial distribution of velocity magnitudes at t = 1.78 s. (**a**) Velocity magnitudes profile at *x* = 0. (**b**) Velocity magnitudes profile at *y* = 0. **Figure A3.** Spatial distribution of velocity magnitudes at t = 1.78 s. (**a**) Velocity magnitudes profile at *x* = 0. (**b**) Velocity magnitudes profile at *y* = 0. magnitudes profile at *y* = 0.

2 and 0.72% for mesh 3.

2 and 0.72% for mesh 3.

to perform our other simulations.

to perform our other simulations.

**Figure A4.** Boundary conditions with **n** denoting the normal unit vector on the surface. **Figure A4.** Boundary conditions with **n** denoting the normal unit vector on the surface.

#### **Figure A4.** Boundary conditions with **n** denoting the normal unit vector on the surface. *Appendix A.3. Analytical Solution of Afonso et al Appendix A.3. Analytical Solution of Afonso et al*

*Appendix A.3. Analytical Solution of Afonso et al* Based on PTT model, Afonso et al. [28] derived the analytical solution of mixed electro-osmotic/pressure driven flows of viscoelastic fluids between parallel plates based on the Debey–Hückel approximation. The analytical solution of the velocity profile across the Based on PTT model, Afonso et al. [28] derived the analytical solution of mixed electro-osmotic/pressure driven flows of viscoelastic fluids between parallel plates based on the Debey–Hückel approximation. The analytical solution of the velocity profile across the height of the channel is given as: Based on PTT model, Afonso et al. [28] derived the analytical solution of mixed electroosmotic/pressure driven flows of viscoelastic fluids between parallel plates based on the Debey–Hückel approximation. The analytical solution of the velocity profile across the height of the channel is given as:

*η*

$$u^{E} = \left(\frac{\epsilon\psi\_{0}E\_{x}}{\eta} - 2\overline{C}\kappa^{2}\epsilon\lambda^{2}\left[\frac{\epsilon\psi\_{0}E\_{x}}{\eta}\right]^{3}\right)(\overline{A}-1) + \frac{2}{3}\kappa^{2}\epsilon\lambda^{2}\left[\frac{\epsilon\psi\_{0}E\_{x}}{\eta}\right]^{3}\times(\overline{A}^{3}-1) \tag{A1}$$

where is the dielectric constant of the solution, ̅= ((cosh ())/(cosh ())), <sup>ଶ</sup> = ((2ଶଶ)/()), is the Boltzmann constant, and is the ionic density. ((2ଶଶ)/()), is the Boltzmann constant, and is the ionic density. where *e* is the dielectric constant of the solution, *A* = ((cosh(*κy*))/(cosh(*κH*))), *κ* <sup>2</sup> = 2*n*0*e* 2 *z* 2 /(*eκBT*) , *κ<sup>B</sup>* is the Boltzmann constant, and *n*<sup>0</sup> is the ionic density. For OB model, *ε* = 0. The Equation (A1) becomes *u <sup>E</sup>* = *eψ*0*E<sup>x</sup> A* − 1 .

*Appendix A.4. Summary of Current Studies*



*Micromachines* **2021**, *12*, 417


#### *Appendix A.5. Results for Zeta potentials of* −*70 mV and* −*150 mV Appendix A.5. Results for Zeta potentials of −70 mV and −150 mV*

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 1 of 29

Two other different zeta potentials (−70 mV and −150 mV) are studied for 150 ppm PAA solution under *E*app = 600 V/cm. Figure A5 shows the results of −70 mV zeta potential. Significant curvatures of the centerline streamlines are observed. However, at other places, no significant disturbance is observed. Figure A6 shows the results of −150 mV zeta potential. Similar to the results of −110 mV zeta potential, strong disturbance is induced in the viscoelastic EOF. Upstream vortices form and disappear. An increase of elastic normal stress is observed within the constriction microchannel. Figure A7 shows the velocity magnitude at the center of constriction channel (i.e., (0,0)). When the zeta potential is low (−70 mV), the velocity magnitude is almost steady state. However, when the zeta potential is high (−150 mV), the velocity magnitude shows strong fluctuation. The results for zeta potentials of −70 mV, −110 mV, and −150 mV show that the elastic instabilities of the EOF of PAA solutions are dependent on the value of zeta potential. Higher zeta potential induces larger electroosmotic velocity, and therefore stronger stretching of the polymers at the entrance of the constriction microchannel. Higher velocity and polymer normal stress lead to stronger instabilities of the viscoelastic flow. Two other different zeta potentials (−70 mV and -150 mV) are studied for 150 ppm PAA solution under ୟ୮୮ = 600 V/cm. Figure A5 shows the results of −70 mV zeta potential. Significant curvatures of the centerline streamlines are observed. However, at other places, no significant disturbance is observed. Figure A6 shows the results of −150 mV zeta potential. Similar to the results of -110 mV zeta potential, strong disturbance is induced in the viscoelastic EOF. Upstream vortices form and disappear. An increase of elastic normal stress is observed within the constriction microchannel. Figure A7 shows the velocity magnitude at the center of constriction channel (i.e., (0,0)). When the zeta potential is low (−70 mV), the velocity magnitude is almost steady state. However, when the zeta potential is high (−150 mV), the velocity magnitude shows strong fluctuation. The results for zeta potentials of −70 mV, −110 mV, and −150 mV show that the elastic instabilities of the EOF of PAA solutions are dependent on the value of zeta potential. Higher zeta potential induces larger electroosmotic velocity, and therefore stronger stretching of the polymers at the entrance of the constriction microchannel. Higher velocity and polymer normal stress lead to stronger instabilities of the viscoelastic flow.

**Figure A5.** Streamlines of 150 ppm PAA solution under ୟ୮୮ = 600 V/cm at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, (**f**) 1.80 s. Zeta potential is −70 mV. The color bar represents the elastic normal stress ௫௫. **Figure A5.** Streamlines of 150 ppm PAA solution under *E*app = 600 V/cm at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, (**f**) 1.80 s. Zeta potential is −70 mV. The color bar represents the elastic normal stress *τxx*.

## *Appendix A.6. Dimensionless Numbers*

Dimensional analysis is a useful tool to fully characterize the flow and identify the dominant forces in complex flows of polymeric materials. Reynolds number (*Re*) is commonly used in rheological studies to determine whether the inertial force or the viscous force is dominating the flow, which is given by: *Re* = *ρul*/*η*, where *ρ* is the fluid density, *u* is the average velocity in the microchannel, *l* is the characteristic length scale (*H*c), and *η* is the fluid viscosity. In addition, the Weissenberg number (*Wi*) is used to assess the flow elasticity of the PAA solutions, which is defined as: *Wi* = *uλ*/*l*, where λ is the relaxation time. Based on our results, the *Re* is nearly zero, indicating that the inertial force of the EOF is negligible. The *Wi* is on the order of 1, which indicates that the elastic effect in the flow is significant.

**References** 

**2015**, *440*, 126–132.

*1060*, 1–16.

*ter.* **2014**, *10*, 9016–9037.

*Matter.* **2018**, *14*, 3541–3549.

*Nanofluid.* **2006**, *2*, 481–492.

*Technol.* **2004**, *157*, 183–187.

**2011**, *166*, 1076–1079.

*Electrophoresis* **1994**, *15*, 623–626.

*Mach. Tool. Manu.* **2008**, *48*, 415–426.

*Colloid Interface Sci.* **2008**, *326*, 503–510.

crofluidics devices. *Electrophoresis* **2003**, *24*, 185–192.

control. *J. Micromech. Microeng.* **2011**, *21*, 085019.

nanopores using electro-osmotic flow. *J. Chem. Phys.* **2018**, *149*, 163311.

protein biomarkers with FraC nanopores. *Nat. Commun.* **2017**, *8*, 1–11.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 2 of 29

**Figure A6.** Streamlines of 150 ppm PAA solution under ୟ୮୮ = 600 V/cm at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, (**f**) 1.80 s. Zeta potential is −150 mV. The color bar represents the elastic normal stress ௫௫. **Figure A6.** Streamlines of 150 ppm PAA solution under *E*app = 600 V/cm at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, (**f**) 1.80 s. Zeta potential is −150 mV. The color bar represents the elastic normal stress *τxx*. **Figure A6.** Streamlines of 150 ppm PAA solution under ୟ୮୮ = 600 V/cm at different times: (**a**) 1.70 s, (**b**) 1.72 s, (**c**) 1.74 s, (**d**) 1.76 s, (**e**) 1.78 s, (**f**) 1.80 s. Zeta potential is −150 mV. The color bar represents the elastic normal stress ௫௫.

*Appendix A.6. Dimensionless Numbers*  Dimensional analysis is a useful tool to fully characterize the flow and identify the **Figure A7.** Velocity magnitudes at (0,0) for 150 ppm PAA solution with zeta potentials of −70 mV and −150 mV under ୟ୮୮ = 600 V/cm. **Figure A7.** Velocity magnitudes at (0, 0) for 150 ppm PAA solution with zeta potentials of −70 mV and −150 mV under *E*app = 600 V/cm. *Micromachines* **2021**, *12*, x FOR PEER REVIEW 3 of 29

**Figure A8.** Reynolds number (**a**) and Weissenberg number (**b**) of the viscoelastic EOF under the conditions of Figure 15. *Heat Fluid Flow* **2012**, *55*, 2044–2051. **Figure A8.** Reynolds number (**a**) and Weissenberg number (**b**) of the viscoelastic EOF under the conditions of Figure 15.39

3. Palyulin, V.V.; Ala-Nissila, T.; Metzler, R. Polymer translocation: The first two decades and the recent diversification. *Soft Mat-*

4. Takamura, Y.; Onoda, H.; Inokuchi, H.; Adachi, S.; Oki, A.; Horiike, Y. Low-voltage electroosmosis pump for stand-alone mi-

5. Li, L.; Wang, X.; Pu, Q.; Liu, S. Advancement of electroosmotic pump in microflow analysis: A review. *Anal. Chim. Acta* **2019**,

6. Jiang, H.; Weng, X.; Chon, C.H.; Wu, X.; Li, D. A microfluidic chip for blood plasma separation using electro-osmotic flow

7. Ermann, N.; Hanikel, N.; Wang, V.; Chen, K.; Weckman, N.E.; Keyser, U.F. Promoting single-file DNA translocations through

8. Huang, G.; Willems, K.; Soskine, M.; Wloka, C.; Maglia, G. Electro-osmotic capture and ionic discrimination of peptide and

9. Buyukdagli, S. Facilitated polymer capture by charge inverted electroosmotic flow in voltage-driven polymer translocation. *Soft* 

10. Bello, M.S.; Besi, P.D.; Rezzonico, R.; Righetti, P.G.; Casiraghi, E. Electroosmosis of polymer solutions in fused silica capillaries.

13. Zimmerman, W.; Rees, J.; Craven, T. Rheometry of non-Newtonian electrokinetic flow in a microchannel T-junction. *Microfluid.* 

15. Das, M.; Jain, V.; Ghoshdastidar, P. Fluid flow analysis of magnetorheological abrasive flow finishing (MRAFF) process. *Int. J.* 

*17.* Thien, N.P.; Tanner, R.I. A new constitutive equation derived from network theory. *J. Non-Newton. Fluid Mech.* **1977**, *2*, 353–365. 18. Koszkul, J.; Nabialek, J. Viscosity models in simulation of the filling stage of the injection molding process. *J. Mater. Process.* 

19. Öztekin, A.; Brown, R.A.; McKinley, G.H. Quantitative prediction of the viscoelastic instability in cone-and-plate flow of a Boger

20. Das, S.; Chakraborty, S. Analytical solutions for velocity, temperature and concentration distribution in electroosmotic micro-

*21.* Zhao, C.; Zholkovskij, E.; Masliyah, J.H.; Yang, C. Analysis of electroosmotic flow of power-law fluids in a slit microchannel. *J.* 

22. Zhao, C.; Yang, C. An exact solution for electroosmosis of non-Newtonian fluids in microchannels. *J. Non-Newton. Fluid Mech.* 

23. Zhao, C.; Yang, C. Joule heating induced heat transfer for electroosmotic flow of power-law fluids in a microcapillary. *Int. J.* 

11. Chang. F.M.; Tsao, H.K. Drag reduction in electro-osmosis of polymer solutions. *Appl. Phys. Lett.* **2007**, *90*, 194105. 12. Kamişli, F. Flow analysis of a power-law fluid confined in an extrusion die. *Int. J. Eng. Sci.* **2003**, *41*, 1059–1083.

14. Hakim, A. Mathematical analysis of viscoelastic fluids of White-Metzner type. *J. Math. Anal. Appl.* **1994**, *185*, 675–705.

16. Oldroyd, J.G. On the formulation of rheological equations of state. *Proc. Math. Phys. Eng. Sci.* **1950**, *200*, 523–541.

fluid using a multi-mode Giesekus model. *J. Non-Newton. Fluid Mech.* **1994**, *54*, 351–377.

channel flows of a non-Newtonian bio-fluid. *Anal. Chim. Acta* **2006**, *559*, 15–24.

## **References**

