2.1 System description

We consider a simple planner surface, as shown in figure 1. The origin is assumed to lie on this surface, while the $x$ and the $z$ axes run parallel to the surface. The vertical axis is taken as the $y$ axis. This surface bears a surface charge ( $\unicode[STIX]{x1D70E}^{\prime }(x^{\prime },z^{\prime })$ ), which is non-uniform in nature, being an arbitrary function of $x$ and $z$ . The axis system can be chosen such that the surface charge has the simplest algebraic form. On top of the surface, a second-order fluid rests and has the following properties: viscosity $\unicode[STIX]{x1D702}$ , permittivity $\unicode[STIX]{x1D700}$ and retarded motion constants $\unicode[STIX]{x1D6FC}_{1}$ and $b\unicode[STIX]{x1D6FC}_{1}$ . This fluid also contains electrolytes which dissociate into ions and hence form an EDL over the solid surface. The reference or, bulk density of the ions are taken to be $c_{0}^{\prime }$ . An external electric field of magnitude $E_{0}$ is applied parallel to the surface, which actuates EOF. We assume that this potential makes an angle $\unicode[STIX]{x1D703}$ with the $x$ axis. We only consider a steady state EOF in the present work. Note that, the direction of the applied external field is denoted by the unit vector $\hat{\boldsymbol{e}}$ , where $\hat{\boldsymbol{e}}=\hat{\boldsymbol{e}}_{x}\cos (\unicode[STIX]{x1D703})+\hat{\boldsymbol{e}}_{z}\sin (\unicode[STIX]{x1D703})$ . Here $\hat{\boldsymbol{e}}_{x}$ and $\hat{\boldsymbol{e}}_{z}$ are unit vectors along the $x$ and $z$ axes respectively.

Figure 1. A schematic representation of the flow geometry. The origin and the axis system rest on the surface. This surface bears a non-uniform charge of the form, $\unicode[STIX]{x1D70E}^{\prime }(x^{\prime },z^{\prime })$ . An electric field of magnitude $E_{0}$ (represented by the thick red line) is applied parallel to the surface in a direction $\hat{\boldsymbol{e}}$ . This unit vector $\hat{\boldsymbol{e}}$ makes an angle $\unicode[STIX]{x1D703}$ with the $x$ axis. The fluid here obeys a second-order constitute relationship.

Since we are interested in finding the resulting electroosmotic slip velocity, we assume the EDL thickness to be much smaller compared to the typical length scale of variation of the surface charge. This assumption is quite realistic and has been applied to a number of previously executed studies (Ajdari Reference Ajdari1996; Khair & Squires Reference Khair and Squires2008; Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012c ). The typical EDL thickness is given by, $\unicode[STIX]{x1D706}_{D}=\sqrt{(\unicode[STIX]{x1D700}kT)/(2\bar{z}^{2}e^{2}c_{0}^{\prime })}$ , where $k$ is the Boltzmann constant, $T$ is the absolute temperature, $\bar{z}$ is the valance of the ions and $e$ is the proton charge. If we assume that the length scale of variation in the surface charge is $d$ , then the thin EDL limit would indicate $\unicode[STIX]{x1D706}_{D}\ll d$ .

Before we move towards the detailed formulation, it is important to go through the major assumptions, apart from the thin EDL approximation as mentioned above, under which the present analysis remains valid. First, the flow is steady and the inertial effects are assumed to be negligible. Second, we assume that there are no bulk-scale imposed gradients in the salt concentration. It has been shown by a number of previous studies (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2001; Khair & Squires Reference Khair and Squires2008; Yariv Reference Yariv2009) that such gradients lead to additional flows inside the thin double layer which eventually alter the slip velocity. This additional slip velocity is sometimes referred to as chemioosmotic slip (Khair & Squires Reference Khair and Squires2008) and it drives a salt flux to make the salt concentration uniform everywhere. Third, we assume that the Dukhin number ( $Du$ ), which characterizes the strength of conduction within the thin EDL (Dukhin Reference Dukhin1993; Squires & Bazant Reference Squires and Bazant2004; Khair & Squires Reference Khair and Squires2008) (also referred to as the ‘surface current’), is small enough so that the conduction within the EDL may be neglected. Previously, Khair & Squires (Reference Khair and Squires2008) have shown that, if the applied external field is strong enough, it will drive an additional conduction current within the EDL which will lead to bulk concentration gradients and hence chemioosmotic slip. Such effects are neglected in the present study. In fact, Khair & Squires (Reference Khair and Squires2008), along with Yariv (Reference Yariv2009), demonstrated that in the absence of any imposed bulk-scale gradients, concentration polarization does not occur for thin EDLs provided the applied field, advection strength and the zeta potential are not too large. The final results are valid for low surface potential (usually represented by, $\unicode[STIX]{x1D701}\leqslant 25~\text{mV}$ ), in the limit, $\unicode[STIX]{x1D6FF}\ll De<O(1)$ ( $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D706}_{D}/d$ ). Finally, we note that the whole formulation breaks down, if the zeta potential is large compared to the thermal voltage ( $kT/\bar{z}e$ ) (Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012a ). The present model also fails if the applied electric field becomes asymptotically large (Schnitzer & Yariv Reference Schnitzer and Yariv2012c ). Yariv and coworkers (Schnitzer & Yariv Reference Schnitzer and Yariv2012a ,Reference Schnitzer and Yariv c ) have previously shown that such extreme cases are often associated with nested boundary layer structures, leading to multiple distinguished limits (Bender & Orszag Reference Bender and Orszag1999). Later, we will assume that the Deborah number ( $De$ ) remains small enough so that a regular perturbation expansion remains possible, taking $De$ as the gauge function (Bird et al. Reference Bird, Armstrong and Hassager1987).

Note that, although in the present study we take a semi-infinite fluid domain, the same analysis remains valid for a channel as well. This can be attributed to the thin EDL approximation, by virtue of which, the far-field conditions would effectively represent conditions at the channel centreline. We end this subsection with a brief discussion on the possible range of values of Deborah number ( $De$ ). For steady flows, this number is usually defined as (Bird et al. Reference Bird, Armstrong and Hassager1987), $De=t_{m}K$ , where $K$ is the characteristic strain rate and $t_{m}$ is the relaxation time (largest). In the present context, the definition becomes, $De=t_{m}u_{ref}/d$ , $u_{ref}$ being the characteristic velocity (here $u_{ref}\sim \unicode[STIX]{x1D700}\unicode[STIX]{x1D701}_{0}^{2}/\unicode[STIX]{x1D702}d$ ), while $t_{m}=\unicode[STIX]{x1D6FC}_{1}/\unicode[STIX]{x1D702}$ ( $\unicode[STIX]{x1D6FC}_{1}$ is a material constant and $\unicode[STIX]{x1D702}$ is the viscosity). For typical viscoelastic fluids, $t_{m}$ has values ranging from milliseconds (ms) to seconds (s) (Ardekani et al. Reference Ardekani, Joseph, Dunn-Rankin and Rangel2009; Brust et al. Reference Brust2013). As an example, Brust et al. (Reference Brust2013) have previously shown that the relaxation time for blood is of the order of 7 ms, at $37\,^{\circ }\text{C}$ . Hence, for $u_{ref}\sim 10^{-4}~\text{m}~\text{s}^{-1}$ , $d\sim 10^{-4}{-}10^{-6}~\text{m}$ , one deduces, $De\sim O(1)-O(10^{-3})$ . However, strictly speaking, a second-order constitutive relation is not applicable for large values of the Deborah number (Bird et al. Reference Bird, Armstrong and Hassager1987).

2.2 Governing equations

As already mentioned, we will closely follow the framework previously derived by Yariv and coworkers (Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012c ). To avoid repetition of the well-known equations in electrokinetics, we will directly start with the non-dimensional ones. The non-dimensionalization has been carried out following the work of Saville (Reference Saville1977) and Yariv and coworkers (Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012a ,Reference Schnitzer and Yariv c , Reference Schnitzer and Yariv2014). To this end, we first select the characteristic scales for the relevant variables. The characteristic scale of any variable $\unicode[STIX]{x1D709}$ will be denoted by $\unicode[STIX]{x1D709}_{C}$ . As the length scale, we chose, $x_{C}=y_{C}=z_{C}=d$ , which is the typical scale over which the surface charge varies. For example, if the charge is periodic in nature, $d$ would be the wavelength of modulation. The potential scales as, $\unicode[STIX]{x1D713}_{C}\sim \unicode[STIX]{x1D701}_{0}$ , where $\unicode[STIX]{x1D701}_{0}$ denotes the typical magnitude of surface potential. We can relate the potential $\unicode[STIX]{x1D701}_{0}$ with the surface charge density. If $\unicode[STIX]{x1D70E}_{c}$ is the typical order of magnitude of the surface charge, then one can define $\unicode[STIX]{x1D701}_{0}$ as, $\unicode[STIX]{x1D701}_{0}\sim (2kT/\bar{z}e)\sinh ^{-1}[(\bar{z}e\unicode[STIX]{x1D70E}_{c}\unicode[STIX]{x1D706}_{D})/(2\unicode[STIX]{x1D700}kT)]$ . For small values of surface charge, the foregoing relation can be linearized to $\unicode[STIX]{x1D701}_{0}\sim \unicode[STIX]{x1D70E}_{c}\unicode[STIX]{x1D706}_{D}/\unicode[STIX]{x1D700}$ , as previously noted by Ajdari (Reference Ajdari1996). Note that the choice of potential here is slightly different to that of Saville’s (Reference Saville1977) work, where the thermal scale was chosen as the reference value for the potential. The connection (or, equivalence) between the present choice for potential and that of the thermal scale has been discussed in more detail in appendix A.

The concentration scales as, $c_{1,C}=c_{2,C}\sim c_{0}^{\prime }$ . The velocity scales as (Saville Reference Saville1977; Yariv Reference Yariv2009), $u_{C}=v_{C}=w_{C}\sim (\unicode[STIX]{x1D700}\unicode[STIX]{x1D701}_{0}^{2}/\unicode[STIX]{x1D702}d)$ , while the stresses scale as, $\unicode[STIX]{x1D70F}_{C}\sim \unicode[STIX]{x1D702}u_{C}/d$ . Finally, the pressure scales as, $p_{C}=\unicode[STIX]{x1D702}u_{C}/d$ . We non-dimensionalize the variables with their characteristic scales, i.e. for any variable $\unicode[STIX]{x1D709}^{\prime }$ it’s non-dimensional version reads, $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}^{\prime }/\unicode[STIX]{x1D709}_{C}$ . We further define total concentration and charge density as follows: (Mortensen et al. Reference Mortensen2005): $\unicode[STIX]{x1D70C}=c_{1}-c_{2}$ and $c=c_{1}+c_{2}$ . Enforcing this non-dimensionalization scheme yields the following governing equations (Saville Reference Saville1977):

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle Pe(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c)=\unicode[STIX]{x1D6FB}^{2}c+\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})-\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D6FD}(\boldsymbol{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}) & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle Pe(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C})=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}+\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(c\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})-\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D6FD}(\boldsymbol{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c) & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}=-\frac{\unicode[STIX]{x1D70C}}{2\bar{\unicode[STIX]{x1D701}}_{0}} & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}-\unicode[STIX]{x1D6FD}\hat{\boldsymbol{e}})=0;\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0. & \displaystyle\end{eqnarray}$$

First note that, in (2.1)–(2.4), the potential $\unicode[STIX]{x1D711}$ is the contribution from only the diffuse charges. The total non-dimensional potential can be written as (Bahga et al. Reference Bahga, Vinogradova and Bazant2010), $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D711}-\unicode[STIX]{x1D6FD}(\hat{\boldsymbol{e}}\boldsymbol{\cdot }\boldsymbol{r})$ . In the foregoing expression, the component $\unicode[STIX]{x1D6FD}(\hat{\boldsymbol{e}}\boldsymbol{\cdot }\boldsymbol{r})$ represents the contribution from the externally applied field. Note that the non-dimensional parameter $\unicode[STIX]{x1D6FD}$ is defined as, $\unicode[STIX]{x1D6FD}=E_{0}d/\unicode[STIX]{x1D701}_{0}$ , which denotes the strength of the applied electric field with respect to the EDL potential. Additionally, in (2.1)–(2.2) $Pe$ is the ionic Péclet number, defined as $Pe=u_{C}d/D=\unicode[STIX]{x1D700}\unicode[STIX]{x1D701}_{0}^{2}/\unicode[STIX]{x1D702}D$ ; while $\bar{\unicode[STIX]{x1D701}}_{0}$ is the strength of the zeta potential defined as, $\bar{\unicode[STIX]{x1D701}}_{0}=\bar{z}e\unicode[STIX]{x1D701}_{0}/kT$ .

In (2.4), the stress tensor for a second-order fluid is defined as (in dimensional form) (Leal Reference Leal1975; Joseph & Feng Reference Joseph and Feng1996), $\unicode[STIX]{x1D749}^{\prime }=2\unicode[STIX]{x1D702}\unicode[STIX]{x1D63C}^{\prime }+\unicode[STIX]{x1D64E}^{\prime }$ . Here $\unicode[STIX]{x1D63C}^{\prime }$ is the rate of strain tensor ( $=[\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime })^{\text{T}}]/2$ ) and $\unicode[STIX]{x1D64E}^{\prime }$ denotes the excess polymeric stresses given by, $\unicode[STIX]{x1D64E}^{\prime }=\unicode[STIX]{x1D6FC}_{1}(\unicode[STIX]{x1D63D}+b\unicode[STIX]{x1D63C}^{2})$ , where $\unicode[STIX]{x1D6FC}_{1}$ and $b$ are material constants (Bird et al. Reference Bird, Armstrong and Hassager1987). The tensor $\unicode[STIX]{x1D63D}^{\prime }$ is related to the strain rate as follows (analogous to a convected derivative) (Joseph & Feng Reference Joseph and Feng1996): $\unicode[STIX]{x1D63D}^{\prime }=\unicode[STIX]{x2202}\unicode[STIX]{x1D63C}^{\prime }/\unicode[STIX]{x2202}t^{\prime }+\boldsymbol{v}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63C}^{\prime }+\unicode[STIX]{x1D63C}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime })^{\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D63C}^{\prime }$ . In non-dimensional form (following the aforementioned scheme), the stresses become: $\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D63C}+De(\unicode[STIX]{x1D64E})$ . Here $De$ is the Deborah number defined as (Bird et al. Reference Bird, Armstrong and Hassager1987), $De=\unicode[STIX]{x1D6FC}_{1}u_{C}/\unicode[STIX]{x1D702}d$ . With the above considerations, the momentum equation (2.4) can be rewritten as follows:

(2.5) $$\begin{eqnarray}-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+De(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E})+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}-\unicode[STIX]{x1D6FD}\hat{\boldsymbol{e}})=0;\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0.\end{eqnarray}$$

Equations (2.1)–(2.5) are subject to the following boundary conditions:

(2.6a-d ) $$\begin{eqnarray}\text{At}\quad y=0,\quad \boldsymbol{v}=0,\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}y}=\unicode[STIX]{x1D6FF}^{-1}\unicode[STIX]{x1D70E}(x,z)\quad \text{and}\quad \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}+\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}y}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}y}+\bar{\unicode[STIX]{x1D701}}_{0}c\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}y}=0;\end{eqnarray}$$

(2.7a-g ) $$\begin{eqnarray}\text{At}\quad y\rightarrow \infty ,\quad \unicode[STIX]{x1D711}\rightarrow 0,\quad c\rightarrow 2,\quad \unicode[STIX]{x1D70C}\rightarrow 0,\quad \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\rightarrow 0,\quad v\rightarrow 0\quad \text{and}\quad \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}y}\rightarrow 0.\end{eqnarray}$$

In (2.6), we define (Schnitzer & Yariv Reference Schnitzer and Yariv2012c ) $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}^{\prime }\unicode[STIX]{x1D706}_{D}/\unicode[STIX]{x1D700}\unicode[STIX]{x1D701}_{0}$ . The above set of equations (2.1)–(2.7) completes the description for the EOF of second-order fluids. These are a formidable set of equations and to the best of our knowledge analytical solutions to these are not possible. However, one can look for asymptotic solutions in the thin EDL regime ( $\unicode[STIX]{x1D6FF}\ll 1$ ) and small $De$ values ( $De<O(1)$ ). We elaborate on the asymptotic solutions in the next section.

3.1 Equations in the outer layer

The leading-order outer layer equations read:

(3.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70C}=0;\quad Pe(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c)=\unicode[STIX]{x1D6FB}^{2}c;\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }(c\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})-\boldsymbol{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=0\end{eqnarray}$$

(3.3) $$\begin{eqnarray}-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+De(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E})+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}-\unicode[STIX]{x1D6FD}\boldsymbol{e})=0;\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0.\end{eqnarray}$$

Similar equations have already been derived previously by Yariv (Reference Yariv2009), albeit for Newtonian fluids. These are subject to the far-field conditions, as mentioned in (2.7). We must also require the outer variables to match with the inner ones at the edge of the EDL. The matching conditions are discussed in the next subsection. However, based on overall charge and salt concentration balance conditions within the EDL, we can derive bulk-scale conditions for the outer layer concentration and potential to the leading order in $\unicode[STIX]{x1D6FF}$ . It can be shown (Yariv Reference Yariv2009) that the said bulk-scale conditions eventually translate to the following form, provided $\unicode[STIX]{x1D6FD}$ and $Pe$ are not asymptotically large:

(3.4a,b ) $$\begin{eqnarray}\lim _{y\rightarrow 0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}y}=\lim _{y\rightarrow 0}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}=0\quad \text{and}\quad \lim _{y\rightarrow 0}v=0.\end{eqnarray}$$

Taking a hint from (3.4) and in the absence of any imposed bulk-scale concentration gradients, equations (3.2) have the simple trivial solutions:

(3.5a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D711}=\unicode[STIX]{x1D70C}=0;\quad c=2.\end{eqnarray}$$

Similar solutions have previously been derived by Yariv (Reference Yariv2009) and Bazant and coworkers (Bazant, Thornton & Ajdari Reference Bazant, Thornton and Ajdari2004; Højgaard Olesen et al. Reference Højgaard Olesen, Bazant and Bruus2010). The above solution essentially indicates that to leading order in $\unicode[STIX]{x1D6FF}$ the bulk is charge free and has uniform concentration everywhere. In light of (3.5), the outer layer momentum equations can be written in the following form:

(3.6) $$\begin{eqnarray}-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+De(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E})=0;\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0.\end{eqnarray}$$

The velocities and the pressure are subject to the matching conditions at the edge of the EDL, discussed below.

3.2 Equations in the inner layer

To derive the inner layer equations, we first need to rescale the variables. To this end, we define the following inner layer variables (rescaled for the thin EDL limit), $Y=y/\unicode[STIX]{x1D6FF}$ ; $X=x$ ; $Z=z$ ; $U=u$ ; $W=w$ ; $V=v/\unicode[STIX]{x1D6FF}$ ; $\tilde{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D711}$ ; $\tilde{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}$ ; $\tilde{c}=c$ . The pressure is rescaled as (Yariv Reference Yariv2009), $P=p/\unicode[STIX]{x1D6FF}^{2}$ . The excess stresses ( $\unicode[STIX]{x1D64E}$ ) are rescaled as follows: $\unicode[STIX]{x1D61A}_{xx}=\tilde{\unicode[STIX]{x1D61A}}_{xx}/\unicode[STIX]{x1D6FF}^{2}$ ; $\unicode[STIX]{x1D61A}_{xy}=\tilde{\unicode[STIX]{x1D61A}}_{xy}/\unicode[STIX]{x1D6FF}$ ; $\unicode[STIX]{x1D61A}_{xz}=\tilde{\unicode[STIX]{x1D61A}}_{xz}/\unicode[STIX]{x1D6FF}^{2}$ ; $\unicode[STIX]{x1D61A}_{yz}=\tilde{\unicode[STIX]{x1D61A}}_{yz}/\unicode[STIX]{x1D6FF}$ ; $\unicode[STIX]{x1D61A}_{yy}=\tilde{\unicode[STIX]{x1D61A}}_{yy}/\unicode[STIX]{x1D6FF}^{2}$ ; and $\unicode[STIX]{x1D61A}_{zz}=\tilde{\unicode[STIX]{x1D61A}}_{zz}/\unicode[STIX]{x1D6FF}^{2}$ . Note that we denote the variables in the inner layer either by uppercase letters or with a ‘tilde’ sign overhead. Similar scaling for the stresses have also been used by Saprykin, Koopmans & Kalliadasis (Reference Saprykin, Koopmans and Kalliadasis2007) for studying free surface flows of viscoelastic fluids. Previously, Yariv and coworkers (Yariv Reference Yariv2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012c ), Squires and coworkers (Squires & Bazant Reference Squires and Bazant2004; Khair & Squires Reference Khair and Squires2008) and many others (Ghosh & Chakraborty Reference Ghosh and Chakraborty2015) have demonstrated that, to leading order in $\unicode[STIX]{x1D6FF}$ , the charge, concentration and the potential within the EDL are independent of the velocities (provided the assumptions stated in the previous sections are met). Hence, we can easily infer that $\tilde{\unicode[STIX]{x1D70C}}$ and $\tilde{c}$ will satisfy the Boltzmann distribution, while the potential would satisfy the Poisson–Boltzmann equation (Hunter Reference Hunter1981). For a detailed derivation of the said equations, one can refer to the following works Mortensen et al. (Reference Mortensen2005), Yariv (Reference Yariv2009), Schnitzer & Yariv (Reference Schnitzer and Yariv2012c ). Previously Schnitzer & Yariv (Reference Schnitzer and Yariv2012c ) have shown that to leading order in $\unicode[STIX]{x1D6FF}$ , the ‘specified surface charge’ condition for the potential can be conveniently replaced by a ‘specified potential’ condition, which reads: $\tilde{\unicode[STIX]{x1D711}}(X,Y=0,Z)=\unicode[STIX]{x1D701}(X,Z)$ , where the zeta potential $\unicode[STIX]{x1D701}(X,Z)$ is an arbitrary function of $x$ and $z$ . The solution for the potential reads (Schnitzer & Yariv Reference Schnitzer and Yariv2012c ):

(3.7) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}=\frac{2}{\bar{\unicode[STIX]{x1D701}}_{0}}\ln \left[\frac{1+\text{e}^{-Y}\tanh (\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D701}/4)}{1-\text{e}^{-Y}\tanh (\bar{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D701}/4)}\right].\end{eqnarray}$$

The above solution remains valid beyond low values of $\bar{\unicode[STIX]{x1D701}}_{0}$ as well. For low surface potential, i.e. within the paradigm of Debye–Huckel linearization, the above equation can be simplified to:

(3.8) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D701}(X,Z)\text{e}^{-Y}+O(\bar{\unicode[STIX]{x1D701}}_{0}^{2}).\end{eqnarray}$$

We can now directly move on to the momentum equations inside the EDL. In presence of excess polymeric stresses, they take the following form:

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}Y^{2}}+De\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xx}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xy}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xz}}{\unicode[STIX]{x2202}Z}\right)+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}X}-\unicode[STIX]{x1D6FD}\cos \unicode[STIX]{x1D703}\right)=0 & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Y}+De\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{yy}}{\unicode[STIX]{x2202}Y}\right)+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y}=0 & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}^{2}W}{\unicode[STIX]{x2202}Y^{2}}+De\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xz}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{yz}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{zz}}{\unicode[STIX]{x2202}Z}\right)+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Z}-\unicode[STIX]{x1D6FD}\sin \unicode[STIX]{x1D703}\right)=0 & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Z}=0. & \displaystyle\end{eqnarray}$$

These equations are subject to the following boundary conditions:

(3.13a,b ) $$\begin{eqnarray}\text{At}\quad Y=0,\quad U=V=W=0.\end{eqnarray}$$

At the edge of the EDL the following matching conditions have to be satisfied (to leading order in $\unicode[STIX]{x1D6FF}$ ) for the velocities and the pressure:

(3.14a-d ) $$\begin{eqnarray}\lim _{Y\rightarrow \infty }U=\lim _{y\rightarrow 0}u;\quad \lim _{Y\rightarrow \infty }W=\lim _{y\rightarrow 0}w;\quad \lim _{y\rightarrow 0}v=0;\quad \lim _{Y\rightarrow \infty }P=0.\end{eqnarray}$$

The different components of the excess stress tensor ( $\tilde{\unicode[STIX]{x1D61A}}_{xx}$ , $\tilde{\unicode[STIX]{x1D61A}}_{xy}\ldots$ etc.) have the following expressions, to leading order in  $\unicode[STIX]{x1D6FF}$ :

(3.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D61A}}_{xx}=\left(1+\frac{b}{4}\right)\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\right)^{2};\quad \tilde{\unicode[STIX]{x1D61A}}_{xz}=\left(1+\frac{b}{4}\right)\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y};\\ \displaystyle \tilde{\unicode[STIX]{x1D61A}}_{yy}=\frac{b}{4}\left[\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)^{2}\right];\quad \tilde{\unicode[STIX]{x1D61A}}_{zz}=\left(1+\frac{b}{4}\right)\left(\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)^{2}\\ \begin{array}{@{}rcl@{}}\displaystyle \tilde{\unicode[STIX]{x1D61A}}_{xy}\ & =\ & \displaystyle \frac{1}{2}\left(U\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}X\unicode[STIX]{x2202}Y}+V\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}Y^{2}}+W\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}Z\unicode[STIX]{x2202}Y}\right)+\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y}+\frac{1}{2}\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Z}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}-\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Z}\right)\\ \ & \ & \displaystyle +\,\frac{b}{4}\left\{\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}X}\right)\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}-2\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Z}\right\}\end{array}\\ \begin{array}{@{}rcl@{}}\displaystyle \tilde{\unicode[STIX]{x1D61A}}_{yz}\ & =\ & \displaystyle \frac{1}{2}\left(U\frac{\unicode[STIX]{x2202}^{2}W}{\unicode[STIX]{x2202}X\unicode[STIX]{x2202}Y}+V\frac{\unicode[STIX]{x2202}^{2}W}{\unicode[STIX]{x2202}Y^{2}}+W\frac{\unicode[STIX]{x2202}^{2}W}{\unicode[STIX]{x2202}Z\unicode[STIX]{x2202}Y}\right)+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y}+\frac{1}{2}\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}X}-\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)\\ \ & \ & \displaystyle +\,\frac{b}{4}\left\{\left(\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}X}\right)\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}-2\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right\}.\end{array}\end{array}\right\}\end{eqnarray}$$

Here, we stick to the limit of $\unicode[STIX]{x1D6FD}\sim O(1)$ , as attributable to the choice of the velocity scale. The most notable feature of the scaling in the inner layer is that, for Newtonian fluids, the shear stress components vary as (non-dimensional) $\unicode[STIX]{x1D70F}_{xy}\sim \unicode[STIX]{x1D70F}_{yz}\sim O(\unicode[STIX]{x1D6FF}^{-1})$ , whereas the normal stress components scale as $\unicode[STIX]{x1D70F}_{xx}\sim \unicode[STIX]{x1D70F}_{yy}\sim \unicode[STIX]{x1D70F}_{zz}\sim \unicode[STIX]{x1D70F}_{xz}\sim O(1)$ inside the EDL. In sharp contrast, for a second-order fluid, the normal stress components ( $\unicode[STIX]{x1D70F}_{xx}$ , $\unicode[STIX]{x1D70F}_{zz}$ , etc.) and the shear stress component $\unicode[STIX]{x1D70F}_{xz}$ show the following variations in the inner layer, $\unicode[STIX]{x1D70F}_{xx}\sim \unicode[STIX]{x1D70F}_{yy}\sim \unicode[STIX]{x1D70F}_{zz}\sim \unicode[STIX]{x1D70F}_{xz}\sim O(\unicode[STIX]{x1D6FF}^{-2})$ . All other stresses, i.e. the remaining shear stresses scale similarly to those of Newtonian fluids, $\unicode[STIX]{x1D70F}_{xy}\sim \unicode[STIX]{x1D70F}_{yz}\sim O(\unicode[STIX]{x1D6FF}^{-1})$ . Therefore, the normal stresses, as well as the shear component $\unicode[STIX]{x1D70F}_{xz}$ , are asymptotically large compared to their counterparts in the bulk. This scaling pattern has been verified through comparison with numerical simulations (§ 5.1). Therefore, these stress components explicitly appear in the momentum equations within the EDL, even to leading order in $\unicode[STIX]{x1D6FF}$ . We later show that the aforementioned stress components play an important role in altering the overall flow dynamics, eventually leading to anisotropic effects. This is one of the major differences between the dynamics of Newtonian and second-order fluids, within the EDL.

Finally, we note that the left-hand side of the first two conditions in (3.14) represents the slip velocity at the edge of the EDL, along the plane surface. In other words, $\lim _{Y\rightarrow \infty }U$ is the slip velocity along the $x$ axis while $\lim _{Y\rightarrow \infty }W$ is the slip velocity along the $z$ axis. These would furnish the bulk-scale boundary conditions for the outer layer velocities. However, the outer layer equations must also be complemented by appropriate boundary conditions along the axial directions. Fortunately, to determine the modified slip velocities we do not require any boundary conditions along the $X$ and $Z$ axes.

3.3 Solutions for the velocity field – Smoluchowski slip velocity

The inner layer equations of fluid motion are given by (3.9)–(3.12), subject to no-slip conditions at $Y=0$ and matching condition (3.14). The potential $\tilde{\unicode[STIX]{x1D711}}$ is given by (3.7). Although these equations are much simpler than the full set of governing equations, they are still nonlinear in nature and hence cannot be solved exactly. Approximate solutions can however, be derived assuming $De<O(1)$ . To this end, we expand all the flow variables in a regular asymptotic series of $De$ as follows (Leal Reference Leal1975; Chan & Leal Reference Chan and Leal1979; Bird et al. Reference Bird, Armstrong and Hassager1987):

(3.16) $$\begin{eqnarray}\text{Im}=\text{Im}_{0}+De\text{Im}_{1}+De^{2}\text{Im}_{2}+\cdots .\end{eqnarray}$$

Here we only deduce the asymptotic solutions in terms of $De$ , up to $O(De^{2})$ . Therefore, the slip velocity deduced will be correct up to $O(De^{2})$ .

3.3.1 Leading-order solutions

The leading-order equations in $De$ can be deduced from (3.9)–(3.12) quite easily. These remain same as those for Newtonian fluids. This is expected, since to leading order in $De$ , second-order fluids behave like Newtonian fluids. The leading-order equations read:

(3.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle -\frac{\unicode[STIX]{x2202}P_{0}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}^{2}U_{0}}{\unicode[STIX]{x2202}Y^{2}}+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}X}-\unicode[STIX]{x1D6FD}\cos \unicode[STIX]{x1D703}\right)=0;\quad -\frac{\unicode[STIX]{x2202}P_{0}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y}=0\\ \displaystyle -\frac{\unicode[STIX]{x2202}P_{0}}{\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}^{2}W_{0}}{\unicode[STIX]{x2202}Y^{2}}+\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Y^{2}}\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Z}-\unicode[STIX]{x1D6FD}\sin \unicode[STIX]{x1D703}\right)=0;\quad \frac{\unicode[STIX]{x2202}U_{0}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}V_{0}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}W_{0}}{\unicode[STIX]{x2202}Z}=0.\end{array}\right\}\end{eqnarray}$$

Equation (3.17) can be solved using a general procedure that we will also follow for higher-order corrections. First we integrate the $Y$ momentum equation to deduce the pressure, which here reads, $P_{0}=(\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}/\unicode[STIX]{x2202}Y^{2})/2=(1/2\bar{\unicode[STIX]{x1D701}}_{0}^{2})[\cosh (\bar{\unicode[STIX]{x1D701}}_{0}\tilde{\unicode[STIX]{x1D711}})-1]$ . This pressure can be separately plugged into the $X$ and $Z$ momentum equations and these can be subsequently integrated twice to obtain the $X$ and $Z$ velocities. These velocities read:

(3.18a,b ) $$\begin{eqnarray}U_{0}=\unicode[STIX]{x1D6FD}\cos (\unicode[STIX]{x1D703})\tilde{\unicode[STIX]{x1D711}}+\unicode[STIX]{x1D712}_{1}(X,Z)Y+\unicode[STIX]{x1D712}_{2}(X,Z);\quad W_{0}=\unicode[STIX]{x1D6FD}\sin (\unicode[STIX]{x1D703})\tilde{\unicode[STIX]{x1D711}}+\unicode[STIX]{x1D712}_{3}(X,Z)Y+\unicode[STIX]{x1D712}_{4}(X,Z).\end{eqnarray}$$

These velocities satisfy the no-slip condition at the wall, i.e. at $Y=0$ . Further, in the limit $Y\rightarrow \infty$ , $U_{0}$ and $W_{0}$ both must be finite and bounded. Therefore, we must have $\unicode[STIX]{x1D712}_{1}=\unicode[STIX]{x1D712}_{3}=0$ . Enforcing the no-slip boundary conditions we get, $\unicode[STIX]{x1D712}_{2}=-\unicode[STIX]{x1D6FD}\cos (\unicode[STIX]{x1D703})\unicode[STIX]{x1D701}(X,Z)$ ; $\unicode[STIX]{x1D712}_{4}=-\unicode[STIX]{x1D6FD}\sin (\unicode[STIX]{x1D703})\unicode[STIX]{x1D701}(X,Z)$ . The final form of the velocities is thus:

(3.19a,b ) $$\begin{eqnarray}U_{0}=-\unicode[STIX]{x1D6FD}\cos (\unicode[STIX]{x1D703})(\unicode[STIX]{x1D701}(X,Z)-\tilde{\unicode[STIX]{x1D711}});\quad W_{0}=-\unicode[STIX]{x1D6FD}\sin (\unicode[STIX]{x1D703})(\unicode[STIX]{x1D701}(X,Z)-\tilde{\unicode[STIX]{x1D711}}).\end{eqnarray}$$

The velocities in (3.19) can now be used in the continuity equation (the last of (3.17)) to deduce the $Y$ -velocity, subject to the condition that at $Y=0$ , $V_{0}=0$ . The $Y$ -velocity thus reads:

(3.20) $$\begin{eqnarray}V_{0}=\unicode[STIX]{x1D6FD}(\hat{\boldsymbol{e}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{1}\unicode[STIX]{x1D701})\left(Y-\int _{0}^{Y}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\,\text{d}\breve{Y}\right),\quad \text{where }\unicode[STIX]{x1D735}_{1}=\hat{\boldsymbol{e}}_{x}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X}+\hat{\boldsymbol{e}}_{z}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}Z}.\end{eqnarray}$$

The leading-order slip velocities in the $x$ and $z$ directions can be deduced respectively by evaluating the limits, $\lim _{Y\rightarrow \infty }U_{0}$ and $\lim _{Y\rightarrow \infty }W_{0}$ . The slip velocities are simply:

(3.21a,b ) $$\begin{eqnarray}u_{s}^{(0)}=\lim _{Y\rightarrow \infty }U_{0}=-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D701}(X,Z)\cos (\unicode[STIX]{x1D703});\quad w_{s}^{(0)}=\lim _{Y\rightarrow \infty }W_{0}=-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D701}(X,Z)\sin (\unicode[STIX]{x1D703}).\end{eqnarray}$$

In (3.21), $u_{s}^{(0)}$ and $w_{s}^{(0)}$ represent leading-order slip velocities along $x$ and $z$ directions respectively. This slip velocity can be represented in a compact vector form as follows:

(3.22) $$\begin{eqnarray}\boldsymbol{v}_{s}^{(0)}=u_{s}^{(0)}\boldsymbol{e}_{x}+w_{s}^{(0)}\boldsymbol{e}_{z}=-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D701}(X,Z)\hat{\boldsymbol{e}}.\end{eqnarray}$$

We recognize that (3.22) is simply the Smoluchowski slip velocity for Newtonian fluids. We further note that this condition remains valid even when the surface potential magnitude ( $\bar{\unicode[STIX]{x1D701}}_{0}$ ) is not small. Other conditions for its validity are $Pe\ll O(\unicode[STIX]{x1D6FF}^{-1})$ and $\unicode[STIX]{x1D6FD}\ll O(\unicode[STIX]{x1D6FF}^{-1})$ , as mentioned in the previous section. Note that to leading order, the slip velocity is exactly parallel to the applied electric field, i.e. there are no cross-flow components.

3.3.2 The $O(De)$ solutions

The $O(De)$ equations are given by:

(3.23) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle -\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}^{2}U_{1}}{\unicode[STIX]{x2202}Y^{2}}+\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xx}^{(0)}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xy}^{(0)}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xz}^{(0)}}{\unicode[STIX]{x2202}Z}\right)=0;\quad -\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{yy}^{(0)}}{\unicode[STIX]{x2202}Y}=0\\ \displaystyle -\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}^{2}W_{1}}{\unicode[STIX]{x2202}Y^{2}}+\left(\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{xz}^{(0)}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}\tilde{S}_{zy}^{(0)}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D61A}}_{zz}^{(0)}}{\unicode[STIX]{x2202}Z}\right)=0;\quad \frac{\unicode[STIX]{x2202}U_{1}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}V_{1}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}W_{1}}{\unicode[STIX]{x2202}Z}=0.\end{array}\right\}\end{eqnarray}$$

These are subjected to the no-slip condition at the wall and the velocities must be bounded as $Y\rightarrow \infty$ ; these conditions are identical to those used for the leading-order evaluations. The excess stress components, i.e.  $\tilde{\unicode[STIX]{x1D61A}}_{xx}$ , $\tilde{\unicode[STIX]{x1D61A}}_{xy}$ , $\tilde{\unicode[STIX]{x1D61A}}_{yz},\ldots$ etc. can be easily deduced from the relations given in (3.15). The solution procedure remains exactly the same as that described for the leading-order velocities, with $P_{1}=\tilde{\unicode[STIX]{x1D61A}}_{yy}^{(0)}$ . We note that, at $O(De)$ and at higher orders as well, we can continue without making any assumptions regarding the magnitude of $\bar{\unicode[STIX]{x1D701}}_{0}$ and evaluate the velocity fields based on (3.8). However, the resulting solutions thus obtained will have extremely complicated algebraic expressions and will in general be very hard to interpret. Therefore, the approximate analytical solutions for arbitrary $\bar{\unicode[STIX]{x1D701}}_{0}$ will be somewhat pyrrhic. Hence, to simplify the algebra, yet without sacrificing the essential physics, we now apply the low surface potential approximation, wherein, $\bar{\unicode[STIX]{x1D701}}_{0}\ll 1$ . Therefore, from here onwards we use (3.8) to represent the potential. The $X$ and $Z$ velocities are then given by:

(3.24) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle U_{1}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}_{x})F_{x}^{(1)}(Y)+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}_{z})F_{z}^{(1)}(Y)\\ \displaystyle W_{1}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}_{x})G_{x}^{(1)}(Y)+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FD}^{2}(\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}_{z})G_{z}^{(1)}(Y).\end{array}\right\}\end{eqnarray}$$

In (3.24), $\unicode[STIX]{x1D701}_{x}=\unicode[STIX]{x2202}\unicode[STIX]{x1D701}/\unicode[STIX]{x2202}X$ ; $\unicode[STIX]{x1D701}_{z}=\unicode[STIX]{x2202}\unicode[STIX]{x1D701}/\unicode[STIX]{x2202}Z$ . The functions $F$ and $G$ are given by:

(3.25) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}F_{x}^{(1)}=(2\cos ^{2}(\unicode[STIX]{x1D703})Y-2\cos ^{2}(\unicode[STIX]{x1D703})+b-b\cos ^{2}(\unicode[STIX]{x1D703}))\text{e}^{-Y}-b+(b+2)\cos ^{2}(\unicode[STIX]{x1D703})\\ F_{z}^{(1)}=(-b+2Y-2)\text{e}^{-Y}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})+(b+2)\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\\ G_{x}^{(1)}=(b+2-2Y)\text{e}^{-Y}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})-(b+2)\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\\ G_{z}^{(1)}=(2-2Y-b\cos ^{2}(\unicode[STIX]{x1D703})-2\cos ^{2}(\unicode[STIX]{x1D703})+2\cos ^{2}(\unicode[STIX]{x1D703})Y)\text{e}^{-Y}-2+(b+2)\cos ^{2}(\unicode[STIX]{x1D703}).\end{array}\right\}\end{eqnarray}$$

The $Y$ -velocity can be deduced by simply plugging the expressions for $W_{1}$ and $U_{1}$ , as mentioned in (3.24), into the continuity equation and integrating once, subject to $V_{1}(X,Y=0,Z)=0$ . However, the expression for $V_{1}$ is algebraically complicated and hence we do not explicitly mention it here. The resulting slip velocities along the surface are:

(3.26) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{s}^{(1)}=\lim _{Y\rightarrow \infty }U_{1}={\textstyle \frac{1}{4}}\unicode[STIX]{x1D701}\unicode[STIX]{x1D6FD}^{2}\{(-b+b\cos ^{2}(\unicode[STIX]{x1D703})+2\cos ^{2}(\unicode[STIX]{x1D703}))\unicode[STIX]{x1D701}_{x}+(b+2)\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\unicode[STIX]{x1D701}_{z}\}\\ \displaystyle w_{s}^{(1)}=\lim _{Y\rightarrow \infty }W_{1}=-{\textstyle \frac{1}{4}}\unicode[STIX]{x1D701}\unicode[STIX]{x1D6FD}^{2}\{(-b-2)\unicode[STIX]{x1D701}_{x}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})-2\unicode[STIX]{x1D701}_{z}+(b+2)\unicode[STIX]{x1D701}_{z}\cos ^{2}(\unicode[STIX]{x1D703})\}.\end{array}\right\}\end{eqnarray}$$

In vector form the slip velocity is thus,

(3.27) $$\begin{eqnarray}\boldsymbol{v}_{s}^{(1)}=u_{s}^{(1)}\boldsymbol{e}_{x}+w_{s}^{(1)}\boldsymbol{e}_{z}.\end{eqnarray}$$

We discuss further the nature of the slip velocity in the forthcoming section (§ 4).

3.3.3 The $O(De^{2})$ solutions

The $O(De^{2})$ equations can be deduced in a similar fashion as done for the $O(De)$ solutions. These velocities are subject the same conditions as mentioned for the previous orders. The solution procedure also remains the same as described for the $O(1)$ and $O(De)$ solutions. The solutions for $U_{2}$ and $W_{2}$ thus obtained have complex algebraic forms and thus we choose not to mention them in this document for the sake of brevity. The MAPLE files containing the expressions can be made available upon request to the authors. However, we do mention the slip velocities occurring at $O(De^{2})$ . These are given by:

(3.28) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{s}^{(2)}=\lim _{Y\rightarrow \infty }U_{2}={\textstyle \frac{1}{64}}\unicode[STIX]{x1D6FD}^{3}\unicode[STIX]{x1D701}[F_{xx}^{(2)}\unicode[STIX]{x1D701}_{xx}+F_{zz}^{(2)}\unicode[STIX]{x1D701}_{zz}+F_{xz}^{(2)}\unicode[STIX]{x1D701}_{xz}+M_{x}\unicode[STIX]{x1D701}_{x}^{2}+M_{z}\unicode[STIX]{x1D701}_{z}^{2}+M_{xz}\unicode[STIX]{x1D701}_{x}\unicode[STIX]{x1D701}_{z}]\\ \displaystyle w_{s}^{(2)}=\lim _{Y\rightarrow \infty }W_{2}={\textstyle \frac{1}{64}}\unicode[STIX]{x1D701}\unicode[STIX]{x1D6FD}^{3}[G_{xx}^{(2)}\unicode[STIX]{x1D701}_{xx}+G_{zz}^{(2)}\unicode[STIX]{x1D701}_{zz}+G_{xz}^{(2)}\unicode[STIX]{x1D701}_{xz}+N_{x}\unicode[STIX]{x1D701}_{x}^{2}+N_{z}\unicode[STIX]{x1D701}_{z}^{2}+N_{xz}\unicode[STIX]{x1D701}_{x}\unicode[STIX]{x1D701}_{z}].\end{array}\right\}\end{eqnarray}$$

The functions $F$ , $G$ , $M$ , $N$ etc. have the following expressions:

(3.29) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}F_{xx}^{(2)}=-[24\cos ^{3}(\unicode[STIX]{x1D703})-4b^{2}\cos (\unicode[STIX]{x1D703})+18\cos ^{3}(\unicode[STIX]{x1D703})b-18b\cos (\unicode[STIX]{x1D703})+4b^{2}\cos ^{3}(\unicode[STIX]{x1D703})]\unicode[STIX]{x1D701}\\ F_{zz}^{(2)}=-[24\cos (\unicode[STIX]{x1D703})-24\cos ^{3}(\unicode[STIX]{x1D703})-b^{2}\cos (\unicode[STIX]{x1D703})+6b\cos (\unicode[STIX]{x1D703})-4\cos ^{3}(\unicode[STIX]{x1D703})b^{2}-18\cos ^{3}(\unicode[STIX]{x1D703})b]\unicode[STIX]{x1D701}\\ F_{xz}^{(2)}=-[(36b+8b^{2}+48)\cos ^{2}(\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})+(b^{2}-6b)\sin (\unicode[STIX]{x1D703})]\unicode[STIX]{x1D701}\\ M_{x}=4[-8b\cos (\unicode[STIX]{x1D703})+2b^{2}\cos ^{3}(\unicode[STIX]{x1D703})+8\cos ^{3}(\unicode[STIX]{x1D703})b-2b^{2}\cos (\unicode[STIX]{x1D703})+8\cos ^{3}(\unicode[STIX]{x1D703})]\\ M_{z}=4[-8\cos ^{3}(\unicode[STIX]{x1D703})b+b^{2}\cos (\unicode[STIX]{x1D703})-2b^{2}\cos ^{3}(\unicode[STIX]{x1D703})+6b\cos (\unicode[STIX]{x1D703})-8\cos ^{3}(\unicode[STIX]{x1D703})+8\cos (\unicode[STIX]{x1D703})]\\ M_{xz}=4\sin (\unicode[STIX]{x1D703})[16b\cos ^{2}(\unicode[STIX]{x1D703})+4\cos ^{2}(\unicode[STIX]{x1D703})b^{2}-6b+16\cos ^{2}(\unicode[STIX]{x1D703})-b^{2}]\\ G_{xx}^{(2)}=-[(24+18b+4b^{2})\cos ^{2}(\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})-(12b+5b^{2})\sin (\unicode[STIX]{x1D703})]\unicode[STIX]{x1D701}\\ G_{zz}^{(2)}=-[24\sin (\unicode[STIX]{x1D703})-(24+18b+4b^{2})\cos (\unicode[STIX]{x1D703})^{2}\sin (\unicode[STIX]{x1D703})]\unicode[STIX]{x1D701}\\ G_{xz}^{(2)}=-[(9b^{2}+30b+48)\cos (\unicode[STIX]{x1D703})-(8b^{2}+36b+48)\cos ^{3}(\unicode[STIX]{x1D703})]\unicode[STIX]{x1D701}\\ N_{x}=-4[(2b^{2}+8b+8)\sin (\unicode[STIX]{x1D703})\cos ^{2}(\unicode[STIX]{x1D703})-(2b+b^{2})\sin (\unicode[STIX]{x1D703})]\\ N_{z}=-4[(-8b-2b^{2}-8)\sin (\unicode[STIX]{x1D703})\cos ^{2}(\unicode[STIX]{x1D703})+8\sin (\unicode[STIX]{x1D703})]\\ N_{xz}=4\cos (\unicode[STIX]{x1D703})[-10b+4\cos ^{2}(\unicode[STIX]{x1D703})b^{2}-3b^{2}+16b\cos ^{2}(\unicode[STIX]{x1D703})+16\cos ^{2}(\unicode[STIX]{x1D703})-16].\end{array}\right\}\end{eqnarray}$$

In vector form, the slip velocity reads,

(3.30) $$\begin{eqnarray}\boldsymbol{v}_{s}^{(2)}=u_{s}^{(2)}\boldsymbol{e}_{x}+w_{s}^{(2)}\boldsymbol{e}_{z}.\end{eqnarray}$$

Therefore, the total slip velocity, correct up to $O(De^{2})$ , is given by:

(3.31) $$\begin{eqnarray}\boldsymbol{v}_{s}=\boldsymbol{v}_{s}^{(0)}+De\boldsymbol{v}_{s}^{(1)}+De^{2}\boldsymbol{v}_{s}^{(2)}+\cdots .\end{eqnarray}$$

Based on (3.31), the outer layer velocities (as mentioned in (3.6)) can be expanded in an asymptotic series of $De$ as, $\boldsymbol{v}=\boldsymbol{v}_{0}+De\boldsymbol{v}_{1}+De^{2}\boldsymbol{v}_{2}+\cdots$ . The $O(De^{k})$ velocity field satisfies the equations, $-\unicode[STIX]{x1D735}p_{k}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{k}+(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{(k-1)})=0$ ; $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{k}=0$ , subject to the following conditions: $\boldsymbol{v}_{k}(x,y=0,z)=\boldsymbol{v}_{s}^{(k)}(x,z)$ ; $\boldsymbol{v}_{k}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{y}=0$ . Of course these asymptotic solutions are only valid for $De<O(1)$ . Note that, in (3.28), $\unicode[STIX]{x1D701}_{xx}=\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D701}/\unicode[STIX]{x2202}X^{2}$ ; $\unicode[STIX]{x1D701}_{zz}=\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D701}/\unicode[STIX]{x2202}Z^{2}$ ; $\unicode[STIX]{x1D701}_{xz}=\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D701}/\unicode[STIX]{x2202}X\unicode[STIX]{x2202}Z$ .

5.1 Example 1: two-dimensional axially periodic potential

Consider an axially periodic potential of the form:

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D701}^{\prime }(x^{\prime })=\unicode[STIX]{x1D701}_{0}(n+m\cos (2\unicode[STIX]{x03C0}x^{\prime }/L)).\end{eqnarray}$$

Here, the natural length scale of variation in the surface potential can be chosen as $d=L/2\unicode[STIX]{x03C0}$ . Therefore, application of the thin EDL approximation and evaluation of the subsequent slip velocity would require, $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D706}_{D}/d\ll 1$ . Since the flow under consideration here is two-dimensional, the slip velocity is given by (4.2). Plugging (5.1) into (4.2), we obtain:

(5.2) $$\begin{eqnarray}\displaystyle u_{s} & = & \displaystyle -\unicode[STIX]{x1D6FD}(n+m\cos (x))+De\left(-{\textstyle \frac{1}{4}}m^{2}\unicode[STIX]{x1D6FD}^{2}\sin (2x)-{\textstyle \frac{1}{2}}n\unicode[STIX]{x1D6FD}^{2}m\sin (x)\right)\nonumber\\ \displaystyle & & \displaystyle +\,De^{2}\unicode[STIX]{x1D6FD}^{3}\left({\textstyle \frac{5}{8}}nm^{2}\cos (2x)+{\textstyle \frac{5}{32}}m^{3}\cos (x)+{\textstyle \frac{3}{8}}m\cos (x)n^{2}+{\textstyle \frac{1}{8}}nm^{2}+{\textstyle \frac{7}{32}}m^{3}\cos (3x)\right)\nonumber\\ \displaystyle & & \displaystyle +\,\cdots .\end{eqnarray}$$

The most interesting point to note from (5.2) is that, at $O(De^{2})$ , there is an axially invariant component of the slip velocity given by $(1/8)nm^{2}$ . This component would alter the net throughput in a channel. For Newtonian fluids, the net throughput would only depend on $n$ ; however, for second-order fluids the net throughput also depends on $m$ , as evident from (5.2).

Note that, if $m=0$ , i.e. if the potential is uniform along the channel, the $O(De)$ , $O(De^{2})$ and consequently all higher-order contributions to the slip velocity vanish. In such cases the Smoluchowski slip velocities for second-order fluids and Newtonian fluids are identical. Figure 2 depicts the typical variation in $u_{s}$ as a function of $x$ for four different values of $De=0$ , 0.1, 0.2 and 0.3, while we have chosen, $n=0.5$ and $m=1$ . We observe that notable deviation from the Newtonian case ( $De=0$ ) only occurs where the potential has comparatively larger values. As we increase the value of $De$ , the slip velocity becomes somewhat asymmetric with respect to  $x$ .

Figure 2. Variation in slip velocity ( $u_{s}$ ) versus $x$ , for four different values of $De=0$ (Newtonian), 0.1, 0.2 and 0.3. The potential here is periodically varying along the wall, with the expression, $\unicode[STIX]{x1D701}^{\prime }(x^{\prime })=\unicode[STIX]{x1D701}_{0}[n+m\cos (2\unicode[STIX]{x03C0}x^{\prime }/L)]$ . Here we have taken $n=0.5$ , $m=1$ and $\unicode[STIX]{x1D6FD}=1.5$ .

5.1.1 Comparison with numerical simulations

In an effort to verify the accuracy of the asymptotic approximations, we compare them with independent numerical simulations. To this end, we simulate purely two-dimensional flows where the surface potential is chosen as (5.1). Figure 3 depicts the whole set-up for the numerical simulations, including the governing equations, the boundary conditions and the mesh. The details of the domain size and number of grid points are stated in the caption. The mesh size has been made smaller near the bottom wall to capture the dynamics within the EDL. The governing equations, as shown in figure 3, are subject to the following assumptions: (i) low surface potential; (ii) low ionic Péclet number ( $Pe$ ) and (iii) the field strength $\unicode[STIX]{x1D6FD}\sim O(1)$ or less. Note that here the potential can be expressed as follows: $\unicode[STIX]{x1D711}=ne^{-y/\unicode[STIX]{x1D6FF}}+me^{-y/\unicode[STIX]{x1D6FF}_{1}}\cos (x)$ , where $\unicode[STIX]{x1D6FF}_{1}^{-1}=\sqrt{1+\unicode[STIX]{x1D6FF}^{-2}}$ . We executed the numerical simulations for the thin EDL limit, implying $\unicode[STIX]{x1D6FF}\ll 1$ , in order to comply with the conditions of our asymptotic analysis. The simulations were performed in ANSYS FLUENT, which uses a finite volume based framework along with the SIMPLE algorithm for pressure–velocity coupling (Patankar Reference Patankar1980). A more detailed formulation for the numerical simulations along with some of the major difficulties in simulating the given system are discussed in appendix B.

Figure 3. Set-up for the numerical computation. The governing equations and the boundary conditions are presented in the set-up for comprehensiveness. The domain size is chosen as $2\unicode[STIX]{x03C0}\times 10$ ensuring lengthwise periodicity and an infinite domain along the transverse direction, as required in our analysis. We perform simulations with grid resolution $50\times 600$ . We consider non-uniformly spaced grids, and the successive ratios between two consecutive grid sizes are 1.05. The arrow shows the direction of non-uniform spacing. Note that the expressions for excess polymeric stresses $\unicode[STIX]{x1D64E}$ have been given in § 2. Here, $\unicode[STIX]{x1D6FF}=0.01$ . Other relevant details have been mentioned within the figure itself.

Figure 4 demonstrates the comparison between numerical and analytical solutions for the slip velocity for four different values of Deborah number ( $De=0.04$ , 0.08, 0.11, 0.15). Specifically, we have compared the quantity, $\unicode[STIX]{x0394}u_{s}=u_{s}-u_{newt.}$ , i.e. the difference between the slip velocities of Newtonian and second-order fluids, with $\unicode[STIX]{x1D6FF}=0.01$ . Note that, analytically $\unicode[STIX]{x0394}u_{s}=De\,u_{s}^{(1)}+De^{2}u_{s}^{(2)}+\cdots$ and hence we essentially compare only the contribution from the higher-order corrections with the corresponding numerical estimates. Therefore, the present comparison is perhaps more rigorous in nature. For the numerical solutions, the slip velocity was estimated just outside the EDL, at a location $y=0.013$ . Values of all other relevant parameters have been stated in the caption. We note that very good comparison between numerical and analytical solutions can be observed for all values of $De$ , which validates our analytical estimates.

Figure 4. Comparison between numerical (marker) and analytical (lines) solutions for a periodic zeta potential for four different values of $De=0.04$ , 0.08, 0.11 and 0.15. We have compared the solutions for the quantity, $\unicode[STIX]{x0394}u_{s}=u_{s}-u_{newt.}$ , i.e. the change in the slip velocity for second-order fluids. Here, $u_{s}$ is the slip velocity for second-order fluids and $u_{newt.}$ is that of Newtonian fluids. Values of other parameters are as follows: $\unicode[STIX]{x1D6FD}=0.5$ , $m=0.2$ , $n=0.5$ , $\unicode[STIX]{x1D6FF}=0.01$ . The numerical solutions have been depicted for the location $y=0.013$ , i.e. just outside the EDL.

Figure 5 demonstrates the novel scaling pattern of normal stresses (as mentioned in § 3.3) inside the EDL for second-order fluids. Here, we plot the numerically evaluated normal stress component ( $\unicode[STIX]{x1D70F}_{xx}$ ) as a function of EDL thickness ( $\unicode[STIX]{x1D6FF}$ ) for both Newtonian and second-order fluids (with $De=0.15$ ). Simultaneously, the predicted scaling patterns have also been plotted for both of the fluids. Values of all other relevant parameters have been stated in the caption. Our scaling indicates that for second-order fluids, $\unicode[STIX]{x1D70F}_{xx}\sim T_{s}\unicode[STIX]{x1D6FF}^{-2}$ , where $T_{s}$ depends on $n$ , $m$ , $\unicode[STIX]{x1D6FD}$ etc. and for Newtonian fluids, $\unicode[STIX]{x1D70F}_{xx}\sim T_{n}=\text{a constant}$ , depending on $n$ , $m$ , $\unicode[STIX]{x1D6FD}$ etc. The figure under consideration clearly verifies the said scaling and acts as another good test for the validity of the asymptotic analysis executed herein. The constants have been chosen as $T_{s}=0.16$ and $T_{n}=0.1$ .

Figure 5. Comparison between analytical scaling (lines) and numerical solutions (markers) for the normal stress ( $\unicode[STIX]{x1D70F}_{xx}$ ) inside the EDL as a function of $\unicode[STIX]{x1D6FF}$ . Our analytical scaling suggests that for second-order fluids, $\unicode[STIX]{x1D70F}_{xx}\sim \unicode[STIX]{x1D6FF}^{-2}$ in the EDL, whereas for a Newtonian fluid, $\unicode[STIX]{x1D70F}_{xx}\sim \unicode[STIX]{x1D6FF}^{0}$ , i.e. has a constant value inside the EDL. The exact same trend is indeed reflected by the numerical simulations. The rest of the parameters have been taken as follows: $\unicode[STIX]{x1D6FD}=0.5$ , $n=0.5$ , $m=0.2$ .

5.2 Example 2: step-change-like variation in zeta potential

A step change in the zeta potential is typically represented by an abrupt discontinuous change in $\unicode[STIX]{x1D701}$ along the wall. A number of studies have previously been executed (Fu Reference Fu2003; Brotherton & Davis Reference Brotherton and Davis2004; Yariv Reference Yariv2004), which investigate the dynamics of Newtonian fluids in the presence of a step change in the surface potential. However, in reality, such abrupt changes are not likely to exist, since they would give rise to very high electric field strengths inside the solid surface. Therefore, a more likely scenario is where the potential varies smoothly over a short distance as compared to the total length of the device. Hence, if our resolution lies in the length scale of the channel/device, the sharp but continuous change in the potential would indeed be perceived as an abrupt step change.

A good analogy for modelling abrupt variations in a certain variable can be found in phase-field methods (Jacqmin Reference Jacqmin1999), where the order parameter $\unicode[STIX]{x1D719}$ changes abruptly across the interface from $-1$ to 1. Although on the system length scale this is a step change in the value of $\unicode[STIX]{x1D719}$ leading to a sharp interface, on a finer scale the interface is actually assumed to be diffuse (Yang, Li & Ding Reference Yang, Li and Ding2013) with a typical thickness $\unicode[STIX]{x1D70D}$ . Therefore, the variation in $\unicode[STIX]{x1D719}$ is actually assumed to be smooth and is taken to be of the form (Wang, Qian & Sheng Reference Wang, Qian and Sheng2008), $\unicode[STIX]{x1D719}=\tanh [(x^{\prime }-x_{0}^{\prime })/\unicode[STIX]{x1D70D}]$ , where $x_{0}^{\prime }$ is the location of the interface. The system length scale (say  $L$ ) is much larger than the typical interface thickness $\unicode[STIX]{x1D70D}$ , with the ratio of the two being termed the Cahn number (Chaudhury, Ghosh & Chakraborty Reference Chaudhury, Ghosh and Chakraborty2013). In the same spirit, here also we take a similar variation in the surface potential, assuming the potential changes from $\unicode[STIX]{x1D701}_{1}^{\prime }$ to $\unicode[STIX]{x1D701}_{2}^{\prime }$ smoothly over a length scale of $d$ . This length scale is assumed to be much smaller than the overall system length scale  $L$ . However, since typical EDL thickness is very small, usually of the order of nanometres, the length scale $d$ can still be large compared to $\unicode[STIX]{x1D706}_{D}$ . Therefore, we assume that although $d\ll L$ , we still have $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D706}_{D}/d\ll 1$ , so that the thin EDL approximation remains valid. Accordingly, the zeta potential is expressed as:

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D701}^{\prime }(x^{\prime })={\textstyle \frac{1}{2}}[\unicode[STIX]{x1D701}_{1}^{\prime }+\unicode[STIX]{x1D701}_{2}^{\prime }+(\unicode[STIX]{x1D701}_{2}^{\prime }-\unicode[STIX]{x1D701}_{1}^{\prime })\tanh (x^{\prime }/d)].\end{eqnarray}$$

We can choose $\unicode[STIX]{x1D701}_{0}=\unicode[STIX]{x1D701}_{1}^{\prime }$ and $d$ as the length scale of variation in the potential. The non-dimensional potential would thus read:

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D701}(x)={\textstyle \frac{1}{2}}[1+\unicode[STIX]{x1D701}_{2}+(\unicode[STIX]{x1D701}_{2}-1)\tanh (x)].\end{eqnarray}$$

We further note that, similar continuous but sharp variations are also witnessed in shock waves (Kundu & Cohen Reference Kundu and Cohen2004), where the fluid properties vary continuously yet extremely rapidly over a very small length (thickness of the shock wave front) so that the change in the properties seems abrupt over a larger length scale.

Equation (5.4) can be plugged into (4.2) to obtain the slip velocity correct up to $O(De^{2})$ . Figure 6 showcases the resulting variation in $u_{s}$ versus $x$ for three different values of $De=0$ (Newtonian), 0.2 and 0.35, while the other parameters are stated in the caption. In the inset we have depicted the variation in the surface potential as a function of axial location to represent the step-change-like variation. We again note that, similar to the previous figure, the deviations from Newtonian behaviour are only notable when the potential, and hence the velocity, have relatively high magnitude. Note that the slip velocities for second-order fluids slightly lag behind those of Newtonian fluids, as attributable to the fading memory of such fluids.

Figure 6. Variation in slip velocity ( $u_{s}$ ) versus $x$ for four different values of $De=0$ (Newtonian), 0.2 and 0.35. The potential here mimics a step-change-like variation, given by, $\unicode[STIX]{x1D701}^{\prime }(x^{\prime })=1/2\cdot ((\unicode[STIX]{x1D701}_{1}^{\prime }+\unicode[STIX]{x1D701}_{2}^{\prime })+(\unicode[STIX]{x1D701}_{2}^{\prime }-\unicode[STIX]{x1D701}_{1}^{\prime })\cdot \tanh (x^{\prime }/d))$ . This kind of variation resembles that of a step change in the potential, from $\unicode[STIX]{x1D701}_{1}$ to $\unicode[STIX]{x1D701}_{2}$ . In the inset we have shown how the potential varies along the surface. Here, we have taken $\unicode[STIX]{x1D701}_{2}=0.2$ and $\unicode[STIX]{x1D6FD}=1.75$ .

5.3 Cases of transverse flow

We take up one last example, where cases of fully three-dimensional transverse flows are investigated. To this end, we would consider a zeta potential of the following form (non-dimensional):

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D701}(x,z)=[n+m\cos (x)\cos (z)].\end{eqnarray}$$

Clearly, this potential varies along both the $x$ and $z$ axes and is analogous to a ‘chess-board’ patterning. We also consider an applied electric field which has components along both the $x$ and $z$ axes. The choices for the characteristic scales here remain the same as those in § 5.1.

Figure 7. Vector plot of slip velocity ( $\boldsymbol{v}_{s}=u_{s}\hat{\boldsymbol{e}}_{\boldsymbol{x}}+w_{s}\hat{\boldsymbol{e}}_{\boldsymbol{z}}$ ) on the $xz$ plane for an applied electric field at an angle $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ with the $x$ axis, denoted by the red arrow. Here, (a) is for Newtonian fluids and (b) is for second-order fluids, with $De=0.35$ , $b=1$ . Other parameters have the values: $\unicode[STIX]{x1D6FD}=3$ , $m=1$ , $n=0$ . The quantities in the vector plot have been scaled with a value 1.8.

Figure 7(a,b) depict the vector plot of the slip velocity ( $\boldsymbol{v}_{s}=u_{s}\hat{\boldsymbol{e}}_{\boldsymbol{x}}+w_{s}\hat{\boldsymbol{e}}_{\boldsymbol{z}}$ ) for Newtonian (in a) and second-order fluids (in b), with the applied field having components in both the $x$ and $z$ directions. In this figure the field (represented with the thick red arrow) acts at an angle $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ with the $x$ axis. Values of all other parameters have been stated in the caption. We first note from figure 7(a) that for Newtonian fluids, the slip velocity ( $\boldsymbol{v}_{s}$ ) is perfectly aligned to, or parallel with, the applied electric field. On the other hand, figure 7(b) clearly demonstrates that for second-order fluids, the velocity is, in general, non-parallel to the applied electric field, as we discussed in § 4 (see (4.1)). Note that in a Newtonian fluid, such anisotropic effects can only come into play if bulk-scale concentration gradients are present (chemioosmotic flow (Khair & Squires Reference Khair and Squires2008; Yariv Reference Yariv2009)). However, in sharp contrast, for second-order fluids, the nonlinearity in the constitutive equations itself guarantees such ‘cross-stream’ flows, without any bulk-scale concentration gradients. The main reason for the said misalignment is the presence of secondary flows originating from $O(De)$ and $O(De^{2})$ terms, as discussed in detail in §§ 3 and 4. On a related note, we emphasize that a series of studies by Vinogradova and coworkers (Belyaev & Vinogradova Reference Belyaev and Vinogradova2010; Feuillebois, Bazant & Vinogradova Reference Feuillebois, Bazant and Vinogradova2010; Belyaev & Vinogradova Reference Belyaev and Vinogradova2011; Vinogradova & Belyaev Reference Vinogradova and Belyaev2011; Schmieschek et al. Reference Schmieschek2012) and others (Ghosh & Chakraborty Reference Ghosh and Chakraborty2013) have demonstrated that similar anisotropic flows for Newotnian fluids might also be generated by employing patterned slip on the channel wall in the presence of simple pressure-driven flows.