1. Introduction
Advective dispersion of solutes in long thin tubes is important to the design and optimization of a wide range of devices, from chemical processing and separations systems, to microfluidic chips (Brenner & Edwards Reference Brenner and Edwards1993). G.I. Taylor (Reference Taylor1953) first reported a closed-form solution for the long-term dispersive behaviour of a solute in a circular cylindrical tube driven by a fully developed Poiseuille flow. His original solution considered a regime where axial molecular diffusion (molecular diffusion in the axial direction) is negligible compared to the dispersion effect, but radial molecular diffusion continues to play a role because of the radial concentration gradient from the convection. In terms of inequalities, this original analysis is applicable for 
$L/a\gg Pe_a\gg 6.9$, where 
$L$ is the characteristic channel length, 
$a$ is the radius of the tube, and 
$Pe_a$ is a Péclet number based on 
$a$. Aris (Reference Aris1956) later included the effects of axial molecular diffusion of the solute and introduced the method of moments, which enabled treatment of all stages of the dispersive process in long thin channels with fairly arbitrary cross-sections. Dispersion analyses, including molecular diffusion and dynamic regimes where radial diffusion times are much shorter than advection times, are typically termed Taylor–Aris analyses. Since these seminal papers, there has been much work analysing the dispersion behaviour in various systems. Brenner & Edwards (Reference Brenner and Edwards1993) includes a broad range example analyses, including dispersion of flows through porous media and effects of chemical reactions and surface adsorption.
A significant emphasis has been placed on the dispersion problem in spatially-periodic channel geometries. The seminal work of H. Brenner (Reference Brenner1980) (which is later referred as Brenner–Aris theory, or the macrotransport paradigm) provided a generalized framework to predict the long-term dispersion behaviour of point-size particles in any spatially periodic channels or porous media. This approach solves an elliptical partial differential equation of a cell field 
$\boldsymbol {B}$ (commonly referred as the 
$\boldsymbol {B}$ field) defined on a periodic unit cell, and uses this to compute a long-term dispersion tensor. Hoagland & Prud'Homme (Reference Hoagland and Prud'Homme1985) applied Brenner–Aris theory to the dispersion in a long thin, axisymmetric channel with a sinusoidal variation of radius along the streamwise direction. They presented numerical solutions of the method of moments for arbitrary wavelengths, and they also derived an analytical expression for the long-term dispersion of solute in the limit of long wavelength, which compared well with their numerical model. Bolster, Dentz & Le Borgne (Reference Bolster, Dentz and Le Borgne2009) performed an analysis similar to that of Hoagland & Prud'Homme (Reference Hoagland and Prud'Homme1985) in axisymmetric channel with a sinusoidal variation of radius, and extended the results by comparing the predictions with Brownian dynamics simulations. Dorfman and his coworkers also used Brenner's method but for the advective component considered electrophoretic and electro-osmotic flows with finite double layers (Yariv & Dorfman Reference Yariv and Dorfman2007; Dorfman Reference Dorfman2008, Reference Dorfman2009). Later, Adrover, Venditti & Giona (Reference Adrover, Venditti and Giona2019) considered dispersion of both point and finite-sized particles in a long thin sinusoidal channel. For point particles, they obtained various asymptotic expressions for the effective long-term dispersion coefficient in the limits of large and small 
$Pe_a$.
Methods other than Brenner's theory were also proposed to analyse the dispersion in periodic channels. For example, Rosencrans (Reference Rosencrans1997) followed the centre manifold theory of Mercer & Roberts (Reference Mercer and Roberts1990) and derived the dispersion coefficient for periodic curved channels. They hypothesized that their analysis extended to strong variations in geometry but did not support this with numerical analyses. In a series of papers, Martens and his coworkers explored the dispersion process in periodic channels from the (Lagrangian) perspectives of individual particles and formulated dispersion theory based on the Fick–Jacobs equation (Martens et al. Reference Martens, Schmid, Schimansky-Geier and Hänggi2011, Reference Martens, Schmid, Straube, Schimansky-Geier and Hänggi2013a,Reference Martens, Straube, Schmid, Schimansky-Geier and Hänggib, Reference Martens, Straube, Schmid, Schimansky-Geier and Hänggi2014).
Several studies also addressed the dispersion problem in spatially non-periodic systems. Early work from Saffman (Reference Saffman1959) analysed the long-term longitudinal dispersion in the direction of mean flow in a random porous media where the flow satisfies Darcy's law. Gill, Ananthakrishnan & Nunge (Reference Gill, Ananthakrishnan and Nunge1968) analysed analytically and numerically the dispersion in a velocity entrance region where the flow field is not fully developed. Gill & Güceri (Reference Gill and Güceri1971) later investigated numerically the dispersion in diverging channels over a wide range of angles of divergence and Péclet number. Smith (Reference Smith1983) described the longitudinal dispersion in a varying channel in terms of the memory effect of the dispersion coefficient, and also investigated the changes in dispersion coefficients in several profiles that are pertinent to geophysics, such as sudden changes in breadth, centrifugally driven secondary flow associated with the curvature in the flow path, and changes in the depth profile. Mercer & Roberts (Reference Mercer and Roberts1990) used centre manifold theory to derive the spatially dependent dispersion coefficient in channels with varying flow properties. Although they claimed that their approach may be applicable to time-dependent flow and variable diffusivity, the dispersion coefficient alone can be misleading in predicting actual dispersion behaviour. Horner, Arce & Locke (Reference Horner, Arce and Locke1995) studied the Taylor dispersion in an axisymmetric system with two linear profiles of radius changes, which was intended to approximate the convergent and divergent sections of stenosis sites in arteries. Bryden & Brenner (Reference Bryden and Brenner1996) analysed the Taylor–Aris dispersion in a diverging conical channel and a flared, axisymmetric Venturi tube geometry. They used multiple time scale analysis to obtain an asymptotic expression for the effective dispersion coefficient, and derived a partial differential equation (PDE) for a conditional probability density function describing the solute. However, the latter work neither provided a solution for this PDE nor analysed the long-term dispersion and growth of a solute zone. Stone & Brenner (Reference Stone and Brenner1999) also investigated the dispersion in a different spatially non-periodic system. This work considered dispersion in a radially expanding flow between large parallel plates. That flow results in a mean velocity that varies in the streamwise direction and therefore exhibits a location-dependent dispersion coefficient.
Note that all the aforementioned literature considers the spatial dispersion problem, as the dispersion process is quantified by the spatial moments of the averaged solute distribution. An alternative framework of temporal dispersion problem, proposed by Danckwerts (Reference Danckwerts1953), quantifies the dispersion process in terms of its temporal moments. This formulation is particularly useful when the solute is monitored at a given distance from the inlet, and can also be useful in analysing and benchmarking results from numerical models (Rodrigues Reference Rodrigues2021). Although the formulations are different, the spatial moments and temporal moments of the dispersion process have been shown to be interrelated (Vikhansky & Ginzburg Reference Vikhansky and Ginzburg2014; Ginzburg & Vikhansky Reference Ginzburg and Vikhansky2018).
To our knowledge, all previous analyses focused on the long-term dispersion behaviour and were not able to provide a simple prediction of both the short-term and long-term spatial evolution of solute. We also know of no analytical work towards Taylor–Aris dispersion in an arbitrarily shaped axisymmetric channel. Such analyses have significant potential to influence the design and analyses of a wide range of systems, including microfluidic devices. In the current study, we analyse the Taylor–Aris dispersion in an arbitrarily shaped axisymmetric channel with a slowly varying radius. We derive an expression for solute dynamics in terms of two coupled ordinary differential equations (ODEs). These two ODEs enable prediction of the time evolution of the solute zone based on channel geometry and the assumed lubrication flow. We compare and benchmark our predictions with Brownian dynamics analysis. We further develop a formulation useful in inverse problems where channels are designed to achieve desired spatial distribution of solute variance. We demonstrate this for the design of two channel geometries of specialized function. The first geometry results in a constant axial variance, and the second results in a sinusoidal axial variance in space.
2. Theory
2.1. Taylor–Aris dispersion in arbitrarily shaped cylinders
We analyse Taylor–Aris dispersion for flow in an axisymmetric channel with a slowly varying, arbitrary radius 
$a = a(x)$ (figure 1). We define variables for the first and second derivatives as 
$\beta (x)=\text {d}a/\text {d}x$ and 
$\gamma (x)=\text {d}^2a/\text {d}x^2$. Given the azimuthal symmetry, the concentration 
$c(x,r,t)$ of a solute evolves according to the following convection–diffusion equation in cylindrical coordinates:
\begin{equation} \dfrac{\partial c}{\partial t} + u_{r}\,\frac{\partial c}{\partial r} + u_{x}\,\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). \end{equation}For the velocity field, the Navier–Stokes equation is solved with the assumptions typical of lubrication theory. That is, we consider a slowly varying radius with a characteristic radius much smaller than the axial distance of interest (Langlois & Deville Reference Langlois and Deville1964, pp. 229–240). Hence we write the velocity field as
\begin{equation} \left. \begin{aligned} u_{x}(x,r) & =2U(x)\left(1-\frac{r^2}{a(x)^2}\right)+{O}(U_0\,Re\,\epsilon^2, U_0\epsilon^2),\\ u_r(x,r) & =2\beta(x)\,U(x)\left(\frac{r}{a(x)}-\frac{r^3}{a(x)^3}\right)+{O}(U_0\,Re\,\epsilon^3, U_0\epsilon^3), \end{aligned} \right\} \end{equation}
where 
$Re$ is the Reynolds number, and 
$\epsilon$ is our primary smallness parameter defined as the ratio between the characteristics radius 
$a_0$ (i.e. the radius at 
$x=0$) and the characteristic axial length scale of radius variation 
$\lambda$ (table 1). The second terms on the right-hand side indicate the order of magnitude of the truncation error associated with this asymptotic approximation. Here, 
$U_0$ is the characteristic mean axial velocity defined as the mean axial velocity at 
$x=0$. We also assume that the characteristic width of the solute zone 
$\sigma$ is much smaller than the characteristic wavelength of radius variation, or 
$\sigma /\lambda \ll 1$. For the Taylor–Aris analysis, we follow the notation of Stone & Brenner (Reference Stone and Brenner1999) and expand each variable into cross-sectional averages of the form 
$\langle ({\cdot })\rangle \equiv ({1}/{{\rm \pi} \, a(x)^2})\int _{0}^{a(x)}2{\rm \pi} r({\cdot })\,\textrm {d}r$ and deviations therefrom defined as 
$({\cdot })^{\prime } \equiv ({\cdot }) - \langle ({\cdot })\rangle$. Whence the area-averaged velocity components are
\begin{gather} \langle u_{x}\rangle=U(x),\quad \langle u_r\rangle=\frac{8}{15}\,\beta(x)\,U(x), \end{gather}
\begin{gather} u_x^{\prime}=U(x)\left(1-\frac{2r^2}{a(x)^2}\right),\quad u_r^{\prime}=2\beta(x)\,U(x)\left(\frac{r}{a(x)}-\frac{r^3}{a(x)^3}-\frac{4}{15}\right). \end{gather}We then expand the convective–diffusion equation as
$$\begin{gather} \frac{\partial \langle c\rangle}{\partial t}+\frac{\partial c^{\prime}}{\partial t} + \langle u_{r}\rangle\,\frac{\partial c^{\prime}}{\partial r} + u_{r}^{\prime}\,\frac{\partial c^{\prime}}{\partial r} + U(x)\,\frac{\partial \langle c\rangle}{\partial x} + U(x)\,\frac{\partial c^{\prime}}{\partial x} + u_{x}^{\prime}\,\frac{\partial \langle c\rangle}{\partial x} + u_{x}^{\prime}\,\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 \langle c\rangle}{\partial x^2} + \frac{\partial^2 c^{\prime}}{\partial x^2}\right) . \end{gather}$$The area average of the latter equation is then
$$\begin{gather} \frac{\partial \langle c\rangle}{\partial t} + \left\langle u_{r}\right\rangle\left\langle\frac{\partial c^{\prime}}{\partial r}\right\rangle + \left\langle u_{r}^{\prime}\,\frac{\partial c^{\prime}}{\partial r}\right\rangle + U(x)\,\frac{\partial \langle c\rangle}{\partial x} + U(x)\left\langle\frac{\partial c^{\prime}}{\partial x}\right\rangle + \left\langle u_{x}^{\prime}\,\frac{\partial c^{\prime}}{\partial x}\right\rangle \nonumber\\ = D\left(\left\langle\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right)\right\rangle + \frac{\partial^2 \langle c\rangle}{\partial x^2} + \left\langle\frac{\partial^2 c^{\prime}}{\partial x^2}\right\rangle\right) . \end{gather}$$We simplify four of the terms in (2.6) using Leibniz's rule, integration by parts and the chain rule, as follows:
\begin{align} \left\langle\frac{\partial c^{\prime}}{\partial r}\right\rangle &= \frac{1}{a(x)^2}\int_{0}^{a(x)}2r\,\frac{\partial c^{\prime}}{\partial r}\,\mathrm{d}r=\frac{2}{a(x)^2}\int_{0}^{a(x)}\left(\frac{\partial}{\partial r}(rc^{\prime})-c^{\prime}\right)\,\mathrm{d}r\nonumber\\ &=\frac{2}{a(x)^2}\left(a\left.c^{\prime}\right\vert_{r=a}-\int_{0}^{a}c^{\prime}\,\mathrm{d}r\right), \end{align}
\begin{align} \left\langle\frac{\partial c^{\prime}}{\partial x}\right\rangle &= \frac{1}{a(x)^2}\int_{0}^{a(x)}2r\,\frac{\partial c^{\prime}}{\partial x}\,\mathrm{d}r = \frac{2}{a(x)^2}\int_{0}^{a(x)}\frac{\partial (rc^{\prime})}{\partial x}\,\mathrm{d}r\nonumber\\ &=\frac{2}{a(x)^2}\left[\frac{\mathrm{d}}{\mathrm{d} x}\int_{0}^{a(x)}rc^{\prime}\,\mathrm{d}r-a\left.c^{\prime}\right\vert_{r=a}\beta(x)\right]=\frac{-2\left.c^{\prime}\right\vert_{r=a}\beta(x)}{a(x)}, \end{align}
\begin{align} \left\langle\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right)\right\rangle &= \frac{2}{a(x)^2}\int_{0}^{a(x)}r\,\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right)\,\mathrm{d}r=\frac{2}{a(x)^2}\left(a\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a} - 0\right)\nonumber\\ &= \frac{2}{a}\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a}, \end{align}
\begin{align} \left\langle\frac{\partial^2 c^{\prime}}{\partial x^2}\right\rangle &= \frac{2}{a(x)^2}\int_{0}^{a(x)}r\,\frac{\partial^2c^{\prime}}{\partial x^2}\,\mathrm{d}r=\frac{2}{a(x)^2}\int_{0}^{a(x)}\frac{\partial^2 (rc^{\prime})}{\partial x^2}\,\mathrm{d}r\nonumber\\ &=\frac{2}{a^2}\left[\frac{\mathrm{d}}{\mathrm{d} x}\int_{0}^{a(x)}\frac{\partial (rc^{\prime})}{\partial x}\,\mathrm{d}r - \beta(x)\,\frac{\partial\left(a\left.c^{\prime}\right\vert_{r=a}\right)}{\partial x}\right]\nonumber\\ &=\frac{2}{a^2}\left[\frac{\mathrm{d}}{\mathrm{d} x}\biggl(\frac{\mathrm{d}}{\mathrm{d} x}\int_{0}^{a(x)}rc^{\prime}\,\mathrm{d}r -\beta(x)\,a\left.c^{\prime}\right\vert_{r=a}\biggr) - \beta(x)\,\frac{\partial\left(a\left.c^{\prime}\right\vert_{r=a}\right)}{\partial x}\right]\nonumber\\ &={-}\frac{2\gamma \left.c^{\prime}\right\vert_{r=a}}{a} - \frac{4\beta}{a^2}\,\frac{\partial}{\partial x}\left(a\left.c^{\prime}\right\vert_{r=a}\right). \end{align}Inserting these terms into (2.6), we have
$$\begin{gather} \frac{\partial \langle c\rangle}{\partial t} + \frac{2\langle u_{r}\rangle}{a^2}\left(a\left.c^{\prime}\right\vert_{r=a}-\int_{0}^{a}c^{\prime}\,\mathrm{d}r\right) + \left\langle u_{r}^{\prime}\,\frac{\partial c^{\prime}}{\partial r}\right\rangle + U(x)\,\frac{\partial \langle c\rangle}{\partial x} -\frac{2U(x)\left.c^{\prime}\right\vert_{r=a}\beta(x)}{a(x)} + \left\langle u_{x}^{\prime}\,\frac{\partial c^{\prime}}{\partial x}\right\rangle \nonumber\\ = \frac{2D}{a}\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a} + D\,\frac{\partial^2 \langle c\rangle}{\partial x^2} - \frac{4\beta D}{a^2}\,\frac{\partial}{\partial x}\left(a\left.c^{\prime}\right\vert_{r=a}\right)-\frac{2\gamma D\left.c^{\prime}\right\vert_{r=a}}{a}. \end{gather}$$Subtracting (2.11) from (2.5) yields the following equation for the deviation of concentration:
\begin{align} &\frac{\partial c^{\prime}}{\partial t} + \langle u_{r}\rangle\,\frac{\partial c^{\prime}}{\partial r} - \frac{2\langle u_{r}\rangle}{a^2}\left(a\left.c^{\prime}\right\vert_{r=a}-\int_{0}^{a}c^{\prime}\,\mathrm{d}r\right) + u_{r}^{\prime}\,\frac{\partial c^{\prime}}{\partial r} - \left\langle u_{r}^{\prime}\,\frac{\partial c^{\prime}}{\partial r}\right\rangle\nonumber\\ &\qquad + U(x)\,\frac{\partial c^{\prime}}{\partial x} + \frac{2U(x)\left.c^{\prime}\right\vert_{r=a}\beta(x)}{a(x)} + u_{x}^{\prime}\,\frac{\partial \langle c\rangle}{\partial x} + u_{x}^{\prime}\,\frac{\partial c^{\prime}}{\partial x} - \left\langle u_{x}^{\prime}\,\frac{\partial c^{\prime}}{\partial x}\right\rangle\nonumber\\ &\quad = D\,\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right) - \frac{2D}{a}\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a} + D\,\frac{\partial^2 c^{\prime}}{\partial x^2} + \frac{4\beta D}{a^2}\,\frac{\partial}{\partial x}\left(a\left.c^{\prime}\right\vert_{r=a}\right)+\frac{2\gamma D\left.c^{\prime}\right\vert_{r=a}}{a}. \end{align}
To determine the dominant terms on the left- and right-hand sides of the equation, we perform a scaling analysis with the following scales: 
$\langle c\rangle =c_0\langle c^{*}\rangle$, 
$c^{\prime }=c_0^{\prime }c^{\prime *}$, 
$x=\sigma x^{*}$, 
$r=a_0r^{*}$, 
$\eta =a_0/\sigma$, 
$t=t_{{obs}}t^{*}=({\sigma }/{U_0})t^{*}$, 
$u_x^{\prime }=U_0u_x^{\prime *}$ and 
$u_r=\epsilon U_0u_r^{*}$. We define 
$t_{{obs}}$ as the characteristic time over which the solute is observed. This observation time scale is analogous to the advective time scale of traditional Taylor–Aris theory. Here, 
$c_0$ and 
$c_0^{\prime }$ are the characteristic scales of the area-averaged solute concentration and its deviation, respectively. As summarized in table 1, our scaling assumes four smallness parameters: 
$\epsilon =a_0/\lambda \ll 1$, 
$\eta =a_0/\sigma \ll 1$, 
$\sigma /\lambda \ll 1$ and 
$c_0^{\prime }/c_0\ll 1$. We summarize the order of magnitude of each term in (2.12) in table 2.
Table 1. Summary of assumptions made in this study.
Table 2. Ranked summary of terms in (2.12).
Consistent with a Taylor–Aris dispersion regime, keeping the dominant advective dispersion term and the dominant radial diffusion term results in an approximate balance as
\begin{equation} u_{x}^{\prime}\,\frac{\partial \langle c\rangle}{\partial x}=D\,\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right) - \frac{2D}{a}\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a}. \end{equation}
Next, we consider the boundary condition on the inner wall of the channel. We impose a no-flux boundary condition at each axial location along the wall, 
$\boldsymbol {\nabla } c\boldsymbol {\cdot } \hat {n}=0$. This requires
\begin{equation} \left.{-\frac{\partial c}{\partial r}}\right\vert_{r=a}+\beta\left.\frac{\partial c}{\partial x}\right\vert_{r=a}={-}\left.\frac{\partial c^{\prime}}{\partial r}\right\vert_{r=a}+\beta\,\frac{\partial \langle c\rangle}{\partial x}+\beta\left.\frac{\partial c^{\prime}}{\partial x}\right\vert_{r=a}=0. \end{equation}Substituting the boundary condition into (2.13) and keeping only terms of the same order as the dominant terms in table 2, we have
\begin{equation} D\,\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\,\frac{\partial c^{\prime}}{\partial r}\right)\approx u_{x}^{\prime}\,\frac{\partial \langle c\rangle}{\partial x} + \frac{2D\beta}{a}\,\frac{\partial \langle c\rangle}{\partial x}. \end{equation}We next directly integrate (2.15) as
\begin{equation} r\,\frac{\partial c^{\prime}}{\partial r} = \frac{\partial \langle c\rangle}{\partial x}\,\frac{U(x)}{D}\left(\frac{r^2}{2} - \frac{r^4}{2a^2}\right) +\frac{\beta r^2}{a}\,\frac{\partial \langle c\rangle}{\partial x} + c_1(x,t). \end{equation}
Here, the function 
$c_1(x,t)$ results from the partial integration with respect to 
$r$. Evaluating this expression at 
$r=0$, we see 
$r\left.{\partial c^{\prime }}/{\partial r}\right \vert _{r=0}=c_1=0$, whence 
$c_1(x,t)=0$.Integrating the equation a second time, we obtain
\begin{equation} c^{\prime} = \frac{\partial \langle c\rangle}{\partial x}\,\frac{U(x)}{D}\left(\frac{r^2}{4} - \frac{r^4}{8a^2} \right) +\frac{\beta}{2a}\,r^2\,\frac{\partial \langle c\rangle}{\partial x} + c_2(x,t). \end{equation}
By the definition of the fluctuating quantity, 
$\int _0^{a}rc^{\prime }\,\mathrm {d}r=0$. We use this to solve for 
$c_2$ and obtain an expression for 
$c^{\prime }$,
\begin{equation} c^{\prime}(x,r,t) = \frac{\partial \langle c\rangle}{\partial x}\,\frac{U(x)}{D}\left(\frac{r^2}{4} - \frac{r^4}{8a^2} -\frac{a^2}{12}\right)+\left(\frac{\beta r^2}{2a} - \frac{1}{4}\,\beta a\right)\frac{\partial \langle c\rangle}{\partial x}, \end{equation}
with error of order 
$c_0\times {O}(Pe_{a_0}\,(c_0^{\prime }/c_0)\eta, Pe_{a_0}\,Re\epsilon ^2\eta, Pe_{a_0}\,\epsilon ^2\eta )$.
We next differentiate (2.18) with respect to 
$x$, multiply by the known function 
$u_x^{\prime }$, and take an area integral to construct the term 
$\langle u_x^{\prime }({\partial c^{\prime }}/{\partial x})\rangle$. We perform similar procedures for each term in (2.11), and we arrive at
$$\begin{gather} \frac{\partial\langle c\rangle}{\partial t} + U(x)\left(1+\frac{\beta^2}{12}\right)\frac{\partial\langle c\rangle}{\partial x}={-}D\left(\frac{2\beta^3}{a}+\frac{3}{2}\,\beta\gamma-\frac{2\beta}{a}\right)\frac{\partial\langle c\rangle}{\partial x}\nonumber\\ {}+D\left(1-\frac{1}{12}\,\frac{U(x)\,a\beta}{D}+\frac{U(x)^2\,a^2}{48D^2}-\beta^2\right)\frac{\partial^2\langle c\rangle}{\partial x^2}. \end{gather}$$We rank the terms in (2.19) and summarize them in table 3. The detailed derivation of (2.19) is summarized in § A of the supplementary material available at https://doi.org/10.1017/jfm.2023.504.
Table 3. Ranked summary of terms in (2.19).
Keeping terms of the order of 
$({Uc_0}/{\sigma })\,{O}(\eta \,Pe_{a_0})$ and 
$({Dc_0}/{\sigma ^2})\,{O}({\sigma }/{\lambda })$ yields the following expression for the development of the area-averaged concentration:
\begin{equation} \frac{\partial\langle c\rangle}{\partial t}+\left(U(x)-\frac{2\beta D}{a}\right)\frac{\partial\langle c\rangle}{\partial x} = D\left(1+\frac{U(x)^2\,a^2}{48D^2}\right)\frac{\partial^2\langle c\rangle}{\partial x^2}, \end{equation}which can also be rearranged and written in the following conservative form (H.A. Stone, personal communication 14 December 2022):
\begin{equation} \frac{\partial\langle c\rangle}{\partial t}+U(x)\,\frac{\partial\langle c\rangle}{\partial x} = \frac{1}{a^2}\,\frac{\partial}{\partial x}\left(\left(D+\frac{U(x)^2\,a^2}{48D}\right)a^2\,\frac{\partial\langle c\rangle}{\partial x}\right). \end{equation}
Note that for 
$\beta =0$, (2.20) reduces to the simple result of Taylor–Aris dispersion for laminar, fully-developed flow within a cylindrical channel (with uniform radius) (Aris Reference Aris1956).
In the next subsection, we will use (2.20) to develop a description of the development of the axial mean position of the solute and its axial variance. For now, we note the right-hand side of (2.20) contains a prefactor for the second axial derivative 
${\partial ^2\langle c\rangle }/{\partial x^2}$, which has the form of a typical Taylor–Aris dispersion coefficient:
\begin{equation} D_{{eff}}(x)\equiv D\left(1+\frac{1}{48}\,\frac{U(x)^2\,a^2}{D^2}\right). \end{equation}
In a cylindrical tube with constant radius, the axial variance grows linearly in time according to 
$2D_{{eff}}t$, as 
$\beta =\gamma =0$ and the mean velocity correction due to 
$u_r\neq 0$ vanishes. However, we note that the coefficient 
$D_{{eff}}$ is by itself not useful in predicting the time evolution of axial variance when the cross-section is varying in space. This is true even for the simple case of a linearly converging or diverging channel (i.e. 
$\beta =\mathrm {Constant}$). The reason for this is that development of the solute zone variance is a strong function of the advective operator on the left-hand side. These advective gradients can stretch or shrink the variance as the solute navigates through contractions and expansions, respectively.
2.2. Time evolution of mean and variance
We next formulate the problem in terms of analytical expressions for two moments of solute distributions. Note that despite this use of spatial moments, our analysis does not follow Aris’ famous method of moments (Aris Reference Aris1956) as here we focus on the mean and variance of the solute distribution in the typical Taylor–Aris limit where observation times are much longer than the radial diffusion time. Our analysis also uses scaling analyses and approximations for key terms that are not part of the method of moments.
We first define an axial average as
\begin{equation} \overline{({\cdot})}\equiv\int_{-\infty}^{\infty}({\cdot})\,a^2(x^{\prime})\,\langle c\rangle(x^{\prime}, t)\,\mathrm{d} x^{\prime}\Big/\int_{-\infty}^{\infty}a^2(x^{\prime})\,\langle c\rangle(x^{\prime}, t)\,\mathrm{d} x^{\prime}, \end{equation}
where 
$x^{\prime }$ is a dummy variable for integration. For more compact notation, we also define 
$\bar {c}$ as the axial integration of the mean concentration 
$\int _{-\infty }^{\infty }a^2(x^{\prime })\,\langle c\rangle (x^{\prime }, t)\,\mathrm {d} x^{\prime }\equiv \bar {c}=\mathrm {Constant}$.
To derive the first moment, we multiply (2.20) by 
$x$ and apply the axial average operation 
$\overline {({\cdot })}$:
$$\begin{gather} \int_{-\infty}^{\infty}xa^2\,\frac{\partial\langle c\rangle}{\partial t}\,\mathrm{d} x+\int_{-\infty}^{\infty}xUa^2\,\frac{\partial\langle c\rangle}{\partial x}\,\mathrm{d} x - \int_{-\infty}^{\infty}xa^2\,\frac{2\beta D}{a}\,\frac{\partial\langle c\rangle}{\partial x}\,\mathrm{d} x\nonumber\\ =\int_{-\infty}^{\infty}xa^2D\left(1+\frac{1}{48}\,\frac{U^2a^2}{D^2}\right)\frac{\partial^2\langle c\rangle}{\partial x^2}\,\mathrm{d} x. \end{gather}$$
We evaluate the integrals using integration by parts, divide both sides by 
$\bar {c}$, and arrive at an ODE for the axial mean position 
$\bar {x}(t)$ as follows:
\begin{equation} \frac{\mathrm{d}\bar{x}}{\mathrm{d}t} =U_0a_0^2\,\overline{\left[\frac{1}{a^2}\right]}+2D\,\overline{\left[\frac{\beta}{a}\right]}. \end{equation}
The two terms on the right-hand side are expressions of the form 
$\overline {f(x)}$. Some analysis of this term is essential for further analytical simplification of the problem. We consider a Taylor expansion of 
$f(x)$ about the mean axial position 
$\bar {x}$ as follows:
\begin{equation} f(x)=f(\bar{x})+(x-\bar{x})\,f'(\bar{x})+\tfrac{1}{2}\,(x-\bar{x})^2\,f''(\bar{x})+\tfrac{1}{6}\,(x-\bar{x})^3\,f'''(\bar{x})+\cdots. \end{equation}
To simplify the expressions, we define 
$\tilde {f}\equiv f(\bar {x})$. Inserting the expression in (2.26) into (2.25) and keeping the dominant term, we derive the following ODE for 
$\bar {x}(t)$:
\begin{equation} \frac{\mathrm{d}\bar{x}}{\mathrm{d}t}= U(\bar{x})+2D\,\widetilde{\left[\frac{\beta}{a}\right]}+{O}\left(U_0\left(\frac{\sigma}{\lambda}\right)^2, \frac{D}{\lambda}\left(\frac{\sigma}{\lambda}\right)^2\right), \end{equation}
which reduces to the standard Taylor result for 
$\beta (x)=0$.
To derive the second moment equation of (2.20), we multiply (2.20) by 
$(x-\bar {x})^2$, apply the axial average operation, perform integration by parts, and divide both sides by 
$\bar {c}$. This then yields the following ODE for the dynamics of 
$\sigma _x^2$:
\begin{equation} \frac{\mathrm{d}\sigma^2_x}{\mathrm{d}t}= 2Ua^2\left[\overline{\frac{x}{a^2}}-\bar{x}\,\overline{\frac{1}{a^2}}\right] +\frac{U^2a^4}{24D}\,\overline{\frac{1}{a^2}}+2D+4D\,\overline{\frac{\beta}{a}\,(x-\bar{x})}. \end{equation}
Note that the right-hand side of (2.28) includes expressions of the form 
$\overline {f(x)}$, 
$\overline {(x-\bar {x})\,f(x)}$, and 
$\overline {x\,f(x)}$. We again simplify the latter second and third terms using Taylor expansion at the mean axial position 
$\bar {x}$ as follows:
\begin{align} \left. \begin{aligned}
(x-\bar{x})\,f(x) &
=(x-\bar{x})\,f(\bar{x})+(x-\bar{x})^2\,f'(\bar{x})+\tfrac{1}{2}\,(x-\bar{x})^3\,f''(\bar{x})+\tfrac{1}{6}\,(x-\bar{x})^4\,f'''(\bar{x})\cdots,\\
x\,f(x) & =\bar{x}\,f(\bar{x}) +
(x-\bar{x})(\bar{x}\,f^{\prime}(\bar{x})+f(\bar{x}))+(x-\bar{x})^2\big(\tfrac{1}{2}\,\bar{x}\,f^{\prime\prime}(\bar{x})+f^{\prime}(\bar{x})\big)\\
& \quad +\tfrac{1}{6}\,(x-\bar{x})^3(\bar{x}\,f'''(\bar{x})+3f''(\bar{x}))+\cdots.
\end{aligned} \right\} \end{align}
We insert the expression in (2.29) into (2.28), keep the terms of order 
$D\times {O}(1)$, 
$D\times {O}(Pe_{a_0}^2)$, 
$D\times {O}(Pe_{a_0}({\sigma ^2}/{a_0\lambda }))$, and arrive at an ODE of the axial variance 
$\sigma _x^2(t)$ as
\begin{align}
\frac{\mathrm{d}\sigma^2_x}{\mathrm{d}t}=2D+\frac{U^2a^4}{24D}\,\widetilde{\left[\frac{1}{a^2}\right]}
-
4Ua^2\sigma^2_x\,\widetilde{\left[\frac{\beta}{a^3}\right]}+D\,{O}\Bigl(\left(\frac{\sigma}{\lambda}\right)^2,
Pe_{a_0}^2\left(\frac{\sigma}{\lambda}\right)^2,
Pe_{a_0}\,\frac{\sigma^2}{a_0\lambda}\,\frac{\sigma}{\lambda}\,\mathrm{Skew}_x
\Bigr).
\end{align}
In this equation, as 
$\sigma /\lambda$ approaches unity, the dominant error term is the 
$D\times {O}(Pe_{a_0}({\sigma ^2}/{a_0\lambda })({\sigma }/{\lambda })\,\mathrm {Skew}_x )$ term. For 
$\sigma /\lambda$ approaching unity, this has the order of magnitude
\begin{equation} 2U_0a_0^2\,\mathrm{Skew}_x\sigma_x^3\,\widetilde{\frac{3\beta^2-a\gamma}{a^4}}, \end{equation}
Here, 
$\mathrm {Skew}_x$ is the axial skewness of the solute zone, defined as 
$\mathrm {Skew}_x=\overline {(x-\bar {x})^3}/\sigma _x^3$. Detailed derivation of (2.31) is summarized in § B of the supplementary material. Note that (2.30) also reduces to standard Taylor results when 
$\beta (x)=0$. Equations (2.27) and (2.30) are two coupled ODEs that predict the time evolution of the axial mean and variance of the solute based on channel geometry and the assumption of lubrication theory in the velocity field. In § 3, we will use these two equations to predict the time evolution of the axial mean and variance of the solute zone for various interesting channel geometries.
2.3. Analytical derivation of the boundary between the regimes of transient positive and negative growth of axial variance
Advective dispersion can cause the axial variance of the solute zone to decrease as the solute travels through a region where the channel is expanding. Such an expanding channel region can be characterized qualitatively by a positive and sufficiently large value of 
$\beta$ and a sufficiently large 
$Pe_a$ (so that molecular diffusion does not completely overpower the advective effect), and a sufficiently large value of the axial variance relative to the local channel radius. We explore this effect here. Enabled by our analytical approach, we will formulate the boundary between the physical regimes of transient positive and negative growth of axial variance of the solute. Rearranging terms in (2.30), we see that
\begin{equation} \frac{1}{D}\,\frac{\mathrm{d}\sigma^2_x}{\mathrm{d}t}=2\left(1+\frac{1}{48}\,Pe_a^2\right)-4\,Pe_a\,\beta\left(\frac{\sigma_x^2}{a^2}\right). \end{equation}
This expression yields a useful description of the boundary between the regimes of transient positive and transient negative growth of axial variance in terms of 
$\beta$, 
$Pe_a$ and 
$\sigma _x^2/a^2$. Note that here the negative growth of variance is a transient phenomenon, and is limited to simply the axial width of the sample. The growth of the three-dimensional extent of the solute cloud can, of course, never be negative.
Figure 2 shows a plot of this boundary as a three-dimensional surface in terms of the aforementioned three variables. That is, the surface (figure 2a) delineates the regions of transient positive and negative growth of axial variance as a function of local Péclet number (
$Pe_a$), local slope of channel radius distribution (
$\beta$), and local ratio between variance and squared radius (
$\sigma _x^2/a^2$). Figure 2(b) shows the contours of horizontal cross-section of the surface for various values of 
$\ln (\sigma _x^2/a^2)$ (labelled on each contour line). Points above the surface exhibit a transient negative rate of axial variance growth, while points below have a positive rate of growth. Accordingly, the surface defines the solution for transient zero growth in variance in this three-parameter space. As expected, the surface asymptotes towards the 
$Pe_a$–
$\beta$ plane for finite 
$Pe_a$ and larger values of 
$\beta$. We also see that for fixed and finite values of (positive) 
$\beta$, setting the left-hand side of (2.32) to zero reveals the critical variance value 
$\sigma _x^2/a^2$ required for transient negative growth of axial variance:
\begin{equation} \frac{\sigma_x^2}{a^2}=\frac{1}{4\beta}\left(\frac{2}{Pe_a}+\frac{Pe_a}{24}\right). \end{equation}
We now differentiate the latter expression with respect to 
$Pe_a$, set it equal to zero, and arrive at the minimal variance required for transient negative axial variance growth:
\begin{equation} \min\left(\frac{\sigma_x^2}{a^2}\right)=\frac{1}{4\sqrt{3}\,\beta}. \end{equation}
The minimum variance occurs at 
$Pe_a=\sqrt {48}$.
Figure 1. Schematic of an axisymmetric channel with a slowly varying, arbitrary distribution of radius 
$a(x)$. A nominal radius is taken as 
$a_0$ at 
$x=0$. The slope and the curvature of the cylinder wall are respectively 
$\beta (x)$ and 
$\gamma (x)$, as shown. Here, 
$\sigma _{x}$ is the characteristic width of a solute zone, and 
$U(x)$ is the area-averaged axial velocity distribution. The (
$x$-direction) length of the channel depicted schematically in the sketch has been compressed relative to its characteristic radius for clarity of presentation.
Figure 2. Regimes of transient positive and negative growth of axial variance of the solute. (a) Surface of zero variance growth in a space of Péclet number, 
$\beta$, and the natural log of the local ratio between variance and squared radius (
$\ln (\sigma _x^2/a^2)$). This surface is computed analytically from (2.32). Shown is a surface where the axial variance rate of growth 
${\mathrm {d}\sigma _x^2}/{\mathrm {d}t}$ is approximately 0. (b) Contour plot showing the horizontal cross-section curves of the zero transient axial variance growth surface in a space of Péclet number, 
$\beta$ and 
$\ln (\sigma _x^2/a^2)$, at different 
$\ln (\sigma _x^2/a^2)$ values (labelled on each contour line).
2.4. Approximations useful for channel design
The analytical expression of (2.32) is useful in engineering channel geometries which will generate desired variance evolutions in the Taylor–Aris dispersion regime. We can simplify such a solution further by implementing the approximation
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\approx\frac{\mathrm{d}\bar{x}}{\mathrm{d}t}\,\frac{\mathrm{d}}{\mathrm{d} x}. \end{equation}
Inserting this approximation into (2.30), and substituting the definition of 
$\beta$, we obtain an ODE for 
$a(x)$ given a desired spatial distribution 
$\sigma _x^2(x)$:
\begin{equation} \frac{\mathrm{d}a}{\mathrm{d} x}\approx\dfrac{-\dfrac{\mathrm{d}\sigma_x^2}{\mathrm{d} x}+\dfrac{2a^2D}{U_0a_0^2}+\dfrac{U_0a_0^2}{24D}}{\dfrac{4\sigma_x^2}{a}+\left(\dfrac{2D}{U_0a_0^2}\,\dfrac{\mathrm{d}\sigma_x^2}{\mathrm{d} x}\right)a}. \end{equation}We can then solve (2.36) numerically to obtain the shape of the engineered channel. In § 3.4, We will use this approximation to design channels that can result in a specific variance evolution pattern in space.
2.5. Brownian dynamics simulations
We used Brownian dynamics simulations to benchmark and evaluate our analytical expressions for variance evolution. Each Brownian dynamics simulation consisted of tracking 5000 point-particles that were initially distributed uniformly along the radius and normally along the streamwise direction. We considered mean axial position zero (by definition) and initial axial variance 
$\sigma _{x,0}^2$. We set 
$\sigma _{x,0}^2 = 300a_0^2$ for the two engineered channels in § 3.4, and we set 
$\sigma _{x,0}^2 = 10a_0^2$ for all other cases. The velocity field follows (2.2). The system evolves according to the forward Euler method.
For the time steps, we set the ratio between diffusion time scale (
$a_0^2/D$) and time step of discretization to 40. Each reflection at the wall was checked twice at each time step. Unless specified, the (constant, initial) Péclet number based on initial radius (
$Pe_{a_0}$) was 10. For the three periodic channels shown in § 3.2, the initial Péclet number was 0.1, 1, 10, 100 and 1000. Five simulations from different random seed initial conditions were computed for each case, and the reported values are their average. The numerical computation of axial averaged quantities 
$\overline {f(x)}$ is computed with 
$\overline {f(x)}=\sum _{i=1}^{N_{{particle}}}f(x_i)/N_{{particle}}$, where 
$x_i$ is the axial location of the 
$i$th particle. The program was written in Python 3 and is available on Github (jrchang612/Taylor_dispersion_arbitrary). The detailed parameters used for the simulation presented in the figures are summarized in § C of the supplementary material.
3. Results and discussion
3.1. Diverging and converging conical channels
We first apply our analysis to a simple case of diverging and converging conical section channels. Figures 3(a,b) show solutions for dispersion of a solute zone in linearly diverging and converging channels, respectively. Figures 3(ai,bi) show particle zones at five non-dimensional times 
$2Dt/a_0^2$ (from 5 to 800) as they migrate through the channels. The axes are non-dimensionalized coordinates 
$r^*=r/a_0$ and 
$x^*=x/a_0$. Note the intensely exaggerated height-to-length aspect ratio of the figures, chosen here for clarity of presentation. The top halves of the solute concentration in the channels are raw data results of the three-dimensional Brownian dynamics simulations. The bottom halves of the channels show the (one-dimensional) area-averaged axial distribution of concentration of solute from the model developed here (i.e. numerical solutions of ODEs given by (2.27) and (2.30)). Note that we present and plot the predicted axial distribution assuming a Gaussian distribution as we have only the information about mean and variance. However, the derivations of (2.27) and (2.30) do not rely on a Gaussian distribution assumption. The solute is initially uniformly distributed over the cross-sectional area of the channel, and the initial standard deviation widths of the axial solute distributions are each set equal to the same value (equal to 
$\sqrt {10}\,a_0=\sqrt {10}\,a(x=0)$ of the diverging channel). The current model has very good comparison with the Brownian dynamics simulations. The characteristic 
$Pe_a$ values here are 
${O}(10)$, and this results in negligible radial gradients of particle density in the Brownian dynamics. Hence the particle clouds from Brownian dynamics are fairly symmetrical about 
$x^*$. Note the general trend of relatively rapid increase of zone variance in the converging case. By comparison, the expansion in the diverging channel limits the growth of variance. For this case, the divergence is not sufficiently pronounced to result in transient negative growth of axial variance (but we will explore such cases below.)
Figure 3. Taylor–Aris dispersion in (a) diverging and (b) converging conical channels. (ai,bi) Results from Brownian dynamic simulation (upper half) along with the predicted axial distribution of area-averaged concentration (lower half). (aii,bii) Solute axial variance as a function of 
$\bar {x}^*$. Plots show axial variance from the current analytical model (subscript 
$curr$, (2.30)) and from Brownian dynamics simulation (subscript 
$B$). Axial variance computed using (2.30) shows excellent agreement with Brownian dynamics simulations.
Figures 3(aii,bii) are plots of the axial variance of the area-averaged axial distribution of solute concentration from the current model (
$\sigma _{x,{curr}}^{*2}$) and Brownian dynamics (
$\sigma _{x,{B}}^{*2}$). These quantities are plotted as a function of the axial mean position of the solute zone, 
$\bar {x}^*$. There is excellent agreement between the current model and the Brownian dynamics simulation for both the diverging and converging cases, demonstrating the ability of the current model to predict the spatial evolution pattern of axial variance.
3.2. Periodic channels
We next apply our model to a periodic channel with sinusoidal radius distribution. Similar to figures 3(ai,bi), figure 4(a) shows a plot of the channel geometry in coordinates 
$r^*$ versus 
$x^*$ in a highly exaggerated aspect ratio for clarity. We show results for 
$Pe_{a_0}=10$, where 
$a_0$ is defined as the radius at 
$x=0$ and is also the axially averaged radius. Again, the channel shows plots of solute zones at non-dimensional times ranging from 5 to 800. The top halves of the solute zones are raw results from the Brownian dynamics simulations, while the bottom halves are plots of the predicted axial distribution of area-averaged concentration from the current model. Again, we see agreement between the axial distributions of Brownian dynamics particle densities and the axial distributions from the model.Note how the model captures the strongly positive and negative growth of axial variance as the solute traverses through contractions and expansions, respectively. For example, note the rapid increases in axial variance due to the contraction just upstream of 
$x^*=1800$, and then the subsequent rapid decrease in axial variance caused by the advective effects of the solute expanding from 
$x^*=1800$ to just upstream of 
$x^*=3000$. Note that the particle clouds from Brownian dynamics are also fairly symmetric about 
$x^*$.
Figure 4. Taylor–Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, 
$a(x)$. (a) Results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (2.30). (b) Distribution of 
$\sigma _{x,{B}}^{*2}$, 
$\sigma _{x,{curr}}^{*2}$, 
$U^*(x)$ and 
$D_{{eff}}/D$, each as a function of the non-dimensional axial location along the channel, 
$\bar {x}/a_0$. Here, 
$U^*(x)$ and 
$D_{{eff}}^*$ are periodic functions, shown as a reference.(c) Predicted (2.31) and observed errors of our current model in the growth rate of axial variance. Also shown for reference is the ratio between the square root of the axial variance and channel wavelength 
$\sigma _x/\lambda$. The variance computed using (2.30) shows excellent agreement with Brownian dynamics simulations. Note that variance averaged along the axial spatial period increases monotonically as expected. The error in axial variance growth rate is small and matches the predicted error from (2.31). This shows the accuracy of our model when 
$\sigma _x/\lambda$ is small.
Figure 4(b) shows the axial variance of the solute zone as a function of solute average axial location 
$\bar {x}^*$. The axial variances from the Brownian dynamics simulation and the current model are plotted. Note the excellent agreement between these two models, showing how the current model captures very well the detailed development of the axial variance. Note also how the axial variance at equal values of the phase increases, and the axial variance averaged over the period increases monotonically. At the end of the current section, we will consider this development further and use our model to analyse period-averaged axial variance over long times.
Figure 4(b) also shows the non-dimensional area-averaged flow velocity and the local normalized effective dispersion coefficient, 
$D_{{eff}}/D$ (cf. two scales on the right-hand side of the plot). Both 
$U$ and 
$D_{{eff}}^*$ are periodic functions, antiphase to the shape of the channel; 
$D_{{eff}}^*$ increases as the channel contracts, and decreases as the channel expands, as we expected. However, even when the variance is actually decreasing (e.g. between 
$x^* = 2000$ and 
$x^*=3000$), the effective dispersion coefficient 
$D_{{eff}}/D$ remains positive and greater than unity. This demonstrates the inability of the effective dispersion coefficient alone to describe the axial variance evolution (including here in a periodic channel).
Figure 4(c) shows the observed and predicted errors in dimensionless growth rate of axial variance of the solute as a function of 
$\bar {x}^*$. Also shown is the ratio between the square root of axial variance and channel wavelength 
$\sigma _x/\lambda$. Consistent with figure 4(b), the observed error in axial variance growth rate is small and stochastic, and it matches with the predicted error from (2.31). The ratio 
$\sigma _x/\lambda$ is much smaller than 1 throughout the simulation, confirming the accuracy of our model when the assumptions listed in table 1 are satisfied.
We next increase the Péclet number of the system to evaluate how the model assumptions become inaccurate and the current model fails. We repeat our simulation in the same channel as presented in figure 4, but increase the initial Péclet number to 100. Similar to figure 4(a), figure 5(a) shows a plot of the channel geometry in coordinates of 
$r^*$ versus 
$x^*$ in a highly exaggerated aspect ratio for clarity. Again, the channel shows plots of solute zones at non-dimensional times ranging from 5 to 700, with non-dimensional time points at 318 and 492 to demonstrate the best- and worst-case scenarios in predictions. The top halves of the solute zones are raw results from the Brownian dynamics simulations, while the bottom halves are plots of the predicted axial distribution of area-averaged concentration from the current model. Again, we see overall agreement between the axial distributions of Brownian dynamics particle densities and the axial distributions from the model. The model captures the strongly positive and transient negative growth of axial variance as the solute traverses contractions and expansions, respectively. Our model tends to overestimate the fluctuation in axial variance, especially when the solute traverses contractions. This overestimation is likely due to the fact that the width of the solute zone becomes comparable to the wavelength of the channel. For example, when the solute mean position is within a contraction, the leading and trailing ends of the solute zone are in expansion zones.
Figure 5. Taylor–Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, 
$a(x)$. The channel is the same as that in figure 4, but the initial Péclet number 
$Pe_{a_{0}}$ is 100. (a) Results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (2.30). (b) Distributions of 
$\sigma _{x,{B}}^{*2}$ and 
$\sigma _{x,{curr}}^{*2}$, each as a function of the non-dimensional axial location along the channel, 
$\bar {x}/a_0$. (c) Predicted (2.31) and observed error of our current model in the growth rate of axial variance. Also shown for reference is the ratio between the square root of axial variance and channel wavelength 
$\sigma _x/\lambda$. Variance computed using (2.30) shows excellent agreement with Brownian dynamics simulations when the channel is expanding, but overestimates the axial variance when the channel is contracting. Note that variance averaged along the axial spatial period increases monotonically as expected. The predicted error in axial variance growth rate using (2.31) shows excellent agreement with the observed error, especially when 
$\sigma _x/\lambda$ is small. As 
$\sigma _x/\lambda$ increases, there is more deviation between observed and predicted error in axial variance growth rate.
Figure 5(b) shows the axial variance of the solute zone as a function of 
$\bar {x}^*$. The variances from the Brownian dynamics simulation and the current model are plotted. There is good agreement between these two models in trends of development of the variance, and excellent agreement when 
$\bar {x}^*$ is less than 5000. The agreement between the two models remains very good when the solute passes through expansions in the channel, but the current model overestimates the axial variance when the solute passes through a contraction.
Figure 5(c) again shows the observed and predicted error in dimensionless growth rate of axial variance of solute as a function of 
$\bar {x}^*$. Also shown is the ratio between the square root of axial variance and channel wavelength 
$\sigma _x/\lambda$. The magnitude of observed error grows as the axial variance grows, and the error is most pronounced when the solute traverses contractions. There is excellent agreement between the predicted error from (2.31) when 
$\bar {x}^*$ is less than 20 000, but there is more deviation between the two predictions as 
$\sigma _x/\lambda$ grows larger.
We further increase the initial Péclet number of the system to 1000 using the same channel of figures 4 and 5. Similar to figure 5(a), figure 6(a) shows a plot of the channel geometry in coordinates of 
$r^*$ versus 
$x^*$ in a highly exaggerated aspect ratio for clarity. Again, solute zones at non-dimensional times ranging from 5 to 80 are shown. The top halves of the solute zones are raw results from the Brownian dynamics simulations, while the bottom halves are Gaussian distributions with mean and axial variances predicted by the current model. The agreement between the two is now merely qualitative. Note that for these conditions, the sample zone axial width quickly becomes of the same order as the wavelength of the channel.
Figure 6. Taylor–Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, 
$a(x)$. The channel is the same as that in figure 4, but the initial Péclet number 
$Pe_{a_{0}}$ is 1000. (a) Results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (2.30). (b) Distributions of 
$\sigma _{x,{B}}^{*2}$ and 
$\sigma _{x,{curr}}^{*2}$, each as a function of the non-dimensional axial location along the channel, 
$\bar {x}/a_0$. (c) Predicted (2.31) and observed error of our current model in the growth rate of axial variance. Also shown for reference is the ratio between the square root of axial variance and channel wavelength 
$\sigma _x/\lambda$. Variance computed using (2.30) shows agreement with Brownian dynamics simulations only in the very beginning. Our model fails to predict the monotonic growth of axial variance when the axial variance is of the same order as the channel wavelength. The predicted error in axial variance growth rate using (2.31) shows fair agreement with the observed error for 
$\bar {x}^{*}$ below 10 000. As the ratio between axial variance and channel wavelength approaches unity, the predicted error also become inaccurate as the basic assumptions of our model are no longer valid.
Figure 6(b) shows the axial variance of the solute zone as a function of 
$\bar {x}^*$. The variances from the Brownian dynamics simulation and the current model are plotted. The agreement between the two is limited to the initial development of the axial variance (
$\bar {x}^*<2000$). The axial variance growth in Brownian dynamics simulation becomes monotonic after 
$\bar {x}^*>30\,000$, but the current model is not able to capture this.
Figure 6(c) again shows the observed and predicted error in dimensionless growth rate of axial variance of the solute as a function of 
$\bar {x}^*$. Also shown is the ratio between the square root of axial variance and channel wavelength 
$\sigma _x/\lambda$. The magnitude of observed error grows as the axial variance grows, and the error is now pronounced in both channel contractions and expansions. The predicted error from (2.31) can capture only the observed error in the initial development of axial variance (
$\bar {x}^*<10\,000$), but becomes inaccurate as 
$\sigma _x/\lambda$ grows above 0.3. This finding suggests that the accuracy of our model depends less on the Péclet number itself but is more sensitive to the ratio 
$\sigma _x/\lambda$. At large Péclet number, since the axial variance growth rate is much faster (as discussed in figure 7) and 
$\sigma _x/\lambda$ quickly approaches and exceeds unity, the current model becomes inaccurate. A zoomed-in comparison of the solute zone subject to Péclet numbers 10, 100 and 1000 is in figure S1 of the supplementary material.
Figure 7. Normalized long-term effective dispersion coefficient 
$D_{{eff}}^{\infty *}$ as a function of 
$Pe_{a_0}$ for three different periodic channels with equal period and similar radius amplitude (
$\delta =0.05$, 
$\lambda =200$). All plots show 
$D_{{eff}}^{\infty *}$ computed using (2.30) and with a Brownian dynamics simulation. Plot (d) also shows a comparison with the expression derived by Adrover et al. (Reference Adrover, Venditti and Giona2019) for a sinusoidal channel. Note that the plots across all three channels are very similar (but not exactly the same) in magnitude and shape. This similarity reflects the fact that the long-term development of the solute in periodic channels is most strongly a function of channel amplitude and a weak function of channel shape.
We next consider the long-term averaged dispersion coefficient in periodic channels. As with previous studies of dispersion (Hoagland & Prud'Homme Reference Hoagland and Prud'Homme1985; Adrover et al. Reference Adrover, Venditti and Giona2019), we will define a new effective dispersion coefficient for the long-term (many-period) growth of the axial variance 
$D_{{eff}}^{\infty }$ defined as
\begin{equation} D_{{eff}}^{\infty}\equiv\lim_{t\rightarrow\infty}\frac{\sigma_x^2}{2t}. \end{equation}
The non-dimensional version of this quantity will be defined as 
$D_{{eff}}^{\infty *}\equiv D_{{eff}}^{\infty }/D$. We analysed long-term dispersion behaviour in periodic channels with three different shapes but with the same period and similar radius amplitude. Informed by our analysis in the previous section, however, the predicted variance growth becomes inaccurate when the ratio 
$\sigma _x/\lambda$ is greater than about 0.3. We thus computed the long-term growth of the axial variance from our current model as
\begin{equation} D_{{eff, curr}}^{\infty}=\frac{1}{\left.t\right|_{\sigma_x=0.3\lambda}}\int_{0}^{\left.t\right|_{\sigma_x=0.3\lambda}}\frac{\sigma_x^2(t)}{2t}\,\mathrm{d}t. \end{equation} Figure 7 summarizes the results of this study. Figures 7(a–c) show plots of the radius distribution 
$a(x)$ normalized by initial radius 
$a_0$ as a function of the normalized axial coordinate 
$x^*$. The three functions 
$a(x)$ are sinusoidal, triangular and exponential sine wave functions of the following forms:
\begin{equation} \left. \begin{aligned} a_1(x) & =1+\delta\sin(2{\rm \pi} x/\lambda),\\ a_2(x) & =1+\frac{3\delta}{2\lambda}\,\mathrm{mod}(x,\lambda)-\frac{9\delta}{2\lambda}\left(\mathrm{mod}(x,\lambda)-\frac{2}{3}\lambda\right)H\left(\mathrm{mod}(x,\lambda)-\frac{2}{3}\lambda\right),\\ a_3(x) & =\frac{\delta}{e}\exp(\sin(2{\rm \pi} x/\lambda)) + \left(1-\frac{\delta}{e}\right), \end{aligned} \right\} \end{equation}
where 
$H(x)$ is the Heaviside step function. Figures 7(d–f) show the effective, long-term dispersion coefficient 
$D_{{eff}}^{\infty *}$ for each case as a function of 
$Pe_{a_0}$. We show 
$D_{{eff}}^{\infty *}$ curves for the current model (data points) and Brownian simulations (solid line). For these conditions, we see excellent agreement between the current model and the Brownian simulations for the long-term dispersion coefficient. The results capture an asymptote of 
$D_{{eff}}^{\infty *}=1$ for vanishing 
$Pe_{a_0}$, as expected. Here, 
$D_{{eff}}^{\infty *}$ increases monotonically with increasing 
$Pe_{a_0}$, and asymptotes to a 
$Pe_{a_0}^2$ dependence for large 
$Pe_{a_0}$.
In figure 7(d), we also show a comparison of our model results with the work of Adrover et al. (Reference Adrover, Venditti and Giona2019). Adrover analysed the long-term dispersion of solutes in a sinusoidal channel and provided an approximate analytical formula for 
$D_{{eff}}^{\infty *}$ as a function of 
$Pe_{a_0}$, and this prediction is plotted along with our model in the figure. Adrover also found that 
$D_{{eff}}^{\infty *}$ tends to unity for small 
$Pe_{a_0}$, and scales with the second power of 
$Pe_{a_0}$ for large 
$Pe_{a_0}$. We note the excellent agreement among the three predictions.
We conclude the current model can be adapted readily to a wide variety of periodic channel shapes. We further note the similarity among the three 
$D_{{eff}}^{\infty *}$ versus 
$Pe_{a_0}$ curves for the three cases. This similarity leads us to hypothesize that the long-term solute dispersion in such periodic channels is driven largely by the spatial frequency and the amplitude of the radius oscillation (and 
$Pe_{a_0}$), and may be insensitive to the details of channel shapes.
3.3. Arbitrarily shaped channels
We next demonstrate a novel application of our model to a particular but arbitrary axisymmetric channel shape. Similar to figure 4(a), figure 8(a) shows a plot of the channel geometry in coordinates of 
$r^*$ versus 
$x^*$ in a highly exaggerated aspect ratio for clarity. We show results for 
$Pe_{a_0}=10$ where 
$a_0$ is the radius at 
$x=0$. The channel shows plots of solute zones at non-dimensional times ranging from 25 to 250. The top halves of the solute zones are raw results from the Brownian dynamics simulations, while the bottom halves are plots of the predicted axial distribution of area-averaged concentration from the current model. Again we see excellent agreement between the axial distributions of Brownian dynamics particle densities and the axial distributions from the model. Note how the model captures the sudden positive and negative growth of axial variance as the solute traverses contractions and expansions, respectively. For example, note the sudden decreases in variance due to the expansion just upstream of 
$x^*=500$, and then the sudden increases in variance caused by the contraction near 
$x^*=700$.
Figure 8. Taylor–Aris dispersion for an example arbitrarily shaped axisymmetric channel. (a) Results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (2.30). (b) Distributions of 
$\sigma _{x,{B}}^{*2}$, 
$\sigma _{x,{curr}}^{*2}$, 
$U^*(x)$ and 
$D_{{eff}}/D$, each as a function of the non-dimensional axial location along the channel. Here, 
$U^*(x)$ and 
$D_{{eff}}/D$ are shown for reference. Variance computed using (2.30) shows excellent agreement with Brownian dynamics simulations. (c) The same example used for demonstration and benchmarking of (2.32), plotting a scaled rate of change of axial variance as a function of mean solute position 
$\bar {x}$ as computed using (2.30) and analytical expression (2.32). Both of these solutions are compared to calculations based on a Brownian dynamics simulation. The blue curve is the Brownian smoothed with a moving average, with window size 
$\Delta t^*=2.5$, while the light blue curve shows the original Brownian simulation results.
Figure 8(b) shows the axial variance of the solute zone as a function of non-dimensional solute average axial location, 
$\bar {x}^*$. The axial variance from the Brownian dynamics simulation and the current model are plotted. Note the excellent agreement between these two models, again showing how the current model captures very well the detailed development of the axial variance. For example, note the sudden decreases of axial variance just upstream of 
$\bar {x}^*=500$ and 
$\bar {x}^*=800$, and the sudden increases of axial variance just upstream of 
$\bar {x}^*=700$.
Figure 8(b) also shows the non-dimensional area-averaged flow velocity and the local normalized effective dispersion coefficient, 
$D_{{eff}}/D$ (cf. two scales on the right-hand side). Similar to the case in figure 4(b), both 
$U$ and 
$D_{{eff}}^*$ show trends opposite to the shape of the channel. Again, 
$D_{{eff}}^*$ increases as the channel constricts and decreases as the channel expands, as we expected. However, even in a region where axial variance is decreasing transiently (e.g. between 
$x^*=400$ and 
$x^*=500$), the effective dispersion coefficient remains positive and greater than unity. This again highlights the inability of the effective dispersion coefficient to describe the variance evolution, here in an arbitrarily shaped axisymmetric channel.
Figure 8(c) shows the scaled rate of change of variance (
$({1}/{D})({\mathrm {d}\sigma _x^2}/{\mathrm {d}t})$) as a function of solute average axial location 
$\bar {x}^*$. The results computed from Brownian simulations and the current model (2.30) are plotted. For the Brownian dynamics simulation, both raw data and the moving averaged results are shown. The axial variance growth rate varies according to the geometry of the channel (shown in figure 8(a), decreasing when the channel expands and increasing when the channel contracts, as we expect. There is also excellent agreement among the three solutions, and the current model captured the negative variance growth rate near 
$\bar {x}^*=400$ and 
$\bar {x}^*=700$. This shows the capability of our model to predict the axial variance evolution, and also provides validation for using (2.32) to identify the regime of transient negative axial variance growth.
3.4. Engineering channel shape for specific variance patterns
We apply our analytical approach to design the channel shape to have a desired spatial evolution of axial variance. As a proof of concept, we design one channel that maintains an approximately constant axial variance (channel A) and a second channel that results in a sinusoidal variation of axial variance as the solute develops in the channel (channel B).
The two channel shapes are designed by solving (2.36). The shape of channel A diverges monotonically, and the rate of diverging increases downstream. For channel B, there is an increased diverging rate between 
$x^*=300$ and 
$x^*=500$, coinciding with the region where it is necessary to reduce the variance according to the targeted sinusoidal pattern. Note that the channel A shape in figure 9(a) has an analytical expression. By setting the 
${\mathrm {d}\sigma _x^2}/{\mathrm {d} x}$ term on the right-hand side of (2.36) to 0, we can simplify (2.36) into the following ODE for 
$a(x)$:
\begin{equation} \frac{\mathrm{d}a}{\mathrm{d} x}\approx \frac{D}{2\sigma_x^2U_0a_0^2}\,a^3+\frac{U_0a_0^2}{96\sigma_x^2D}\,a, \end{equation}which has the analytical solution
\begin{equation} a^2(x)=\frac{\dfrac{U_0a_0^2}{96\sigma_x^2D}\exp\left({\dfrac{U_0a_0^2}{48\sigma_x^2D}\,(c_1+x)}\right)} {1-\dfrac{D}{2\sigma_x^2U_0a_0^2}\exp\left({\dfrac{U_0a_0^2}{48\sigma_x^2D}\,(c_1+x)}\right)}, \end{equation}
where 
$c_1$ is the constant of integration related to the initial variance. We know of no analytical solution for the equation for the case of channel B.
Figure 9. Engineering the axial variance evolution in Taylor–Aris dispersion. Using (2.36), we designed two channels, A and B, which (a) maintain an approximately constant axial variance, and (b) result in a sinusoidal (axial) variation of axial variance as the solute develops in the channel, respectively. (ai,bi) Results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (2.30). (aii,bii) Distributions of 
$\sigma _{x,{B}}^{*2}$, 
$\sigma _{x,{curr}}^{*2}$ and 
$\sigma _{x,{target}}^{*2}$, each as a function of the non-dimensional axial location along the channel. Axial variance computed using (2.30), and axial variance computed from Brownian dynamic simulation, both show excellent agreement with the targeted variance evolution pattern.
Figures 9(a,b) show solutions for dispersion of the solute zone in the two engineered channels A and B, respectively. Channel A in figure 9(a) is designed to maintain approximately constant axial variance (
$\sigma _x^{2*}(\bar {x}^*)=300$), while channel B in figure 9(b) is designed to yield a sinusoidal axial variation of axial variance (
$\sigma _x^{2*}(\bar {x}^*)=300+50\sin (2{\rm \pi} \bar {x}^*/600)$). Similar to figure 3(a), figures 9(ai,bi) show the two engineered channel geometries in coordinates of 
$r^*$ versus 
$x^*$. We show results for 
$Pe_{a_0}=10$ where 
$a_0$ is taken as the radius at 
$x=0$. The top halves of the solute concentration in the channel are raw data results from Brownian dynamics simulations, while the bottom halves of the channels show the area-averaged axial distribution of solute concentration predicted from the current model.
Figures 9(aii,bii) are plots of the axial variance of the area-averaged axial distribution of solute concentration from the current model (
$\sigma _{x, {curr}}^{*2}$), Brownian dynamics (
$\sigma _{x, {B}}^{*2}$) and the targeted axial variance spatial evolution pattern (
$\sigma _{x, {target}}^{*2}$). These quantities are plotted as functions of the axial mean position of the solute zone, 
$\bar {x}^*$. For channel A, there is a near-perfect agreement among the current model, Brownian dynamics simulation and targeted pattern. All three lines stay flat at variance 300, as expected. For channel B, there is a good agreement among the three quantities. Although the amplitude of the sinusoidal variation is slightly smaller compared to the targeted pattern, the Brownian dynamics simulation results show the desired sinusoidal spatial evolution pattern. To our knowledge, our study is the first to demonstrate the design of a channel shape that yields a desired spatial evolution pattern of variance. Without the analytical approach described in §§ 2.2–2.4, this process would require repetitive simulation efforts (as with shape optimization) to compute channel shapes to obtain the desired pattern. We hypothesize that the proposed engineered channel shape can be produced by additive manufacturing methods.
Also, the most common application of microfluidics is the chemical and biochemical analyses of species (Wu et al. Reference Wu, He, Chen and Lin2016). The most common method of detection of analytes is an optical technique such as fluorescence, colorimetric or UV adsorption (Wu & Gu Reference Wu and Gu2011). These optical techniques perform line-of-sight averaging of solutes in long thin channels. Hence the detection of species is typically proportional to the area average of some analyte solute zone – and this is the main characteristic of interest of the current work. The ability to tailor a channel shape so as to preserve the area-averaged concentration therefore offers the opportunity to tailor microchannel shapes that can transport species via pressure-driven flow while also preserving and/or tailoring the depth- or area-averaged signal strength.
4. Summary and conclusions
We demonstrated a Taylor–Aris dispersion analysis for axisymmetric channels with slowly varying, arbitrary radius distributions, with assumptions and smallness parameters summarized in table 1. We first derived a PDE for the development of the area-averaged concentration, including an explicit local dispersion coefficient. We then derived equations for the dynamics of the axial mean (first moment) and variance (second moment) of the solute distribution. We proposed a heuristic for relations describing the moments of local geometry experienced by the solute. Our analysis allowed us to simplify the solution for the complex dynamics of solute zones to two coupled ODEs for the mean and variance. These ODEs provide a description of solute position and variance from the channel geometry, given the assumptions of lubrication flow. To our knowledge, this is the first time a full prediction of the time evolution of axial variance is possible using only two ODEs (including for short time scales) for this type of problem. Our method can also be applied to a long time scale, as long as the axial variance remains small compared to the characteristic wavelength of channel variations.
We further derived an analytical expression that delineates the regimes of transient positive and negative axial variance growth. This expression quantifies the solution of transient zero growth axial variance as a function of the local Péclet number, the local slope in the channel, and the ratio between axial variance and the square of the local radius. This analysis demonstrates clearly the conditions required for channel expansion to yield decreases in solute axial variance. We also developed further simplifications of our model that yield a single, first-order nonlinear ODE describing the relation between the axial radius distribution and axial variance (spatial) distribution. This relation is very useful in the design of channel shapes that yield specific (desired) dynamics for solute axial variance.
We applied our model to several interesting test cases, and benchmarked its performance relative to Brownian dynamics simulations. First, we demonstrated that our model yields accurate predictions (relative to Brownian dynamics) of area-averaged solute dynamics in diverging and converging (conical) channels. Second, we applied our model to predict solute dynamics for various periodic channels, and again demonstrated excellent agreement with Brownian dynamics simulations. For the latter, we considered an initial condition and regime where channel expansions result in substantial decreases in axial variance (i.e. transient negative solute axial variance growth). We increased the Péclet number of the channel from 10 to 1000, and studied how the performance of our model fails as the key assumptions are violated. We further analysed the long-term (many-period) developments of solute in periodic channels by defining a long-term effective dispersion coefficient. We analysed three separate periodic channel shapes, and showed that this long-term behaviour is affected mostly by channel (axial) period and the magnitude of the fluctuation of the radius function. For the case of sinusoidal channels, our model demonstrated excellent comparison with a previously published model.
The third example application of our model was for an arbitrarily shaped channel. We here selected some complex channel shape that demonstrated strong advective dispersion effects. Again, our model showed excellent comparison with Brownian dynamics simulations, including the capture of the positive and transient negative growth in axial variance when the solute zone experiences constriction or expansion in the channel.
Finally, we demonstrated the power of our analytical approach by designing two channels that can control the spatial evolution of the solute so as to produce a desired spatial distribution of the solute axial variance. For the latter work, we first designed a channel that can maintain an approximately constant solute axial variance. We then demonstrated a channel geometry that produces a sinusoidal axial distribution of axial variance as the solute develops in the channel.
Overall, our analysis provides a fairly accurate (according to Brownian dynamics simulation), fast and easy-to-use model for solute dynamics in axisymmetric channels of arbitrary variance.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2023.504.
Acknowledgements
R.C. gratefully acknowledges support from the Stanford University Bio-X SIGF Fellows Program and from the Ministry of Education in Taiwan.
Declaration of interests
The authors report no conflict of interest.
Open access
