2. Formulation of the basic flow and solute transport problem

2.1. Flow field and basic problem description

We study the flow and dispersion of a neutral (uncharged) solute subject to simultaneous PDF and EOF in a cylindrical tube of inner radius $a\,.$ Figure 1 shows schematics of an initial condition of the solute and example flow profiles. The flow is generated from applied axial pressure differentials (resulting in a Poiseuille-type flow component, see, e.g. Langlois & Deville Reference Langlois and Deville2014) and/or by EOF caused by the presence of a finite sized EDL at the wetted wall and an electric field applied along the axis of the tube. The EDL may or may not overlap with itself at the centre of the tube (see Appendix A). We will denote the characteristic axial width of the solute zone as $\sigma _{x}$ and assume the solute is of uniform molecular diffusivity ( $D$ ) throughout the tube. The Debye length is denoted as $\lambda _{D}$ . Let $r^{*}=r/a$ denote the dimensionless radial coordinate and $\phi =a/\lambda _{D}$ be the tube radius scaled by Debye length. For EOF, we consider a general electrolyte and a wall zeta potential sufficiently weak such that the Debye–Hückel approximation holds (Probstein Reference Probstein1994). We will further neglect the effects of streaming potential on the fluid velocity field. Note the latter effects should not create new shapes of the velocity field but may introduce a coupling between EOF and PDF which is not explicitly treated here. We provide a derivation of the electric potential field and EOF velocity profile, valid for an EDL of arbitrary thickness relative to the channel radius, in Appendix A. From (A11), the electric potential field within the tube is given by $\psi (r^{*})=\zeta I_{0}(\phi r^{*})/I_{0}(\phi )$ , where $I_{0}(z)$ is a zeroth-order modified Bessel function of the first kind and $\zeta$ is the zeta potential. Further, from (A14), the EOF velocity profile can be expressed as $u_{e}(r^{*})=u_{\textit{HS}}(1-\psi (r^{*})/\zeta )$ where $u_{\textit{HS}}=-\varepsilon _{e}E\zeta /\mu$ is the Helmholtz–Smoluchowski velocity scale. Here, $\varepsilon _{e}$ is the permittivity of the fluid (assumed to be uniform throughout the tube), $E$ is the magnitude of the applied electric field and $\mu$ is the dynamic viscosity of the fluid. For a compact notation, we define the following dimensionless cross-sectional average of the radial electric potential:

(2.1) \begin{equation} \eta =\int \limits_{0}^{1}2r^{*}\frac{\psi\!\left(r^{*}\right)}{\zeta }\mathrm{d}r^{*}=\frac{2}{\phi }\frac{I_{1}\!\left(\phi \right)}{I_{0}\!\left(\phi \right)} .\end{equation}

Taking the area average of the EOF profile

(2.2) \begin{equation} \left\langle u_{e}\right\rangle =u_{\textit{HS}}\!\left(1-\int\limits _{0}^{1}2r^{*}\frac{\psi\!\left(r^{*}\right)}{\zeta }\mathrm{d}r^{*}\right)=u_{\textit{HS}}\!\left(1-\eta \right) .\end{equation}

Further, for the case of a thick EDL relative to the tube radius ( $\phi \ll 1$ ), we can expand $I_{0}$ using a Taylor series expansion. Neglecting terms of order $O[\phi ^{4}]$ and higher, we can approximate

(2.3) \begin{equation} u_{e} (r^{*})\approx u_{\textit{HS}}\frac{\phi ^{2}}{4} \big(1-r^{*2}\big) .\end{equation}

Thus, as the EDL becomes thicker relative to the radius of the channel, EOF tends to a parabolic profile, and $\phi$ informs the bulk velocity of the flow.

The PDF profile is given as $u_{p}(r^{*})=2\langle u_{p}\rangle (1-r^{*2})$ with $\langle u_{p}\rangle$ denoting the area-averaged PDF profile. A tube subject to both pressure gradients and EOF has a net velocity profile given by $u(r^{*})=u_{e}(r^{*})+u_{p}(r^{*})$ . Combining with (2.2), we can express the full velocity profile as

(2.4) \begin{equation} u(r^{*})=2\langle u_{p}\rangle \big(1-r^{*2}\big)+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(1-\frac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right) .\end{equation}

The area average of (2.4) is simply $\langle u\rangle =\langle u_{e}\rangle +\langle u_{p}\rangle$ . Lastly, we take the deviation from this area average as

(2.5) \begin{equation} u^{\prime} (r^{*})=u(r^{*})- \langle u \rangle = \langle u_{p} \rangle \big(1-2r^{*2}\big)+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta -\frac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right) .\end{equation}

Figure 1(a) shows a flow profile where PDF and EOF are in the same direction, corresponding to values of $\langle u_{p}\rangle \gt 0$ and $\langle u_{e}\rangle \gt 0$ . Figure 1(b) shows a flow profile where EOF is opposed by PDF, corresponding to values of $\langle u_{p}\rangle \lt 0\lt \langle u_{e}\rangle$ .

Figure 2 shows examples of the well-known axial flow profiles of combined PDF and EOF in coordinates of $u(r^{*})$ versus the dimensionless inner radius, $r^{*}$ . We show velocity profiles for three combinations of $\langle u_{p}\rangle$ and $\langle u_{e}\rangle$ (three line colours) such that $\langle u_{e}\rangle +\langle u_{p}\rangle =1$ . For each combination of $\langle u_{p}\rangle$ and $\langle u_{e}\rangle$ , we present profiles for $\phi$ values of 10 and 30 (different line types). We shall later show and discuss EOF profiles where $\langle u_{p}\rangle \lt 0\lt \langle u_{e}\rangle$ which are useful for minimising the variance of the solute zone. That is, PDF in opposition to bulk EOF can be used to reduce dispersion for a fixed value of overall bulk velocity.

Figure 2.

Example axial flow profiles for three combinations of PDF and EOF, and two values of the tube radius scaled by the Debye length, $\phi$ . The ordinate shows dimensionless radius, $r^{*}$ , and the abscissa is the flow velocity, $u(r^{*})$ . The three line colours correspond to cases of EOF in opposition to bulk PDF, EOF and PDF in the same direction and PDF in opposition to bulk EOF. For each combination of EOF and PDF, profiles are shown for $\phi$ values of 10 and 30 (different line types). All curves were produced with unit bulk velocity.

2.2. Convective diffusion of solute

We will consider the advective diffusion of a neutral solute subject to the aforementioned velocity profiles. The governing equation, boundary conditions and initial condition for the concentration $c(r,x,t)$ of a solute are

(2.6) \begin{align} &\frac{\partial c}{\partial t}+u\!\left(r\right)\frac{\partial c}{\partial x}=D\left(\frac{1}{r}\frac{\partial }{\partial r}\!\left(r\frac{\partial c}{\partial r}\right)+\frac{\partial ^{2}c}{\partial x^{2}}\right); \left.\frac{\partial c}{\partial r}\right| _{r\,=\,a,0}=0,\nonumber\\ &x^{n}\frac{\partial ^{m}c}{\partial x^{m}}\rightarrow 0\text{ as }\!\left| x\right| \rightarrow \infty \text{ for }n, m\geq 0,c\!\left(r,x,t\,=\,0\right) = c\!\left(r,x,0\right). \end{align}

Here, the boundary condition at $r=0$ is a result of the assumed azimuthal symmetry and $c(r,x,0)$ denotes the initial solute distribution. We shall consider two analytical solution methods for this equation in §§ 3 and 4.

3. Taylor dispersion solution of the quasi-steady regime

We here present a derivation of the effective Taylor dispersion (Taylor Reference Taylor1953) for combined PDF and EOF. We will use an area-averaging and deviation formulation (and notation) similar to that of Stone & Brenner (Reference Stone and Brenner1999). Accordingly, we seek to derive the two-dimensional solute concentration field in regimes where the observation time scale is large with respect to the characteristic time of transverse diffusion, for flows described in § 2.1. We expand the concentration and flow fields into cross-sectional averages, defined as $\langle (.)\rangle = (1/\pi a^{2})\int _{0}^{a}2\pi r(.)\mathrm{d}r$ , and deviations therefrom, $(.)^{\prime}$ , such that $(.)^{\prime} = (.)-\langle (.)\rangle$ . First, we expand (2.6) as

(3.1) \begin{align} &\frac{\partial \!\left\langle c\right\rangle }{\partial t}+\frac{\partial c^{\prime}}{\partial t}+\left\langle u\right\rangle \frac{\partial \!\left\langle c\right\rangle }{\partial x}+u^{\prime}\!\left(r\right)\frac{\partial \!\left\langle c\right\rangle }{\partial x}+\left\langle u\right\rangle \frac{\partial c^{\prime}}{\partial x}+u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}=\nonumber\\ &D\left(\frac{1}{r}\frac{\partial }{\partial r}\!\left(r\frac{\partial c^{\prime}}{\partial r}\right)+\frac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}+\frac{\partial ^{2}c^{\prime}}{\partial x^{2}}\right);\left.\frac{\partial c^{\prime}}{\partial r}\right| _{r\, =\, a,0}=0. \end{align}

Taking the area average of this equation

(3.2) \begin{equation} \frac{\partial \!\left\langle c\right\rangle }{\partial t}+\left\langle u\right\rangle \frac{\partial \!\left\langle c\right\rangle }{\partial x}=D\left(\frac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\right)-\left\langle u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}\right\rangle .\end{equation}

Subtracting (3.2) from (3.1), we derive an expression for the deviation variable $c^{\prime}$

(3.3) \begin{equation} \frac{\partial c^{\prime}}{\partial t}+u^{\prime}\!\left(r\right)\frac{\partial \!\left\langle c\right\rangle }{\partial x}+\left\langle u\right\rangle \frac{\partial c^{\prime}}{\partial x}+u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}=D\left(\frac{1}{r}\frac{\partial }{\partial r}\!\left(r\frac{\partial c^{\prime}}{\partial r}\right)+\frac{\partial ^{2}c^{\prime}}{\partial x^{2}}\right)+\left\langle u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}\right\rangle .\end{equation}

The latter equation is again subject to the condition that $\partial c^{\prime}/\partial r| _{r=a,0}=0$ . Let $c_{o}$ denote the characteristic scale (magnitude) of the area-averaged concentration of solute and $c_{o}{}^{\prime}$ denote the (much smaller) characteristic scale for the deviation from this average. We non-dimensionalise (3.3) using

(3.4) \begin{equation} r^{*}=\frac{r}{a} ;\, x_{\sigma^{*}}=\frac{x}{\sigma _{x}} ;\, t^{*}=\frac{t\!\left\langle u\right\rangle }{\sigma _{x}} ;\, {\hat u^{\prime}}=\frac{u^{\prime}\!\left(r\right)}{\left\langle u\right\rangle } ;\, {\left\langle \hat u\right\rangle }=\frac{\left\langle u\right\rangle }{\left\langle u\right\rangle } ;\, c^{\prime*}=\frac{c^{\prime}}{c_{o}{}^{\prime}} ;\, \left\langle c\right\rangle ^{*}=\frac{\left\langle c\right\rangle }{c_{o}} .\end{equation}

Note that we choose the width of the solute cloud as the characteristic axial scale. This choice is different than that of Taylor (Reference Taylor1953) and Stone & Brenner (Reference Stone and Brenner1999), who chose a macroscopic, dimensional length $L$ along the channel to a detector. This is an important distinction as combined (and opposed) PDF and EOF profiles with zero area-averaged net (bulk) velocity result in a dynamic condition wherein the solute zone disperses but its centroid does not move. In the latter case, a distance to a detector L has no meaning. As in Taylor’s (Reference Taylor1953) original analysis, we assume that the magnitude of the area-averaged concentration field is much larger than the magnitude of the deviation therefrom. Next, we assume that the characteristic axial width of the solute cloud is much greater than the inner radius of the tube ( $\sigma _{x}\gg a$ ). Further assuming the characteristic advection time of the cloud over $\sigma _{x}$ is large compared with radial diffusion time, we consider three smallness parameters as follows:

(3.5) \begin{equation} \varepsilon =\frac{a}{\sigma _{x}}\lt \lt 1 ;\, \frac{a}{\sigma _{x}}\frac{\left\langle u\right\rangle a}{D}=\varepsilon \textit{Pe}\lt \lt 1 ;\, \frac{c_{o}{}^{\prime}}{c_{o}}\lt \lt 1 .\end{equation}

Here, $\textit{Pe}=a\langle u\rangle /D$ is the flow Péclet number. Our non-dimensionalisation of (3.3) then becomes

(3.6) \begin{align} &\varepsilon \textit{Pe}\frac{c_{o}{}^{\prime}}{c_{o}}\frac{\partial c^{\prime*}}{\partial t^{*}}+\varepsilon \textit{Pe} \,{\hat u^{\prime}}\frac{\partial \!\left\langle c\right\rangle^{*}}{\partial x_{\sigma ^{\ast}}}+\varepsilon \textit{Pe}\frac{c_{o}{}^{\prime}}{c_{o}}{\langle \hat u\rangle }\frac{\partial c^{\prime*}}{\partial x_{\sigma^{\ast}}}+\varepsilon \textit{Pe}\frac{c_{o}{}^{\prime}}{c_{o}}{\hat u^{\prime}}\frac{\partial c^{\prime*}}{\partial x_{\sigma^{\ast}}}=\nonumber\\ &\frac{c_{o}{}^{\prime}}{c_{o}}\!\left(\frac{1}{r^{*}}\frac{\partial }{\partial r^{*}}\!\left(r^{*}\frac{\partial c^{\prime*}}{\partial r^{*}}\right)\right)+\varepsilon ^{2}\frac{c_{o}{}^{\prime}}{c_{o}}\frac{\partial ^{2}c^{\prime*}}{\partial x_{\sigma^{\ast}}^{2}}+\varepsilon \textit{Pe}\frac{c_{o}{}^{\prime}}{c_{o}}\!\left\langle {\hat u^{\prime}}\frac{\partial c^{\prime*}}{\partial x_{\sigma ^{\ast}}}\right\rangle . \end{align}

Keeping only terms with one smallness parameter, and keeping only the radial term in dimensionless form, we find

(3.7) \begin{equation} u^{\prime}\!\left(r\right)\frac{\partial \!\left\langle c\right\rangle }{\partial x}\approx \frac{D}{a^{2}}\!\left(\frac{1}{r^{*}}\frac{\partial }{\partial r^{*}}\!\left(r^{*}\frac{\partial c^{\prime}}{\partial r^{*}}\right)\right) .\end{equation}

This approximate balance between the non-uniform advective dispersion and the radial diffusion is a classic result for Taylor dispersion (Taylor Reference Taylor1953). Next, we multiply (3.7) by $r^{*}$ and integrate,

(3.8) \begin{align} &\frac{\partial \!\left\langle c\right\rangle }{\partial x}\int \left[r^{*}\!\left(\langle u_{p}\rangle \!\left(1-2r^{*2}\right)+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta -\frac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right)\right)\right]\mathrm{d}r^{*}= \nonumber \\&\quad\frac{\partial \!\left\langle c\right\rangle }{\partial x}\!\left(\langle u_{p} \rangle \!\left(\frac{r^{*2}}{2}-\frac{r^{*4}}{2}\right) +\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta \frac{r^{*2}}{2}-r^{*}\frac{I_{1}\!\left(\phi r^{*}\right)}{\phi I_{0}\!\left(\phi \right)}\right)\right)\approx \frac{D}{a^{2}}\!\left(r^{*}\frac{\partial c^{\prime}}{\partial r^{*}}\right)+C. \end{align}

We use the boundary condition at $r^{*}=1$ to see that the constant of integration $C$ must vanish. We divide by $r^{*}$ , multiply both sides by $a^{2}/D$ , and integrate from 0 to $r^{*}$ , arriving at

(3.9) \begin{align} c^{\prime}\!\left(r,x,t\right)&=c^{\prime}\!\left(0,x,t\right)+\frac{a^{2}}{D}\frac{\partial \!\left\langle c\right\rangle }{\partial x}\!\left(\frac{ \langle u_{p}\rangle }{4}\!\left(r^{*2}-\frac{r^{*4}}{2}\right)\right.\nonumber \\&\quad+\left.\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta \frac{r^{*2}}{4}-\frac{I_{0}\!\left(\phi r^{*}\right)-1}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right). \end{align}

Applying cross-sectional averaging, we are left with

(3.10) \begin{equation} 0=c^{\prime}\!\left(0,x,t\right)+\frac{a^{2}}{D}\frac{\partial \!\left\langle c\right\rangle }{\partial x}\!\left(\frac{ \langle u_{p} \rangle }{12}+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\frac{\eta }{8}-\frac{\eta }{\phi ^{2}}+\frac{1}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right) .\end{equation}

Subtracting (3.10) from (3.9), we derive

(3.11) \begin{align} c^{\prime}\!\left(r,x,t\right)&=\frac{a^{2}}{D}\frac{\partial \!\left\langle c\right\rangle }{\partial x}\!\left(\frac{ \langle u_{p} \rangle }{4}\!\left(r^{*2}-\frac{1}{3}-\frac{r^{*4}}{2}\right)\right.\nonumber \\&\quad+\left.\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(-\frac{\eta }{8}+\frac{\eta }{\phi ^{2}}+\eta \frac{r^{*2}}{4}-\frac{I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right) . \end{align}

Equation (3.11), combined with our solution for $\partial \langle c\rangle /\partial x$ , will provide an analytical expression for the unsteady, two-dimensional (radial and axial) distribution of the solute. To our knowledge, this is the first reported expression for the radial concentration of solute for combined EOF and PDF. We will use (3.11) to show solute concentration profiles versus axial position for several values of $r^{*}$ in figure 3. We will then visualise the two-dimensional solute concentration field in figures 4 and 5, and visualise the radial solute distribution alone in figure 6. Returning to (3.9), we take the derivative with respect to $x$ , multiply by $u^{\prime}(r)$ , and apply cross-sectional averaging

Figure 3.

The four figure panels show axial solute concentration profiles corresponding to four combinations of the Péclet number based on EOF, $\textit{Pe}_{e}$ , the Péclet number based on PDF, $\textit{Pe}_{p}$ , and the tube radius scaled by Debye length, $\phi$ . The ordinate shows dimensionless solute concentration, $c^{*}(x^{\prime*}|r^{*},\tau )=\pi a^{3}c(x^{\prime*}|r^{*},\tau )/N$ , and the abscissa is dimensionless axial position in a frame moving at the solute bulk velocity, $x^{\prime*}=x^{*}-\tau P\hspace{0pt}e$ . In each panel, solute concentration curves are shown at three values of dimensionless time (different line colours) for three distinct dimensionless radial positions (different line types), resulting in nine total concentration profiles per panel. The dimensionless solute concentration is obtained by evaluating (3.40) for an initial ‘top hat’ of solute with a width negligible to the final axial width of the solute and then normalising the solute concentration as $c^{*}(x^{\prime*}|r^{*},\tau )=\pi a^{3}c(x^{\prime*}|r^{*},\tau )/N$ .

Figure 4.

Benchmark comparisons between Brownian dynamics simulations and the analytical solution of the quasi-steady solute concentration field at three values of dimensionless time, $\tau =\textit{tD}/a^{2}$ . The figure shows dimensionless radius, $r^{*}$ , on the ordinate and dimensionless axial position, $x^{*}$ , on the abscissa. The top half of each panel shows individual particles from the Brownian dynamics simulations and the bottom half of each panel shows the concentration field predicted by the analytical solution in (3.40). Panels show four combinations of the Péclet number based on PDF, $P\hspace{0pt}e_{p}$ , the Péclet number based on EOF, $P\hspace{0pt}e_{e}$ , and the tube radius scaled by Debye length, $\phi$ . The far left of the panels shows the flow profile used to generate the data for each row. Figure was produced with a Péclet number of $\textit{Pe}=100$ (§ B of the SM has similar figures for $\textit{Pe}=$ 20 and 1000).

Figure 5.

Benchmark comparisons between Brownian dynamics simulations and the analytical solution of the quasi-steady solute concentration field at three values of dimensionless time, $\tau =\textit{tD}/a^{2}$ . Solute zones were subject to flows wherein EOF is perfectly opposed by PDF to achieve a net flow with zero area-averaged bulk velocity. Panels show solute zones for four values of the tube radius scaled by Debye length, $\phi$ . The figure shows dimensionless radius, $r^{*}$ , on the ordinate and dimensionless axial position, $x^{*}$ , on the abscissa. The top half of each panel shows individual particles from the Brownian dynamics simulations and the bottom half of each panel shows the concentration field predicted by the analytical solution in (3.40). The far left of the panels shows the flow profile used to generate the data for each row. Figure was produced with Péclet numbers based on EOF bulk velocity and PDF bulk velocity of $\textit{Pe}_{e}=100$ and $\textit{Pe}_{p}=-100$ , respectively.

Figure 6.

Comparisons between the analytical solution and Brownian dynamics simulations of the normalised quasi-steady radial distribution of solute at three values of dimensionless time, $\tau$ . The radial distribution is normalised as $c^{\prime}(r^{*},x^{*},\tau )/\langle c\rangle (x^{*}=\tau \textit{Pe},\tau )$ where $\langle c\rangle (x^{*}=\tau \textit{Pe},\tau )$ is the area-averaged concentration at the mean axial position of solute. Dimensionless radius, $r^{*}$ , is shown on the ordinate and dimensionless axial position, $x^{*}$ , is shown on the abscissa. The top half of each panel shows the radial solute distribution calculated from the Brownian dynamics simulations and the bottom half of each panel shows the analytical solution (given by the second term on the right-hand side of (3.40)). Panels show four combinations of the fraction of bulk velocity caused by pressure, $\beta$ , and the tube radius scaled by Debye length, $\phi$ . The dotted black lines denote two axial standard deviations of the Brownian particles’ position about their mean position. Figure was produced with a Péclet number of $\textit{Pe}=100$ .

(3.12) \begin{equation} \left\langle u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}\right\rangle =\frac{2a^{2}}{D}\frac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\!\int \limits_{0}^{1}\!\left[\!\!\begin{array}{l} r^{*}\!\left( \langle u_{p} \rangle \!\left(1-2r^{*2}\right)+\dfrac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta -\dfrac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right)\right) \times \\[6pt]\left(\dfrac{ \langle u_{p} \rangle }{4}\!\left(r^{*2}-\dfrac{r^{*4}}{2}\right)+\dfrac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\eta \dfrac{r^{*2}}{4}-\dfrac{I_{0}\!\left(\phi r^{*}\right)-1}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right) \end{array}\!\!\!\right]\!\mathrm{d}r^{*}. \end{equation}

Expanding the right-hand side of (3.12) and simplifying yields 17 terms. We collect these into three groups as follows:

(3.13) \begin{align} &\left\langle u^{\prime}\!\left(r\right)\dfrac{\partial c^{\prime}}{\partial x}\right\rangle =\dfrac{2a^{2}}{D}\dfrac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\nonumber\\&\quad \times \int \limits_{0}^{1} \left[\begin{array}{l} \dfrac{ \langle u_{p} \rangle ^{2}}{4}\!\left(r^{*3}-\dfrac{5r^{*5}}{2}+r^{*7}\right)\\[9pt] +\dfrac{ \langle u_{p} \rangle \!\left\langle u_{e}\right\rangle }{1-\eta }\!\left(\begin{array}{l} \dfrac{\eta r^{*3}}{2}-\dfrac{5\eta r^{*5}}{8}-\dfrac{r^{*}I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}+\dfrac{2r^{*3}I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}\\[9pt] -\dfrac{r^{*3}I_{0}\!\left(\phi r^{*}\right)}{4I_{0}\!\left(\phi \right)}+\dfrac{r^{*5}I_{0}\!\left(\phi r^{*}\right)}{8I_{0}\!\left(\phi \right)}+\dfrac{r^{*}}{\phi ^{2}I_{0}\!\left(\phi \right)}-\dfrac{2r^{*3}}{\phi ^{2}I_{0}\!\left(\phi \right)} \end{array}\right)\\[26pt] +\left(\dfrac{\left\langle u_{e}\right\rangle }{1-\eta }\right)^{2}\!\left(\begin{array}{l} \dfrac{\eta ^{2}r^{*3}}{4}-\dfrac{\eta r^{*}I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}-\dfrac{\eta r^{*3}I_{0}\!\left(\phi r^{*}\right)}{4I_{0}\!\left(\phi \right)}\\[9pt] +\dfrac{r^{*}I_{0}\!\left(\phi r^{*}\right)^{2}}{\phi ^{2}I_{0}\!\left(\phi \right)^{2}}+\dfrac{\eta r^{*}}{\phi ^{2}I_{0}\!\left(\phi \right)}-\dfrac{r^{*}I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)^{2}} \end{array}\right) \end{array}\right]\mathrm{d}r^{*} .\nonumber\\[5pt] \end{align}

Evaluating the integral in (3.13), we obtain

(3.14) \begin{equation} \left\langle u^{\prime}\!\left(r\right)\dfrac{\partial c^{\prime}}{\partial x}\right\rangle =-\dfrac{a^{2}}{D}\dfrac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\!\left(\begin{array}{l} \dfrac{ \langle u_{p} \rangle ^{2}}{48}+ \langle u_{p} \rangle \!\left\langle u_{e}\right\rangle \dfrac{\eta }{1-\eta }\!\left(\dfrac{1}{12}-\dfrac{2}{\phi ^{2}}+\dfrac{16}{\phi ^{4}}\!\left(\dfrac{1-\eta }{\eta }\right)\right)\\ +\left\langle u_{e}\right\rangle ^{2}\!\left(\dfrac{\eta }{1-\eta }\right)^{2}\!\left(\dfrac{3}{8}+\dfrac{2}{\phi ^{2}}-\dfrac{1}{\eta \phi ^{2}}-\dfrac{1}{\eta ^{2}\phi ^{2}}\right) \end{array}\right) .\end{equation}

Table 1 lists the relevant integrals involving modified Bessel functions of the first kind which were used in this evaluation. Next, for a compact notation, we define

(3.15) \begin{align} \gamma _{p}&=\frac{1}{48}\nonumber\\ \gamma _{pe}&=\frac{\eta }{1-\eta }\!\left(\frac{1}{12}-\frac{2}{\phi ^{2}}+\frac{16}{\phi ^{4}}\!\left(\frac{1-\eta }{\eta }\right)\right)\nonumber\\ \gamma _{e}&=\left(\frac{\eta }{1-\eta }\right)^{2}\!\left(\frac{3}{8}+\frac{2}{\phi ^{2}}-\frac{1}{\eta \phi ^{2}}-\frac{1}{\eta ^{2}\phi ^{2}}\right), \end{align}

which allows us to express (3.14) as

(3.16) \begin{equation} \left\langle u^{\prime}\!\left(r\right)\frac{\partial c^{\prime}}{\partial x}\right\rangle =-\frac{a^{2}}{D}\frac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\!\left( \langle u_{p} \rangle ^{2}\gamma _{p}+ \langle u_{p}\rangle \!\left\langle u_{e}\right\rangle \gamma _{pe}+\left\langle u_{e}\right\rangle ^{2}\gamma _{e}\right) .\end{equation}

Substituting (3.16) into (3.2), we derive the following one-dimensional, unsteady partial differential equation describing the area-averaged concentration field:

(3.17) \begin{equation} \frac{\partial \!\left\langle c\right\rangle }{\partial t}+\left\langle u\right\rangle \frac{\partial \!\left\langle c\right\rangle }{\partial x}=\frac{\partial ^{2}\!\left\langle c\right\rangle }{\partial x^{2}}\!\left(D+\frac{a^{2}}{D}\!\left( \langle u_{p}\rangle ^{2}\gamma _{p}+ \langle u_{p} \rangle \!\left\langle u_{e}\right\rangle \gamma _{pe}+\left\langle u_{e}\right\rangle ^{2}\gamma _{e}\right)\right) .\end{equation}

Table 1.

Integrals applicable to the evaluation of (3.13).

The first, third and second parenthetic terms on the right-hand side respectively capture the effects of dispersion caused by PDF alone, by EOF alone and by the coupled, non-superposable effects of simultaneous PDF and EOF. From the latter equation, the effective dispersion coefficient valid for the Taylor-type limit of the problem is

(3.18) \begin{equation} D_{\textit{eff}}=D+\frac{a^{2}}{D}\!\left( \langle u_{p} \rangle ^{2}\gamma _{p} + \langle u_{p} \rangle \!\left\langle u_{e}\right\rangle \gamma _{pe}+\left\langle u_{e}\right\rangle ^{2}\gamma _{e}\right) .\end{equation}

Note that, in analogy to Taylor’s (Reference Taylor1953) original analysis, (3.17) has the form of a one-dimensional unsteady diffusion equation with the molecular diffusivity replaced by an effective dispersion coefficient given by (3.18). As such, the solution is an initial value problem that requires an initial axial distribution of the area-averaged solute. However, the solution of (3.17) does not accurately predict the evolution of solute until sufficient time has passed relative to the characteristic time of transverse diffusion, $a^{2}/D$ . In other words, this solution is valid as long as the solute variance gained in early times (during the pre-asymptotic dispersion regime) is small compared with either the initial variance or the variance gained due to the (long-term) Taylor dispersion. We will soon explore predictions by this solution given specific initial conditions. In § 4, we will derive dispersion solutions (given various initial conditions) which are valid for both short and long times. Next, we note that (3.18) is equivalent to the solution reported by Datta & Kotamarthi (Reference Datta and Kotamarthi1990) for the dimensional dispersion coefficient (see their equation (37)). Datta & Kotamarthi reported deriving their solution by applying a Taylor-type analysis but, as mentioned earlier, did not present a derivation in any detail. Further, we note that (3.18), for pure PDF ( $\langle u_{p}\rangle =\langle u\rangle$ ), reduces to Aris’ (Reference Aris1956) classic solution of $D_{\textit{eff}}=D+D^{-1}a^{2}\langle u_{p}\rangle ^{2}\gamma _{p}$ . Additionally, for pure EOF ( $\langle u_{e}\rangle =\langle u\rangle$ ), (3.18) becomes $D_{\textit{eff}}=D+D^{-1}a^{2}\langle u_{e}\rangle ^{2}\gamma _{e}$ , which agrees analytically with Griffiths & Nilson’s (Reference Griffiths and Nilson1999) solution of the coefficient of effective dispersion (their equation (38); see § A of the supplementary material (SM) for proof).

Next, we scale the expression for the dispersion coefficient. To this end, we note that $\langle u_{p}\rangle$ and $\langle u_{e}\rangle$ are signed quantities. They are also variables which are fully independent of each other (e.g. they can be independently controlled in an experiment). Hence, we define here two corresponding (signed) Péclet numbers of the form $P\hspace{0pt}e_{p}=a\langle u_{p}\rangle /D$ and $P\hspace{0pt}e_{e}=a\langle u_{e}\rangle /D$ . For both, we use the characteristic length scale of the cylinder radius. The non-dimensional dispersion coefficient is then

(3.19) \begin{equation} D_{\textit{eff}}^{*}=\frac{D_{\textit{eff}}}{D}=1+\textit{Pe}_{p}^{2}\gamma _{p}+\textit{Pe}_{p}\textit{Pe}_{e} \gamma _{pe}+\textit{Pe}_{e}^{2}\gamma _{e} .\end{equation}

The last three terms on the right-hand side again capture the effects of pure PDF, coupled PDF and EOF and pure EOF. Because the values of $\gamma$ are functions of $\phi$ (cf. (3.15)), $D_{\textit{eff}}^{*}$ is a function of only three dimensionless parameters: $\textit{Pe}_{p}, \textit{Pe}_{e}$ and $\phi$ .

Before continuing our presentation, we briefly compare our (3.19) with the non-dimensional dispersion coefficient presented by Datta & Kotamarthi (Reference Datta and Kotamarthi1990). The latter authors chose to combine the signed quantities $\langle u_{p}\rangle$ and $\langle u_{e}\rangle$ as a sum within a single global Péclet number of the form $\textit{Pe}=(\langle u_{p}\rangle +\langle u_{e}\rangle )a/D$ . Further, they chose to use the ratio $\beta =\langle u_{p}\rangle /\langle u\rangle$ to distinguish between the effects of EOF and PDF. Their approach yields the following result:

(3.20) \begin{equation} D_{\textit{eff}}^{*}=1+P\hspace{0pt}e^{2}\!\left(\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right) .\end{equation}

In this formulation, $D_{\textit{eff}}^{*}$ is again a function of three parameters, but these are now $\phi, \beta$ and $P\hspace{0pt}e$ . We note this combination of two velocity scales into a single Péclet number is most useful to characterise dispersion when EOF and PDF act in the same direction. We recommend the formulation of (3.19) for cases wherein EOF and PDF are in opposition. When EOF and PDF are opposed and nearly equal in magnitude, $P\hspace{0pt}e$ approaches zero while $\beta$ approaches infinity at the same rate and evaluation of (3.20) is more cumbersome than (3.19).

We can rewrite (3.20) as follows:

(3.21) \begin{equation} \frac{48}{\textit{Pe}^{2}}\!\left(D_{\textit{eff}}^{*}-1\right)=48\left(\beta ^{2}\gamma _{p}+\beta \!\left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right) .\end{equation}

This formulation enables exploration of the full solution and in a manner most useful when EOF and PDF are not perfectly opposed. We will use this formulation to plot and benchmark the grouping $48\textit{Pe}^{-2}(D_{\textit{eff}}^{*}-1)$ in heatmaps where the independent variables are $\beta$ and $\phi$ (see figure 7). Further, we will benchmark $D_{\textit{eff}}^{*}$ against Brownian dynamics simulations for the case of perfectly opposed PDF and EOF in figure 8. Lastly, we note (3.21) conveniently reduces to unity for the case of parabolic flow. Parabolic flow is achieved by pure PDF or for combinations of PDF and EOF with thick EDLs (since thick-EDL solutions also lead to parabolic EOF).

Figure 7.

Contour plots of our quasi-steady analytical solution $(a)$ versus Brownian dynamics simulations $(b)$ of the normalised effective dispersion coefficient. In both cases, the effective dispersion coefficient is normalised as $48\textit{Pe}^{-2}(D_{\textit{eff}}^{*}-1)$ where $D_{\textit{eff}}^{*}$ is the effective dispersion coefficient. This non-dimensional quantity is plotted versus the fraction of bulk velocity caused by pressure, $\beta$ , and the radius scaled by Debye length, $\phi$ . The panels cover the entire dynamics of quasi-steady dispersion for all relative EDL thicknesses and many velocity profiles. The white dashed lines show the variance-minimising contour for values of $\beta$ as a function of $\phi$ . These lines are obtained from (3.32) for the analytical solution and numerically for the Brownian dynamics simulations.

Figure 8.

Contour plots of our quasi-steady analytical solution $(a)$ versus Brownian dynamics simulations $(b)$ of the dimensionless effective dispersion coefficient for the case of perfectly opposed EOF and PDF. The dimensionless dispersion coefficient $D_{\textit{eff}}^{*}$ is plotted versus the Péclet number based on PDF, $P\hspace{0pt}e_{p}$ , and the radius scaled by Debye length, $\phi$ . The Péclet number based on PDF is fixed to be equal in magnitude but opposite in sign to the Péclet number based on EOF, $P\hspace{0pt}e_{e}$ . For perfectly opposed PDF and EOF, $D_{\textit{eff}}^{*}$ decreases as $\phi$ decreases.

We now continue with our derivation for the two-dimensional concentration field. We first evaluate the solution of (3.17) for an initial delta function condition for the solute of the form

(3.22) \begin{equation} \left\langle c\right\rangle \!\left(0,0\right)=\frac{N}{A}\delta (0), \end{equation}

where $N/A$ is the moles of solute per unit cross-sectional area of the channel, and the delta function $\delta (x)$ has dimensions of inverse length. Defining $x^{\prime}=x-\langle u\rangle t$ as an axial coordinate for a reference frame moving at net bulk velocity, the area-averaged solute distribution for this fundamental case is

(3.23) \begin{equation} \langle c\rangle (x^{\prime},t)=\frac{N}{A}\frac{1}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right) .\end{equation}

This solution can be interpreted as the Green’s function (also know as heat kernel) for (3.17) which, as we discuss in the next section, can be used to construct a wide variety of solutions for more complex initial conditions. Further, the first derivative of (3.23) with respect to axial dimension is

(3.24) \begin{equation} \frac{\partial \!\left\langle c\right\rangle }{\partial x^{\prime}}=-\frac{N}{A}\frac{x^{\prime}}{4\sqrt{\pi \!\left(D_{\textit{eff}}\,t\right)^{3}}}\exp\!\left(\frac{-(x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)=-\frac{x^{\prime}}{2D_{\textit{eff}}\,t}\!\left\langle c\right\rangle \! (x^{\prime},t) .\end{equation}

Equations (3.23) and (3.24) can be combined with (3.11) to obtain a solution of the full two-dimensional, unsteady solute concentration for a delta function initial condition as follows:

(3.25) \begin{align} c\!\left(r^{*},x^{\prime},t\right)&=\left\langle c\right\rangle +c^{\prime}=\frac{N}{A}\frac{1}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)\nonumber\\&\quad - \frac{Na^{2}x^{\prime}}{4AD\sqrt{\pi \!\left(D_{\textit{eff}}\,t\right)^{3}}}\exp\!\left(\frac{-(x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)\left(\frac{ \langle u_{p}\rangle }{4}\!\left(r^{*2}-\frac{1}{3}-\frac{r^{*4}}{2}\right)\right.\nonumber\\&\quad +\left.\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(-\frac{\eta }{8}+\frac{\eta }{\phi ^{2}}+\eta \frac{r^{*2}}{4}-\frac{I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right). \end{align}

To our knowledge, (3.25) is the first reported expression which provides both the axial and radial distribution of solute for combined EOF and PDF.

3.1. Unsteady, two-dimensional solutions for arbitrary initial conditions

Given the linearity of (3.17), we leverage a Green’s function to present the concentration field as a function of arbitrary initial conditions. The Green’s function for this case is the solution of (3.17) for a delta function initial condition (Cole Reference Cole2011). From (3.23), the Green’s function is

(3.26) \begin{equation} G\!\left(x^{\prime},t|x_{i}\right)=\frac{1}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-\left(x^{\prime}-x_{i}\right)^{2}}{4D_{\textit{eff}}\,t}\right) .\end{equation}

Whence, combining with (3.11), we construct the solution to some arbitrary initial condition $\langle c\rangle (x^{\prime},0)$ as follows:

(3.27) \begin{align} &c\!\left(r^{*},x^{\prime},t\right)\nonumber \\&\quad =\int\limits_{-\infty }^{\infty }G\!\left(x^{\prime},t|x_{i}\right)\!\left\langle c\right\rangle \!\left(x_{i},0\right)\!\mathrm{d}x_{i} +f\!\left(\frac{\partial }{\partial x^{\prime}}\int\limits_{-\infty }^{\infty }G\!\left(x^{\prime},t|x_{i}\right)\!\left\langle c\right\rangle \!\left(x_{i},0\right)\! \mathrm{d}x_{i},r^{*}\right)\nonumber\\&\quad = \int\limits_{-\infty }^{\infty }G\!\left(x^{\prime},t|x_{i}\right)\!\left\langle c\right\rangle \!\left(x_{i},0\right)\! \mathrm{d}x_{i} -\frac{a^{2}}{2D_{\textit{eff}}Dt}\!\left(\int\limits_{\,-\infty }^{\infty }\!\left(x^{\prime}-x_{i}\right)\!G\!\left(x^{\prime},t|x_{i}\right)\! \left\langle c\right\rangle \!\left(x_{i},0\right)\! \mathrm{d}x_{i}\!\right)\nonumber\\&\qquad \times \left(\frac{ \langle u_{p} \rangle }{4}\!\left(r^{*2}-\frac{1}{3}-\frac{r^{*4}}{2}\right)+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(-\frac{\eta }{8}+\frac{\eta }{\phi ^{2}}+\eta \frac{r^{*2}}{4}-\frac{I_{0}\!\left(\phi r^{*}\right)\!}{\phi ^{2}I_{0}\!\left(\phi \right)\!}\right)\!\right)\!. \end{align}

3.2. Variance growth valid for any initial condition

We derive a general expression for the growth of the solute area-averaged variance valid for any initial condition. Let $v_{2}$ denote the axial variance of the solute about the axial mean defined such that

(3.28) \begin{equation} v_{2}(\langle c \rangle \kern2pt \!(x',t)) = \frac{1}{N}\int\limits _{-\infty }^{\infty }x^{\prime2}\!\left\langle c\right\rangle \! (x^{\prime},t)\mathrm{d}x^{\prime} .\end{equation}

Equation (3.28) can be expressed in terms of $G(x^{\prime},t|x_{i})$ as follows:

(3.29) \begin{align} v_{2}\!\left(\left\langle c\right\rangle \!(x^{\prime},t)\right)&=\frac{1}{N}\int\limits_{-\infty }^{\infty }x^{\prime2}\int\limits_{-\infty }^{\infty }G (x^{\prime},t|x_{i})\!\left\langle c\right\rangle \!\left(x_{i},0\right)\!\mathrm{d}x_{i}\mathrm{d}x^{\prime}\nonumber\\&= \frac{1}{N}\int\limits_{-\infty }^{\infty }\!\left\langle c\right\rangle \!\left(x_{i},0\right)\int\limits_{-\infty }^{\infty }x^{\prime2}G (x^{\prime},t|x_{i}) \mathrm{d}x^{\prime}\mathrm{d}x_{i}\nonumber\\&= \frac{1}{N}\int\limits_{-\infty }^{\infty }\!\left\langle c\right\rangle \!\left(x_{i},0\right)\int\limits_{-\infty }^{\infty }x^{\prime2}\!\left(\frac{1}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime}-x_{i})^{2}}{4D_{\textit{eff}}\,t}\right)\right)\!\mathrm{d}x^{\prime} \mathrm{d}x_{i}. \end{align}

We recognise the last interior integral as the second moment of a Gaussian with mean $x_{i}$ and variance $2D_{\textit{eff}}\,t$ . This implies

(3.30) \begin{align} &\frac{1}{N}\int\limits_{-\infty }^{\infty }\!\left\langle c\right\rangle \!\left(x_{i},0\right)\int\limits_{-\infty }^{\infty }x^{\prime2}\!\left(\frac{1}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime}-x_{i})^{2}}{4D_{\textit{eff}}\,t}\right)\right)\!\mathrm{d}x^{\prime}\mathrm{d}x_{i}=\nonumber\\ & \frac{2D_{\textit{eff}}\,t}{N}\int\limits_{-\infty }^{\infty }\!\left\langle c\right\rangle \!\left(x_{i},0\right)\!\mathrm{d}x_{i}+\frac{1}{N}\int\limits_{-\infty }^{\infty }{x_{i}}^{2}\!\left\langle c\right\rangle \!\left(x_{i},0\right)\!\mathrm{d}x_{i}. \end{align}

Therefore, variance grows as

(3.31) \begin{equation} v_{2}(\left\langle c\right\rangle \! (x^{\prime},t))=2D_{\textit{eff}}\,t+v_{2}(\langle c\rangle \kern2pt \!(x^{\prime},0)) .\end{equation}

That is, the solute variance for all time is simply the sum of its initial value and $2D_{\textit{eff}}\,t$ .

3.3. Minimisation of variance for transporting solute over arbitrary axial distances

We here derive analytical expressions and present methods to minimise the axial variance of solute as it is transported over an arbitrary, non-zero axial distance. As such, we do not consider the case of $\langle u\rangle =0$ (e.g. wherein EOF and PDF are perfectly opposed), ensuring the parameter $\beta =\langle u_{p}\rangle /\langle u\rangle$ remains finite.

First, we consider the relatively rare case where the practitioner can vary the relative EDL thickness ( $\phi$ ) as is possible in some nanochannel systems. To minimise variance growth, the EDL should be made as relatively thin as possible ( $\phi \rightarrow \infty$ ) as this minimises dispersion.

Figure 9.

The fraction of bulk velocity associated with pressure required to minimise the growth of the solute variance, $\beta _{o}$ , as a function of the tube radius scaled by Debye length, $\phi$ . The inset shows the variance-minimising normalised flow profiles, $\hat{u}_{o}(r^{*})$ , for select values of $\phi$ of 1, 5, 10, 25 and 50. Opposing EOF with PDF can be used to suppress dispersion.

Next, we consider a more practical situation. Assume the channel geometry and Debye length of a system are fixed. This fixes the relative EDL thickness $\phi$ . To minimise variance, (3.20) shows that we should consider both $\beta$ and $\textit{Pe}$ . From (3.31), we see that minimisation of variance in terms of the parameter $\beta$ implies a minimisation of $D_{\textit{eff}}$ . We proceed by taking the partial derivative of (3.20) with respect to $\beta$ and setting the resulting expression equal to zero. The variance-minimising value of $\beta$ , denoted here as $\beta _{o}$ , is then

(3.32) \begin{equation} \beta _{o}=\frac{\gamma _{pe}-2\gamma _{e}}{2\left(\gamma _{pe}-\gamma _{e}-\gamma _{p}\right)} .\end{equation}

Interestingly, and somewhat counterintuitively, we find that $\beta _{o}$ is determined completely by $\phi$ and is independent of $P\hspace{0pt}e$ . For any one value $\phi$ , there is a unique combination of EOF and PDF which minimises variance and is applicable to all bulk velocities, tube radii and solute diffusivities (which determine $P\hspace{0pt}e$ ). We shall see that this optimum velocity profile shape is achieved by opposing EOF with progressively stronger PDF as the EDL becomes thicker (we will discuss this further when presenting figure 9). Datta & Kotamarthi (Reference Datta and Kotamarthi1990) also identified this optimum $\beta _{o}$ (see equation (50) of Datta & Kotamarthi Reference Datta and Kotamarthi1990). However, we note that Datta & Kotamarthi focused their presentation on the minimisation of the theoretical plate height. Theoretical plate heights are traditional figures of merit for species separation systems, including electrophoresis and chromatography, which aim to capture the number of peaks detectable using a certain column (i.e. channel) length (see Huang Reference Huang2021). We here will instead focus on the more fundamental problem of minimising $D_{\textit{eff}}$ and axial variance. Lastly, we henceforth denote the flow profile based on a value of $\beta =\beta _{o}$ as $u_{o}(r^{*})$ (and define $\hat{u}_{o}(r^{*})=u_{o}(r^{*})/\langle u\rangle$ ).

Next, we explore the effects of flow Péclet number on the minimum dispersion conditions. From (3.31), assuming some negligible initial width of the solute (compared with its final value), the final variance of the solute cloud reduces simply to $2D_{\textit{eff}}\,t$ . First, we non-dimensionalise the solute variance as $2D_{\textit{eff}}^{*}\tau =2D_{\textit{eff}}\,t/a^{2}$ , where $\tau =tD/a^{2}$ is the dimensionless time. We wish to minimise $2D_{\textit{eff}}^{*}\tau$ for some finite distance travelled $L$ or, in dimensionless form, $L^{*}=L/a$ . To proceed, we first multiply (3.20) by the dimensionless time $\tau =L^{*}/\textit{Pe}$

(3.33) \begin{equation} \frac{L^{*}}{\textit{Pe}}D_{\textit{eff}}^{*}=\frac{L^{*}}{\textit{Pe}}+L^{*}P\hspace{0pt}e\left(\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right) .\end{equation}

We differentiate the right-hand side of (3.33) with respect to $P\hspace{0pt}e$ and set this equal to zero

(3.34) \begin{equation} -\frac{L^{*}}{\textit{Pe}^{2}}+L^{*}\!\left(\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right)=0 .\end{equation}

Recall we here assume $\textit{Pe}\neq 0$ . Multiplying (3.34) by $P\hspace{0pt}e^{2}$ , the roots of the resulting equation can be expressed as

(3.35) \begin{equation} P\hspace{0pt}e_{\mathrm{o}}=\pm \sqrt{\frac{1}{\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}}} .\end{equation}

This compact expression yields the Péclet number that minimises the rate of growth of the variance per unit of axial distance travelled. Since we assumed a negligible initial width of the solute, the resulting variance-minimising $\textit{Pe}$ is independent of $L^{*}$ . Note the $\pm$ simply reflects the symmetry of the system for bulk velocities in the ‘positive’ or ‘negative’ direction and can therefore be discarded. Equation (3.35) then provides the optimum Péclet number to minimise variance for transporting solute for any values of the parameters $\beta$ and $\phi$ . If one considers variations of $\beta$ (i.e. adjustments of the relative strength of EOF and PDF), then the substitution of (3.32) into (3.35) conveniently reduces $P\hspace{0pt}e_{\mathrm{o}}$ to a function of only the fixed parameter $\phi$ . If instead $\beta$ is taken to be fixed, then (3.35) provides $P\hspace{0pt}e_{\mathrm{o}}$ as a function of both $\phi$ and the value of $\beta$ . In either case, the relations $\beta \textit{Pe}=\textit{Pe}_{p}$ and $(1-\beta )\textit{Pe}=\textit{Pe}_{e}$ can be used to find the variance-minimising PDF and EOF components once $P\hspace{0pt}e_{\mathrm{o}}$ is known.

We note that (3.35) may be interpreted as the variance-minimising ratio of advective to diffusive flux for transporting solute. In regimes where advection only slightly increases dispersion (e.g. dominant EOF with a thin EDL), the dispersion is dominated by molecular diffusion. In such regimes, the time to transport solute should then be minimised by increasing the EOF velocity and $\textit{Pe}$ . Hence, $P\hspace{0pt}e_{\mathrm{o}}$ will be large for such velocity fields. We shall show (in the discussion of figure 10) that, for finite $\phi, P\hspace{0pt}e_{\mathrm{o}}$ achieves a maximum (with respect to $\beta$ ) at the variance-minimising combination of EOF and PDF, $\beta _{o}$ . To demonstrate this, we differentiate (3.35) with respect to $\beta$ and set the resulting quantity equal to zero

(3.36) \begin{align} \frac{\partial \textit{Pe}_{\mathrm{o}}}{\partial \beta }&=-\frac{1}{2}\!\left(\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right)^{-\frac{3}{2}}\nonumber\\&\quad\times\frac{\partial }{\partial \beta }\!\left(\beta ^{2}\gamma _{p}+\beta \left(1-\beta \right)\gamma _{pe}+\left(1-\beta \right)^{2}\gamma _{e}\right)=0 .\end{align}

Figure 10.

Surface plots for analytical solutions of the quasi-steady state, variance-minimising Péclet number, $\textit{Pe}_{o}$ . In $(a)$ , a surface of $\textit{Pe}_{o}$ is plotted as a function of the fraction of bulk velocity caused by pressure, $\beta$ , and the tube radius scaled by Debye length, $\phi$ . For constant values of $\phi, \textit{Pe}_{o}$ exhibits a maximum at the aforementioned variance-minimising fraction of bulk velocity caused by pressure, $\beta _{o}$ . In $(b)$ , a surface of $2D_{\textit{eff}}^{*}\tau$ is plotted for solute travelling an example dimensionless axial distance of $L^{*}=\tau \textit{Pe}=50$ as a function of $\beta$ and $P\hspace{0pt}e$ . Here, the white circle indicates the minimum of $2D_{\textit{eff}}^{*}\tau$ as determined analytically by (3.32) and (3.35). Figure 10(b) was produced with an example relative EDL thickness of $\phi =50$ .

The prefactor parenthetic term on the right-hand side (raised to $-3/2$ ) is always non-zero. Therefore, the right-hand side of (3.36) vanishes only when $\beta =\beta _{o}$ . Hence, we conclude that (at a fixed value of $\phi$ ) the $P\hspace{0pt}e_{\mathrm{o}}$ function achieves a maximum at $\beta _{o}$ .

Lastly, we note that our (3.35) is mathematically equivalent to Datta & Kotamarthi’s (Reference Datta and Kotamarthi1990) equation (51). However, we again point out that Datta & Kotamarthi derived their optimal $\textit{Pe}$ by minimising the theoretical plate height (and we here do not need to invoke such figures of merit).

3.4. Example solution for a ‘top hat’ initial condition

In the case of an initial ‘top hat’ of solute, centred at $x^{\prime}=0$ , the initial conditions are

(3.37) \begin{equation} \langle c\rangle (x^{\prime},0)=\left\{\begin{array}{ll} c_{i} & | x^{\prime}| \leq L_{o}\\[4pt] 0 & | x^{\prime}| \gt L_{o} \end{array}\right. .\end{equation}

Here, $c_{i}$ denotes the initial concentration of solute (with units of moles per litre) and $L_{o}$ denotes half of the initial width of the ‘top hat’. First, we compute the convolution of this initial condition with Green’s function as

(3.38) \begin{equation} \int\limits_{-L_{o}}^{L_{o}}\frac{c_{i}}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime}-x_{i})^{2}}{4D_{\textit{eff}}\,t}\right)\!\mathrm{d}x_{i}=\frac{c_{i}}{2}\!\left(\mathrm{erf}\!\left(\frac{L_{o}-x^{\prime}}{\sqrt{4D_{\textit{eff}}\,t}}\right)+\mathrm{erf}\!\left(\frac{L_{o}+x^{\prime}}{\sqrt{4D_{\textit{eff}}\,t}}\right)\right) .\end{equation}

For the radial distribution of solute, we calculate

(3.39) \begin{align} \int\limits_{-L_{o}}^{L_{o}}\frac{c_{i} (x^{\prime}-x_{i})}{\sqrt{4\pi D_{\textit{eff}}\,t}}\exp\!\left(\frac{-(x^{\prime}-x_{i})^{2}}{4D_{\textit{eff}}\,t}\right)\!\mathrm{d}x_{i}&=\frac{c_{i}\sqrt{D_{\textit{eff}}\,t}}{\sqrt{\pi }}\!\left(\exp\!\left(-\frac{(L_{o}-x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)\right.\nonumber\\&\quad-\left.\exp\!\left(-\frac{(L_{o}+x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)\right) .\end{align}

Substituting (3.38) and (3.39) into (3.27), we see that, in the ‘top hat’ case, the overall concentration field can be expressed as

(3.40) \begin{align} c(r^{*},x^{\prime},t)&= \frac{c_{i}}{2}\!\left(\mathrm{erf}\!\left(\frac{L_{o}-x^{\prime}}{\sqrt{4D_{\textit{eff}}\,t}}\right)+\mathrm{erf}\!\left(\frac{L_{o}+x^{\prime}}{\sqrt{4D_{\textit{eff}}\,t}}\right)\right)\nonumber\\ &\quad-\frac{a^{2}c_{i}}{2D\sqrt{\pi D_{\textit{eff}}\,t}}\!\left(\exp\!\left(-\frac{(L_{o}-x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)-\exp\!\left(-\frac{(L_{o}+x^{\prime})^{2}}{4D_{\textit{eff}}\,t}\right)\right)\nonumber\\ &\quad\times \left(\frac{\langle u_{p}\rangle }{4}\!\left(r^{*2}-\frac{1}{3}-\frac{r^{*4}}{2}\right)+\frac{\left\langle u_{e}\right\rangle }{1-\eta }\!\left(-\frac{\eta }{8}+\frac{\eta }{\phi ^{2}}+\eta \frac{r^{*2}}{4}-\frac{I_{0}\!\left(\phi r^{*}\right)}{\phi ^{2}I_{0}\!\left(\phi \right)}\right)\right). \end{align}

Again, (3.40) is valid for the Taylor-type limit of the problem (cf. (3.4) and (3.5)).

4. Method of moments solution for transient dispersive behaviour

We present a solution for the integral moments of the solute zone which is valid for all times. Our solution will also be valid for solute zones with arbitrary axial widths. We follow the MoM solution method originated by Aris (Reference Aris1956), as significantly expanded and developed by Barton (Reference Barton1983). Barton advanced the MoM by employing eigenfunction expansions from Sturm–Liouville theory to derive general expressions for the integral moments of the concentration field. Barton’s analysis applies to arbitrary flow profiles and arbitrary channel cross-sections. We will here use Barton’s general framework to analyse the behaviour of solute subject to coupled EOF and PDF in a cylindrical tube. We note there are several key differences between our MoM formulation and our quasi-steady, Taylor-type analysis from § 3. First, the conditions $\textit{Pe}\lt \lt \sigma _{x}/a$ and $a/\sigma _{x}\lt \lt 1$ are not required for our MoM formulation. Second, consistent with the MoM, we non-dimensionalise axial position with the inner radius of the tube and not the characteristic width of the solute cloud ( $x^{*}=x/a$ ). Third, we non-dimensionalise the concentration field with the total moles of solute in the tube and the characteristic volume of the tube ( $c^{*}(r^{*},x^{*},\tau )=\pi a^{3}c(r^{*},x^{*},\tau )/N$ ). Lastly, we non-dimensionalise our velocity profile as follows:

(4.1) \begin{equation} \chi (r^{*})=\frac{a}{D}u (r^{*})=2\textit{Pe}_{p}\big(1-r^{*2}\big)+\frac{\textit{Pe}_{e}}{1-\eta }\!\left(1-\frac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right) .\end{equation}

The latter can be interpreted as a Péclet number in terms of the (radially varying) velocity profile. In terms of these newly defined variables, the dimensionless convection–diffusion equation and boundary conditions are

(4.2) \begin{align} &\frac{\partial c^{*}}{\partial \tau }+\chi (r^{*})\frac{\partial c^{*}}{\partial x^{*}}=\frac{1}{r^{*}}\frac{\partial }{\partial r^{*}}\!\left(r^{*}\frac{\partial c^{*}}{\partial r^{*}}\right)+\frac{\partial ^{2}c^{*}}{\partial x^{*2}}; \left.\frac{\partial c^{*}}{\partial r^{*}}\right| _{r^{*}=1,0}=0,\nonumber\\ & (x^{*})^{n}\frac{\partial ^{m}c^{*}}{\partial x^{*m}}\rightarrow 0\text{ as }\!\left| x^{*}\right| \rightarrow \infty \text{ for }n, m\geq 0,c^{*}(r^{*},x^{*},\tau =0)=c^{*}(r^{*},x^{*},0). \end{align}

Defining $c _{n}^{*}(r^{*},\tau )=\int _{-\infty }^{\infty }(x^{*})^{n}c^{*}(r^{*},x^{*},\tau )\mathrm{d}x^{*}$ as the nth moment of solute along the axial direction, we can multiply (4.2) by $(x^*)^{n}$ and integrate with respect to $x^{*}$ , yielding

(4.3) \begin{align} &\frac{\partial c _{n}^{*}}{\partial \tau }=\frac{1}{r^{*}}\frac{\partial }{\partial r^{*}}\!\left(r^{*}\frac{\partial c _{n}^{*}}{\partial r^{*}}\right)+n\left(n-1\right)c _{n-2}^{*}+n \chi (r^{*})c _{n-1}^{*};\left.\frac{\partial c _{n}^{*}}{\partial r^{*}}\right| _{r^{*}=1,0}=0,\nonumber\\ & c _{n}^{*}(r^{*},\tau =0)=c _{n}^{*}(r^{*},0). \end{align}

Taking an area average of (4.3)

(4.4) \begin{equation} \frac{dM_{n}}{d\tau }=n\left(n-1\right)\langle c _{n-2}^{*} \rangle +n \langle \chi (r^{*})c _{n-1}^{*} \rangle ;\, M_{n}\!\left(\tau =0\right)=M_{n}\!\left(0\right) .\end{equation}

Here, $M_{n}(\tau )=\int _{0}^{1}2r^{*}c _{n}^{*}(r^{*},\tau )\mathrm{d}r^{*}$ and $M_{n}(0)$ are respectively the nth moment of the solute zone and the initial value of the nth moment of solute. Equation (4.4) is then the governing initial value problem for the nth moment of the concentration field. We note that solutions from (4.4) for the nth moment of solute depend upon the solutions of $c _{n-1}^{*}(r^{*},\tau )$ and $c _{n-2}^{*}(r^{*},\tau )$ from (4.3). In this way, solutions for the higher-order moments of solute must be ‘built up’ from solutions of lower-order moments of solute along the axial direction. To this end, we will solve (4.3) with the method of separation of variables. Given the cylindrical geometry of our formulation, the Sturm–Liouville eigenvalue problem (Amrein, Hinz & Pearson Reference Amrein, Hinz and Pearson2005) associated with our solution method is

(4.5) \begin{align} \left(\frac{d}{dr^{*}}\!\left(r^{*}\frac{d}{dr^{*}}\right)+\mu _{i}r^{*}\right)f_{i}=0;\left.\frac{df_{i}}{dr^{*}}\right| _{r^{*}=1}=0,f_{i}\,{\mathrm {is\,finite}}. \end{align}

The orthonormal eigenfunctions and discrete eigenvalues of (4.5) are

(4.6) \begin{equation} f_{i} (r^{*})=\frac{J_{0}\!\left(\alpha _{i}r^{*}\right)}{J_{0}\!\left(\alpha _{i}\right)};\mu _{i}=\alpha _{i}^{2} .\end{equation}

Here, $J_{v}(z)$ denotes a Bessel function of the first kind of order $v$ , and $\alpha _{i}$ is the ith root of $J_{1}(z)$ . As previously mentioned, Barton (Reference Barton1983) provided a general formulation of the first moment of the concentration field in terms of the eigenfunctions associated with arbitrary tube geometries and arbitrary velocity profiles. Importantly, Barton showed that the first moment of the concentration field satisfies the equation

(4.7) \begin{equation} M_{1}(\tau )=\tau \textit{Pe}+\sum _{i}^{\infty }\!\left\langle c _{0}^{*}(r^{*},0)f_{i}\right\rangle \,\left\langle \chi f_{i}\right\rangle \,\left(\frac{1-e^{-{\mu _{i}}\tau }}{\mu _{i}}\right) .\end{equation}

Without loss of generality, we have assumed that $M_{1}(0)=0$ . To apply Barton’s approach to the current geometry and velocity fields, we must evaluate the following eigenfunction-weighted area average of the velocity profile in (4.1)

(4.8) \begin{align} \left\langle \chi f_{i}\right\rangle &=2\int\limits_{0}^{1} \!\left[\frac{r^{*}J_{0} \big(\alpha _{i}r^{*}\big)}{J_{0}\!\left(\alpha _{i}\right)}\!\left(2\textit{Pe}_{p} \big(1-r^{*2}\big)+\frac{\textit{Pe}_{e}}{1-\eta }\!\left(1-\frac{I_{0}\!\left(\phi r^{*}\right)}{I_{0}\!\left(\phi \right)}\right)\right)\right]\mathrm{d}r^{*}\nonumber\\&=\frac{-8\textit{Pe}_{p}}{{\alpha _{i}}^{2}}-\frac{\textit{Pe}_{e}}{1-\eta }\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right) .\end{align}

The substitution of (4.8) into (4.7) then provides the axial mean of the solute zone as a function of the eigenfunction-weighted area average of arbitrary initial conditions, $\langle c _{0}^{*}(r^{*},0)f_{i}\rangle$ .

In this work, we discuss the term $\langle c_{0}^{*}(r^{*},0)f_{i}\rangle$ for two interesting initial conditions. The first initial condition we will consider is the case of solute that is distributed uniformly across the cross-section of the tube. For such initial conditions

(4.9) \begin{equation} \langle c _{0}^{*}(r^{*},0)f_{i}\rangle \propto \int\limits_{0}^{1}r^{*}\frac{J_{0}\!\left(\alpha _{i}r^{*}\right)}{J_{0}\!\left(\alpha _{i}\right)}\mathrm{d} r^{*}=0, \end{equation}

so (4.7) reduces to simply $M_{1}(\tau )=\tau \textit{Pe}$ (where $P\hspace{0pt}e$ is a signed quantity or zero). The latter reflects the fact that the solute zone (eventually) moves at the net bulk velocity of the flow. The second initial condition we will consider is an initial delta distribution of solute at $r^{*}=0$ . From the definition of the delta distribution, $\langle c _{0}^{*}(r^{*},0)f_{i}\rangle =2J_{0}(0)/J_{0}(\alpha _{i})$ , whence

(4.10) \begin{equation} M_{1}(\tau )=\tau \textit{Pe}+\sum _{i}^{\infty }\frac{2}{J_{0}\!\left(\alpha _{i}\right)}\!\left(\frac{-8\textit{Pe}_{p}}{{\alpha _{i}}^{2}}-\frac{\textit{Pe}_{e}}{1-\eta }\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)\right)\left(\frac{1-e^{-{{\alpha _{i}}^{2}}\tau }}{{\alpha _{i}}^{2}}\right) .\end{equation}

In the limit as $\tau \rightarrow \infty$ , (4.10) again reflects the fact that the solute zone ultimately moves at bulk velocity. Also, (4.10) shows that the mean of the solute cloud will in general deviate from that predicted by our quasi-steady solution. In most practical applications, this deviation is small and negligible. Further, in contrast to our quasi-steady solution, the mean axial location of the solute as predicted by this MoM formulation depends explicitly on the initial condition of solute. As per (4.9), this dependence vanishes for initial distributions of solute which are uniform across the cross-section of the tube, simplifying the problem of solving for higher-order moments of the solute zone. Given this, we shall only consider such initial conditions in our solution of the variance of the solute zone.

Before continuing we explicitly write here the quasi-steady limits of (4.7) and (4.10) for the case of pure PDF. Substituting $\textit{Pe}_{p}=\textit{Pe}$ and transforming into a moving frame given by $M_{1}^{\prime}(\tau )=M_{1}(\tau )-\tau \textit{Pe}$ , the values of $M_{1}(\tau )$ specified by (4.7) and (4.10) each reduce, respectively, to

(4.11) \begin{align} \lim _{\tau \rightarrow \infty }M_{1}^{\prime}(\tau )&=-8P\hspace{0pt}e\,\lim _{\tau \rightarrow \infty }\sum _{i}^{\infty }{\alpha _{i}}^{-4}\!\left\langle c _{0}^{*}(r^{*},0)f_{i}\right\rangle \!\big(1-e^{-{{\alpha _{i}}^{2}}\tau }\big)\nonumber\\&=-8\textit{Pe}\sum _{i}^{\infty }{\alpha _{i}}^{-4}\!\left\langle c _{0}^{*}(r^{*},0)f_{i}\right\rangle, \end{align}

and

(4.12) \begin{equation} \lim _{\tau \rightarrow \infty }M_{1}^{\prime}(\tau )=-16P\hspace{0pt}e\underset{\tau \rightarrow \infty }{\enspace \lim }\sum _{i}^{\infty }{\alpha _{i}}^{-4}\frac{1}{J_{0}\!\left(\alpha _{i}\right)}\!\Big(1-e^{-{{\alpha _{i}}^{2}}\tau }\Big)=-16\textit{Pe}\sum _{i}^{\infty }{\alpha _{i}}^{-4}\frac{1}{J_{0}\!\left(\alpha _{i}\right)} .\end{equation}

Equations (4.11) and (4.12) agree with Aris’ (Reference Aris1956) solution for the mean axial location of the solute. Specifically, (4.11) is equivalent to Aris’ (21) and (4.12) is equivalent to Aris’ (32) for an initial delta distribution of solute at $r^{*}=0$ .

Now that we have a solution for the first moment of solute, we can proceed with the separation of variables solution to ‘build up’ our solution for the second moment of solute, as previously mentioned. Accordingly, Barton (Reference Barton1983) further showed that the second moment of the solute distribution about the axial mean (again, for arbitrary geometries and velocity fields) behaves as follows:

(4.13) \begin{equation} v_{2}(\tau )=M_{2}-\left(M_{1}\right)^{2}=v_{2}\!\left(0\right)+2\left(1+\sum _{i}^{\infty }\frac{\left\langle \chi f_{i}\right\rangle ^{2}}{\mu _{i}}\right)\tau -2\sum _{i}^{\infty }\!\left\langle \chi f_{i}\right\rangle ^{2}\!\left(\frac{1-e^{-{\mu _{i}}\tau }}{{\mu _{i}}^{2}}\right) .\end{equation}

For the current work, $\langle \chi f_{i}\rangle$ and $\mu _{i}$ are given by (4.8) and (4.6), respectively. Equation (4.13) applies to solute that is initially distributed uniformly across the cross-section of the tube. Useful in comparing our MoM and quasi-steady formulations, we also define the transient effective dispersion coefficient to be half the rate of change of variance

(4.14) \begin{equation} \frac{1}{2}\frac{\mathrm{d}\!\left(v_{2}(\tau )\right)}{\mathrm{d}\tau }=D_{\textit{eff}}^{*\textit{trns}}=1+\sum _{i}^{\infty }\!\left(\frac{\left\langle \chi f_{i}\right\rangle ^{2}}{\mu _{i}}-\frac{\left\langle \chi f_{i}\right\rangle ^{2}e^{-{\mu _{i}}\tau }}{\mu _{i}}\right) .\end{equation}

(We note that Barton neither discusses this formulation nor makes this comparison.) Equation (4.14) provides the transient effective dispersion coefficient as a function of four dimensionless parameters ( $\textit{Pe}_{p}, \textit{Pe}_{e}, \tau$ and $\phi$ ), and is also useful to the case of perfectly opposed EOF and PDF. However, because (4.14) is a function of four parameters, it is difficult to present clear visualisations of different regimes of the solution. To remedy this, we will again formulate the dispersion coefficient in terms of the parameter $\beta$ and a Péclet number based on the sum of the EOF and PDF bulk velocities, $\textit{Pe}=(\langle u_{p}\rangle +\langle u_{e}\rangle )a/D$ . Then, letting $\hat{u}(r^{*})=u(r^{*})/\langle u\rangle$ , we write

(4.15) \begin{equation} \left\langle \chi f_{i}\right\rangle =\textit{Pe}\left(\frac{-8\beta }{{\alpha _{i}}^{2}}-\frac{(1-\beta)}{1-\eta }\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)\right)=\textit{Pe} \langle \hat{u}\!f_{i}\rangle .\end{equation}

Substituting (4.15) into (4.14), we may now again define a bulk Péclet-independent dispersion parameter as

(4.16) \begin{equation} \frac{48}{\textit{Pe}^{2}}\!\big(D_{\textit{eff}}^{*\textit{trns}}-1\big)=48\sum _{i}^{\infty }\!\left(\frac{\langle \hat{u}\!f_{i} \rangle ^{2}}{\mu _{i}}-\frac{ \langle \hat{u}\!f_{i} \rangle ^{2}e^{-{\mu _{i}}\tau }}{\mu _{i}}\right) .\end{equation}

The right-hand side of (4.16) is a function only of the dimensionless parameters $\beta, \phi$ and $\tau$ . This makes it easier to visualise than (4.14). We use the effective dispersion coefficient provided by (4.16) in most subsequent figures. Further, (4.16) asymptotes to our quasi-steady solution of the normalised dispersion coefficient (given by (3.21)) as $\tau \rightarrow \infty$ (we will note this in the discussion of figure 11).

Figure 11.

Benchmark comparisons between Brownian dynamics simulations and the analytical solution of the early, transient growth of the effective dispersion coefficient normalised as $48\textit{Pe}^{-2}(D_{\textit{eff}}^{*\textit{trns}}-1)$ (non-solid curves) versus dimensionless time, $\tau =\textit{tD}/a^{2}$ . The non-solid curves show the analytical solutions as derived from the MoM. The open markers show Brownian dynamics simulations benchmarking of the solutions. Curves in $(a)$ are indexed by varying values of the tube radius scaled by Debye length, $\phi$ . Curves in $(b)$ are indexed by varying values of the fraction of bulk velocity caused by pressure, $\beta$ . The translucent horizontal line segments denote the analytical, quasi-steady state solutions.

We note that for the special case of pure PDF ( $\textit{Pe}_{p}=\textit{Pe}$ ), (4.13) reduces to Barton’s (Reference Barton1983) solution for the variance of solute in a cylindrical tube (Barton’s equation (5.4)). Further, again for pure PDF, (4.14) analytically reduces to Taghizadeh et al.’s (Reference Taghizadeh, Valdés-Parada and Wood2020) solution for the transient effective dispersion coefficient (Taghizadeh et al.’s equation (7.3)). Taghizadeh et al. analysed the dispersion of solute subject to purely PDF in a cylindrical tube across all time regimes. They applied a cross-sectional-averaging method and derived an analytical expression for the aforementioned transient effective dispersion coefficient, as well as an associated one-dimensional convection–diffusion equation (Taghizadeh et al.’s equation (6.6)). To find solutions to this area-averaged convection–diffusion equation, Taghizadeh et al. relied on a (numerical) finite element method. Lastly, for $\textit{Pe}_{p}=\textit{Pe}$ and in the limit as $\tau \rightarrow \infty$ , (4.14) reduces to Aris’ (Reference Aris1956) solution for the steady effective dispersion coefficient (Aris’ equation (26)). Hence, our solutions for the solute variance and effective dispersion coefficient agree with and reduce to those of Aris (Reference Aris1956), Barton (Reference Barton1983) and Taghizadeh et al. (Reference Taghizadeh, Valdés-Parada and Wood2020) for the special case of pure PDF.

Lastly, we offer a note regarding the second or third moments about the mean for the current velocity fields. For arbitrary initial conditions, such a solution requires the evaluation of the following integral:

(4.17) \begin{equation} \int\limits _{0}^{1}r^{*}I_{0}\!\left(\phi r^{*}\right)J_{0} (\alpha _{i}r^{*})^{2}\mathrm{d}r^{*} .\end{equation}

We know of no analytical solution to this integral.

4.1. Minimisation of variance for transporting solute over a fixed axial distance

We analytically explore the minimisation of variance for transporting solute over some fixed axial distance of interest. This treatment considers, for example, an engineering situation in a microfluidic device where a cloud of solute is transported a distance $L$ from one location in a microchannel to a second location over some characteristic time $\tau _{o}$ . The following analysis allows the user of this device to maximise the area-averaged concentration of the arriving solute by minimising dispersion. Importantly, and in contrast to the optimisation in § 3.3, this analysis applies to values of the characteristic transport time $\tau _{o}$ that are any duration relative to the characteristic time of transverse diffusion, $a^{2}/D$ . In particular, the following analysis is applicable to processes which are brief (in terms of $\tau _{o}$ ) compared with the transverse diffusion time. The fixed value of solute transport distance is here relevant (whereas the precise distance of transport was irrelevant in § 3.3). We are only interested in minimising variance for moving solute, so we preclude the case of $\langle u\rangle \approx 0$ . First, if $\phi$ may be varied, then it should be made as large as possible to help minimise variance (as mentioned in § 3.3). However, we henceforth assume $\phi$ to be finite and fixed. Because we assumed $\langle u\rangle \neq 0$ , it is convenient to optimise the velocity profile in terms of the parameters $\beta =\langle u_{p}\rangle /\langle u\rangle$ and $\textit{Pe}$ . Substituting (4.15) into (4.13), we first present the variance of the solute zone in terms of these parameters

(4.18) \begin{equation} v_{2}(\tau )=v_{2}\!\left(0\right)+2\left(1+\textit{Pe}^{2}\sum _{i}^{\infty }\frac{ \langle \hat{u}f_{i} \rangle ^{2}}{\mu _{i}}\right)\tau -2\textit{Pe}^{2}\sum _{i}^{\infty } \langle \hat{u}f_{i}\rangle ^{2}\!\left(\frac{1-e^{-{\mu _{i}}\tau }}{{\mu _{i}}^{2}}\right) .\end{equation}

Given this formulation, we will show that the value of $\beta$ required to minimise variance for transporting solute over a fixed distance (henceforth denoted $\beta _{o}^{\textit{trns}}(\tau _{o})$ ) depends only on the characteristic dimensionless time to transport solute over that given distance, $\tau _{o}$ , and the dimensionless parameter $\phi$ . Now, we solved the steady Stokes equation for our velocity profile. Hence, upon injection of solute into the tube, all parameters related to the velocity field (e.g. $\beta, \phi$ and $\textit{Pe}$ ) are immediately ‘set’ and cannot vary with time. For the minimisation, we differentiate (4.18) with respect to $\beta$ , substitute $\tau =\tau _{o}$ , and then set the resulting quantity equal to zero

(4.19) \begin{equation} \frac{\partial }{\partial \beta }\!\left(v_{2}\!\left(\tau _{o}\right)\right)=\sum _{i}^{\infty }\frac{\partial }{\partial \beta } \langle \hat{u}\!f_{i} \rangle ^{2}\!\left(\frac{\tau _{o}}{\mu _{i}}-\frac{1-e^{-{\mu _{i}}{\tau _{o}}}}{\mu _{i}^{2}}\right)=0 .\end{equation}

The only value dependent on $\beta$ in (4.19) is $\langle \hat{u}f_{i}\rangle ^{2}$ , so we evaluate as follows:

(4.20) \begin{align} \frac{\partial }{\partial \beta }\!\left(\langle \hat{u}f_{i} \rangle ^{2}\right)&=2\left(\frac{-8\beta }{{\alpha _{i}}^{2}}-\frac{(1-\beta)}{1-\eta }\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)\right)\left(\frac{-8}{{\alpha _{i}}^{2}}+\frac{1}{1-\eta }\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)\right)\nonumber\\&= 2\beta \left(\frac{64}{{\alpha _{i}}^{4}}-\frac{16}{{\alpha _{i}}^{2}\!\left(1-\eta \right)}\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)+\frac{1}{\left(1-\eta \right)^{2}}\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)^{2}\right)\nonumber\\&\quad +\frac{16}{{\alpha _{i}}^{2}\!\left(1-\eta \right)}\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)-\frac{2}{\left(1-\eta \right)^{2}}\!\left(\frac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)^{2}. \end{align}

We substitute (4.20) into (4.19) and rearrange to find

(4.21) \begin{align} & \beta^{trns}_{o}(\tau_{o}) \nonumber \\& \quad =\frac{\sum\limits _{i}^{\infty }\dfrac{\left(1-\mu _{i}\tau _{o}-e^{-{\mu _{i}}{\tau _{o}}}\right)}{\mu _{i}^{2}}\!\left(\dfrac{8}{{\alpha _{i}}^{2}}\!\left(\dfrac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)-\dfrac{1}{\left(1-\eta \right)}\!\left(\dfrac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)^{2}\right)}{\sum\limits _{i}^{\infty }\dfrac{\left(\mu _{i}\tau _{o}-1+e^{-{\mu _{i}}{\tau _{o}}}\right)}{\mu _{i}^{2}}\!\left(\dfrac{64\left(1-\eta \right)}{{\alpha _{i}}^{4}}-\dfrac{16}{{\alpha _{i}}^{2}}\!\left(\dfrac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)+\dfrac{1}{\left(1-\eta \right)}\!\left(\dfrac{\eta \phi ^{2}}{{\alpha _{i}}^{2}+\phi ^{2}}\right)^{2}\right)} .\end{align}

The latter expression provides the unique value of $\beta$ that minimises the dispersion across all time regimes. Equation (4.21) depends only on $\tau _{o}$ and $\phi$ . For example, if we specify the value of $\tau _{o}$ as the characteristic time to transport solute through a tube, then (4.21) can be used to find the fixed value of $\beta _{o}^{\textit{trns}}(\tau _{o})$ that minimises variance for such transport. We note that (4.21) asymptotes to $\beta _{o}$ of our quasi-steady solution as $\tau _{o}\rightarrow \infty$ (see figure 12).

Figure 12.

Deviation of the transient value of the variance-minimising fraction of bulk velocity caused by pressure, $\beta _{o}^{\textit{trns}}(\tau _{o})$ , from the quasi-steady state value, $\beta _{o}$ , versus the characteristic dimensionless time of interest, $\tau _{o}$ . The curves are indexed by the tube radius scaled by Debye length, $\phi$ . The inset is a closeup for the interval of $0\lt \tau _{o}\lt 0.1$ . Dispersion caused by predominantly EOF with very thin EDLs (e.g. $\phi \gt 50$ ) very quickly asymptotes to quasi-steady state dynamics.

Next, we are also interested in deriving the optimum Péclet number $P\hspace{0pt}e_{o}^{\textit{trns}}$ to minimise the variance for transporting the solute over some fixed dimensionless axial distance of interest, $L^{*}$ . To this end, we substitute $\tau =L^{*}/\textit{Pe}$ into (4.18), then set the derivative with respect to $P\hspace{0pt}e$ equal to zero

(4.22) \begin{align} &\frac{\partial \left(v_{2}\right)}{\partial \left(\textit{Pe}\right)}=-\frac{L^{*}}{\left(\textit{Pe}_{o}^{\textit{trns}}\right)^{2}}\nonumber\\&\quad+\sum _{i}^{\infty }\!\left(\frac{L^{*}\!\langle \hat{u}f_{i}\rangle ^{2}}{\mu _{i}}-\frac{2\textit{Pe}_{o}^{\textit{trns}} \langle \hat{u}f_{i}\rangle ^{2}}{\mu _{i}^{2}}+\frac{2\textit{Pe}_{o}^{\textit{trns}} \langle \hat{u}f_{i}\rangle ^{2}e^{\frac{-\mu _{i}L^{*}}{\textit{Pe}_{o}^{\textit{trns}}}}}{\mu _{i}^{2}}+\frac{L^{*} \langle \hat{u}f_{i}\rangle ^{2}e^{\frac{-\mu _{i}L^{*}}{\textit{Pe}_{o}^{\textit{trns}}}}}{\mu _{i}}\right)=0 .\end{align}

Multiplying by $(P\hspace{0pt}e_{o}^{\textit{trns}})^{2}$ and rearranging, (4.22) can be expressed as

(4.23) \begin{equation} -L^{*}+\sum _{i}^{\infty }\!\left(\!L^{*}\!\left(P\hspace{0pt}e_{o}^{\textit{trns}}\right)^{2}\frac{\langle \hat{u}f_{i}\rangle ^{2}}{\mu _{i}}\!\left(\!1+e^{\frac{-\mu _{i}L^{*}}{\textit{Pe}_{o}^{\textit{trns}}}}\right)+2\left(P\hspace{0pt}e_{o}^{\textit{trns}}\right)^{3}\frac{ \langle \hat{u}f_{i}\rangle ^{2}}{\mu _{i}^{2}}\!\left(e^{\frac{-\mu _{i}L^{*}}{\textit{Pe}_{o}^{\textit{trns}}}}-1\right)\!\right)=0 .\end{equation}

We will now explain the process to find the optimum combination of $\beta$ and $\textit{Pe}$ to minimise variance for transporting solute over a fixed axial distance (again assuming $\phi$ to be fixed). Recall from (4.15) that $\langle \hat{u}f_{i}\rangle$ and therefore (4.23) are functions of $\beta$ . Thus, to find the global optimum value of $P\hspace{0pt}e_{o}^{\textit{trns}}$ for a fixed value of $L^{*}$ , we first substitute the $\beta _{o}^{\textit{trns}}(\tau _{o})$ function as given by (4.21) into (4.23). We then substitute $\tau _{o}=L^{*}/P\hspace{0pt}e_{o}^{\textit{trns}}$ into (4.23). Next, we solve the resulting equation for $P\hspace{0pt}e_{o}^{\textit{trns}}$ using the Newton–Raphson method (Engquist Reference Engquist2015, pp. 1023–1028) or some other root-finding method. Once we have $P\hspace{0pt}e_{o}^{\textit{trns}}$ for a given value of $L^{*}$ , we can again use the relationship $\tau _{o}=L^{*}/P\hspace{0pt}e_{o}^{\textit{trns}}$ to find the corresponding value of $\tau _{o}$ . Next, using (4.21), we compute the specific value of $\beta _{o}^{\textit{trns}}(\tau _{o})$ at that value of $\tau _{o}$ . Thus, this process gives the fixed global optimum values of $\textit{Pe}$ and $\beta$ to minimise variance for transporting solute over a dimensionless distance $L^{*}$ . Lastly, per figure 13, we note that $P\hspace{0pt}e_{o}^{\textit{trns}}$ asymptotes to $\textit{Pe}_{o}$ as $L^{*}\rightarrow \infty$ .

Figure 13.

Solutions related to the transient variance-minimising Péclet number, $\textit{Pe}_{o}^{\textit{trns}}$ . In $(a), \textit{Pe}_{o}^{\textit{trns}}$ is plotted as a function of the dimensionless axial distance of travel, $L^{*}$ . The non-solid curves indicate transient solutions from the MoM. The translucent horizontal line segments denote the associated quasi-steady state solutions. The quasi-steady and transient solutions agree asymptotically after approximately $L^{*}=100$ . $(b)$ Variance, $v_{2}(\tau _{o}=L^{*}/\textit{Pe})$ , of the solute zone after travelling a dimensionless axial distance of $L^{*}=100$ as a function of Péclet number, $P\hspace{0pt}e$ . The red dots indicate minima of $v_{2}(\tau _{o})$ as determined from the solution of $\textit{Pe}_{o}^{\textit{trns}}$ from (4.23). For both panels, curves are indexed by the tube radius scaled by Debye length, $\phi$ . Figure 13 was produced by setting $\tau _{o}=L^{*}/\textit{Pe}$ and $\beta =\beta _{o}^{\textit{trns}}(\tau _{o})$ .

6. Summary and conclusions

We performed a Taylor-type analysis for the long-term dispersive behaviour of a neutral solute subjected to simultaneous PDF and EOF in a long, thin cylindrical tube. We first derived a formulation for the EOF velocity profile, valid for an EDL of arbitrary thickness relative to the channel radius. We then performed a classic scaling analysis of the convection diffusion equation, demonstrating that axial dispersion due to the non-uniform velocity field must be in approximate balance with radial diffusion, which allowed us to derive the coefficient of effective dispersion. In the process, we also derived the radial deviation from the area-averaged concentration field. This allowed us to describe the full, two-dimensional concentration field (including the radial component) as a function of the initial conditions of the solute in long-time regimes. To our knowledge, this is the first full description of the concentration field, with the radial component, for combined EOF and PDF. Ours is also, to our knowledge, the first published derivation of the dispersion coefficient for these flows.

Next, we described the behaviour of variance growth and presented analytical methods to optimise the flow Péclet number and combination of EOF and PDF required to minimise the rate of growth of the variance (for some finite axial length travelled). We demonstrated that, to minimise the rate of growth of the variance, it is optimal to oppose EOF with progressively stronger PDF as the EDLs become thicker. Additionally, we showed that the optimum Péclet number tends to be inversely related to the effective dispersion of a normalised flow regime (e.g. thin-EDL EOF with unit bulk velocity tends to have a large optimum Péclet number whereas unit PDF has a low optimum Péclet number). Lastly, we presented the unsteady, two-dimensional concentration field for the specific example initial conditions of a ‘top hat’ of solute and an initial delta distribution of solute. The major features of our analytical solutions were successfully benchmarked with Brownian dynamics simulations. The latter included the dispersion coefficient, the two-dimensional concentration field, the radial solute distribution and the optimum value of $\beta$ to minimise dispersion for a wide range of cases.

Next, we analysed the transient dispersive behaviour of the solute using the MoM. The MoM maintains two key advantages over Taylor-type analysis in that solutions are valid across all time regimes and the width of the solute cloud (including its initial width) can be any size relative to the inner radius of the tube. We derived analytical solutions for the axial mean and variance of the solute zone, allowing us to analyse the early behaviour of the solute. We demonstrated that the transient effective dispersion coefficient asymptotes to the quasi-steady effective dispersion coefficient. We presented analytical methods to minimise the variance for transporting the solute over a fixed axial distance. As expected, the transient parameters (including the $P\hspace{0pt}e_{o}^{\textit{trns}}$ and $\beta _{o}^{\textit{trns}}(\tau _{o})$ functions) required to minimise dispersion asymptote to the corresponding asymptotic values predicted by the quasi-steady, Taylor-type analysis. The transient dispersion coefficient, $D_{\textit{eff}}^{*\textit{trns}}$ , was successfully benchmarked with Brownian dynamics simulations for a wide range of cases. In all tested cases, our analytical models (both Taylor-type and MoM formulations) showed excellent agreement with Brownian dynamics simulations.

One limitation of this work is that it is restricted to physical regimes with a sufficiently low zeta potential such that the Debye–Hückel approximation applies. Further work might treat the case of high zeta potential numerically (or using an analytical approximation). An additional area for future work might also analyse early (transient) variance growth for solutes whose initial distribution is non-uniform across the cross-section of the tube.

In conclusion, this work presented benchmarked (compared with Brownian dynamics simulations) predictions of the behaviour of a neutral solute subjected to simultaneous EOF and PDF in all time regimes. Additionally, we presented analytical methods to minimise hydrodynamic dispersion in a wide range of physical regimes.