Taylor dispersion in arbitrarily shaped axisymmetric channels

Abstract Advective dispersion of solutes in long thin axisymmetric channels is important to the analysis and design of a wide range of devices, including chemical separation systems and microfluidic chips. Despite extensive analysis of Taylor dispersion in various scenarios, most studies focus on long-term dispersion behaviour and cannot capture the transient evolution of the solute zone across the spatial variations in the channel. In the current study, we analyse the Taylor–Aris dispersion for arbitrarily shaped axisymmetric channels. We derive an expression for solute dynamics in terms of two coupled ordinary differential equations, which allow prediction of the time evolution of the mean location and axial (standard deviation) width of the solute zone as a function of the channel geometry. We compare and benchmark our predictions with Brownian dynamics simulations for a variety of cases, including linearly expanding/converging channels and periodic channels. We also present an analytical description of the physical regimes of transient positive versus negative axial growth of solute width. Finally, to further demonstrate the utility of the analysis, we demonstrate a method to engineer channel geometries to achieve desired solute width distributions over space and time. We apply the latter analysis to generate a geometry that results in a constant axial width and a second geometry that results in a sinusoidal axial variance in space.


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 1993).G. I. Taylor first reported a closedform solution for the long-term dispersive behavior of a solute in a circular cylindrical tube driven by a fully-developed Poiseuille flow (Taylor 1953).Taylor's 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 / ≫   ≫ 6.9, where  is the characteristic channel length,  is the radius of the tube, and   is a Peclet number based on .Aris later included the effects of axial molecular diffusion of the solute and introduced the method of moments which enabled treatment of all stages † Email address for correspondence: juan.santiago@stanford.eduarXiv:2211.09255v2[physics.flu-dyn]20 Nov 2023 of the dispersive process in long-thin channels with fairly arbitrary cross-sections (Aris 1956).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 analyzing the dispersion behavior in various systems.The book of Brenner & Edwards includes a broad range example analyses including dispersion of flows through porous media and effects of chemical reactions and surface adsorption (Brenner & Edwards 1993).
A significant emphasis has been placed on the dispersion problem in spatially-periodic channel geometries.The seminal work of H. Brenner (which is later referred as Brenner-Aris theory, or the macrotransport paradigm) (Brenner 1980) provided a generalized framework to predict the long-term dispersion behavior of point-size particles in any spatially periodic channels or porous media.This approach solves an elliptical partial differential equation of a cell field B (commonly referred as B field) defined on a periodic unit cell and uses this to compute a long-term dispersion tensor.Hoagland & Prud'homme (Hoagland & Prud'Homme 1985) applied Brenner-Aris theory to the dispersion in a longthin, 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 et al. performed a similar analysis to Hoagland & Prud'Homme in axisymmetric channel with a sinusoidal variation of radius, and extended the results by comparing the predictions with Brownian dynamics simulations (Bolster et al. 2009).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 2007;Dorfman 2008Dorfman , 2009)).Adrover et al. later considered dispersion of both point and finite-sized particles in a long-thin sinusoidal channel (Adrover et al. 2019).For point particles, they obtained various asymptotic expressions for the effective, long-term dispersion coefficient in the limits of large and small   .
Methods other than Brenner's theory were also proposed to analyze the dispersion in periodic channels.For example, Rosencrans (Rosencrans 1997) followed the center manifold theory of Mercer and Roberts (Mercer & Roberts 1990) 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. 2011(Martens et al. , 2013b(Martens et al. ,a, 2014)).
Several studies also addressed the dispersion problem in spatially non-periodic systems.Early work from Saffman analyzed the long-term longitudinal dispersion in the direction of mean flow in a random porous media where the flow satisfies Darcy's law (Saffman 1959).Gill, Ananthakrishnan and Nunge analytically and numerically analyzed the dispersion in a velocity entrance region where the flow field is not fully developed (Gill et al. 1968).Gill and Guceri later investigated numerically the dispersion in diverging channels over a wide range of angles of divergence and Peclet number (Gill & Güceri 1971).Smith described the longitudinal dispersion in a varying channel in terms of the memory effect of the dispersion coefficient.Smith 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 (Smith 1983).Mercer and Roberts used center 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 behavior (Mercer & Roberts 1990).Horner, Arce and Locke 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 (Horner et al. 1995).
Bryden & Brenner analyzed the Taylor-Aris dispersion in a diverging conical channel and a flared, axisymmetric Venturi tube geometry (Bryden & Brenner 1996).They used multiple-timescale 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 analyzed the long-term dispersion and growth of a solute zone.Stone and Brenner also investigated the dispersion in a different spatially non-periodic system.The latter work considered dispersion in a radially expanding flow between large parallel plates.The latter flow results in a mean velocity which varies in the streamwise direction and therefore exhibits a location-dependent dispersion coefficient (Stone & Brenner 1999).
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 Dankwerts (Danckwerts 1953), 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 analyzing and benchmarking results from numerical models (Rodrigues 2021).Although the formulations are different, the spatial moments and temporal moments of the dispersion process have been shown to be interrelated (Vikhansky & Ginzburg 2014;Ginzburg & Vikhansky 2018).
To our knowledge, all previous analyses focused on the long-term dispersion behavior and were not able to provide a simple prediction of both the short-term and longterm spatial evolution of solute.We also know of no analytical work toward 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 analyze 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.

Taylor-Aris dispersion in arbitrarily shaped cylinders
We analyze Taylor-Aris dispersion for flow in an axisymmetric channel with a slowly varying, arbitrary radius  = ().We define variables for the first and second derivatives as () = d/d and () = d 2 /d 2 .Given the azimuthal symmetry, the concentration (, , ) of a solute evolves according to the following convection-diffusion equation in Key parameters used for scaling analysis:  0 : radius of the channel at  = 0, taken as the characteristic radius. 0 : cross-sectional averaged axial velocity at  = 0, taken as the characteristic axial velocity.: characteristic wavelength of spatial variation of the channel.
Additional assumptions and asymptotics needed to derive Eq. 2.26 and 2.29 from Eq. 2.19: (1):  we write the velocity field as where  is the Reynolds number,  is our primary smallness parameter defined as the ratio between the characteristics radius  0 (defined as the radius at  = 0) and the characteristic axial length scale of radius variation .The second terms on the righthand side (RHS) indicate the order of magnitude of the truncation error associated with this asymptotic approximation. 0 is the characteristic mean axial velocity defined here as the mean axial velocity at  = 0. We also assume that the characteristic width of the solute zone  is much smaller than the characteristic wavelength of radius variation, or / ≪ 1.For the Taylor-Aris analysis, we follow the notation of Stone and Brenner (Stone & Brenner 1999) and expand each variable into cross-sectional averages of the form Whence, the area-averaged velocity components are (2.4) We then expand the convective-diffusion equation as (2.5) The area-average of the latter equation is then We simplify four of these terms in Eq. 2.6 using Leibniz's rule, integration by parts, and chain rule as follows: (2.9) (2.10) Inserting these terms to equation 2.6, we have Subtracting Equation 2.11 from Equation 2.5 yields the following equation for the deviation of concentration: (2.12) Table 2. Ranked summary of terms in Equation 2.12.
Left-hand side (LHS) Right-hand side (RHS) To determine the dominant terms on the left and right hand side of the equation, we perform a scaling analysis with the following scales: , and   =  0  *  .We define  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. 0 and  ′ 0 are the characteristic scales of the area-averaged solute concentration and its deviation, respectively.Our scaling assumes 4 smallness parameters:  =  0 / ≪ 1,  =  0 / ≪ 1, / ≪ 1, and  ′ 0 / 0 ≪ 1.We summarize the order of magnitude of each term in Equation 2.12 in Table 2.
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 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, ▽ • n = 0.This requires Substituting the boundary condition into Equation 2.13 and keeping only terms of the same order as the dominant terms in Table 2, we have (2.14) We next directly integrate Equation 2.14 as Here, the function  1 (, ) results from the partial integration with respect to .Evaluating this expression at  = 0, we see   ′  =0 =  1 = 0, whence  1 (, ) = 0. Integrating the equation a second time we obtain By the definition of the fluctuating quantity, ∫  0  ′ d = 0. We use this to solve for  2 and obtain an expression of We next differentiate Eq. 2.17 with respect to , multiply by the known function  ′  , and take an area integral to construct the term  ′   ′  .We perform similar procedures for each term in Eq. 2.11, and we arrive at Equation 2.18.We rank the terms in Eq. 2.18 and summarize in Table 3. (2.18) The detailed derivation of Eq. 2.18 is summarized in Supplementary Section A. Keeping terms on the order of  0  O   0 and  0  2 O   yields the following expression for the development of the area-averaged concentration: which can also be rearranged and written in this conservative form (Stone December 14, 2022): (2.20) Note that, for  = 0, Equation 2.19 reduces to the simple result of Taylor-Aris dispersion for laminar, fully-developed flow within cylindrical channel (with uniform-radius) (Aris 1956).
In the next section, we will use Equation 2.19 to develop a description of the development of the axial mean position of the solute and its axial variance.For now, we note the RHS of Equation 2.19 contains a prefactor for the second axial derivative  2 ⟨⟩  2 which has the form of a typical Taylor-Aris dispersion coefficient: In a cylindrical tube with constant radius, the axial variance grows linearly in time according to 2 eff , as  =  = 0 and the mean velocity correction due to   ≠ 0 vanishes.However, we note that the coefficient  eff is by itself not useful in predicting the time Table 3. Ranked summary of terms in Equation 2.18.
Left-hand side (LHS) ⟨⟩  evolution of axial variance when the cross-section is varying in space.This is true even for the simple case of linearly converging or diverging channel (i.e.,  = 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 (LHS).These advective gradients can stretch or shrink the variance as the solute navigates through contractions and expansions, respectively.

Time evolution of mean and variance
We next formulate the problem in terms of analytical expressions of two moments of solute distributions.Note that, despite this use of spatial moments, our analysis does not follow Aris' famous method of moments (Aris 1956) as we here 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 which are not part of the method of moments.
We first define an axial average as where  ′ is a dummy variable for integration.For more compact notation, we also define  as the axial integration of the mean concentration To derive the first moment, we multiply Equation 2.19 with  and apply the axial We evaluate integrals using integration by parts, divide both sides with , and arrive at an ODE for the axial mean position () as follows: The two terms on the RHS are expressions of the form  ().Some analysis of this term is essential for further analytical simplification of the problem.We consider a Taylor expansion of  () about the mean axial position  as follows: To simplify the expressions, we define  ≡  ().Insert the expression in Eq. 2.25 into Eq.2.24 and keep the dominant term, we derive the following ODE for (): which reduces to the standard Taylor result for () = 0.
To derive the second moment equation of Equation 2.19, we multiply Equation 2.19 by ( − ) 2 , apply the axial average operation, perform integration by parts, and divide both sides by .This then yields the following ODE for the dynamics of  2  : Note that the RHS of Equation 2.27 includes expressions of the form  (), ( − )  (), and   ().We again simplify the latter second and third terms using Taylor expansion at the mean axial position  as follows: (2.28) We insert the expression in Eq. 2.28 into Eq.2.27, keep the term of order  × O (1), 0  , and arrive at an ODE of the axial variance  2  () as (2.29) In the latter equation, as / approaches unity, the dominant error term is For / approaching unity, this has the following order of magnitude: Here Skew  is the axial skewness of the solute zone, defined as Skew  = ( − ) 3 / 3  .Detailed derivation of Eq. 2.30 is summarized in Supplementary Section B. Note that Equation 2.29 also reduces to standard Taylor results when () = 0. Equation 2.26 and 2.29 are two coupled ODEs that predict the time evolution of axial mean and variance of solute based on channel geometry and the assumption of lubrication theory in velocity field.In Section 3, we will use these two equations to predict the time evolution of the axial mean and variance of solute zone for various interesting channel geometries.

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, qualitatively, characterized by a positive and sufficiently large value of  and a sufficiently large   (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 here explore this effect.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 Eq. 2.29, we see that 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 ,   , and  2  / 2 .Note that here the negative growth of variance is a transient phenomenon, and is limited simply to 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 (Fig. 2A) delineates the regions of transient positive and negative growth of axial variance as a function of local Peclet number (  ), local slope of channel radius distribution (), and local ratio between variance and squared radius ( 2  / 2 ). Figure 2B shows the contours of horizontal crosssection of the surface for various values of ln( 2  / 2 )(labeled 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 toward the   - plane for finite   and larger values of .We also see that for fixed and finite values of (positive) , setting the LHS of Eq. 2.31 to zero reveals the critical variance value  2  / 2 required for transient negative growth of axial variance as follows: We now differentiate the latter expression with respect to   , set it equal to zero, and arrive at the minimal variance required for transient negative axial variance growth: (2.33) The minimum variance occurs at   = √ 48.

Approximations useful for channel design
The analytical expression of Equation 2.31 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 following approximation: Inserting this approximation into Equation 2.29, and substituting the definition of , we obtain an ODE for () given a desired spatial distribution  2  (): We can then numerically solve Equation 2.35 to obtain the shape of the engineered channel.In Section 3.4, We will use this approximation to design channels that can result in specific variance evolution pattern in space.

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 5,000 point-particles which were initially uniformly distributed along the radius and normally distributed along the streamwise direction.We considered a mean axial position of zero (by definition) and an initial axial variance of  2 ,0 .We set  2 ,0 = 300 2 0 for the two engineered channels in Section 3.4, and we set  2 ,0 = 10 2 0 for all other cases.The velocity field follows Equation 2.2.The system evolves according to the forward Euler method.
For the time steps, we set the ratio between diffusion time scale ( 2 0 /) to time step of discretization to 40.Each reflection at the wall was checked twice at each time step.Unless specified, the (constant, initial) Peclet number based on initial radius (  0 ) was 10.For the three periodic channels shown in Section 3.2, the initial Peclet 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  () is computed with  () =  particle =1  (  )/ particle , where   is the axial location of the -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 Supplementary Section C.

Diverging and converging conical channels
We first apply our analysis to a simple case of diverging and converging conical section channels.Figures 3A and 3B shows solutions for dispersion of a solute zone in linearly diverging and converging channels, respectively.The top two figures in Figure 3 show particle zones at five non-dimensional times 2/ 2 0 (from 5 to 800) as they migrate through the channels.The axes are non-dimensionalized coordinates  * = / 0 and  * = / 0 .Note the intensely exaggerated height-to-length aspect ratio of the figures, chosen here for clarity of presentation.The top half of the solute concentration in the channel are raw data results of the three-dimensional Brownian dynamics simulations.The bottom half 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 Eqs.2.26 and 2.29).Note that we present and plot the predicted axial distribution assuming a Gaussian distribution as we only have the information from mean and variance.However, the derivation of Eqs.2.26 and 2.29 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 width of the axial solute distribution are each set equal to the same value (equal to √ 10 0 = √ 10( = 0) of the diverging channel).The current model has very good comparison with the Brownian dynamics simulations.The characteristic   values here are order 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 symmetric about  * .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 shall explore such cases later below.) The bottom two plots in Figure 3 are plots of the axial variance of the area-averaged axial distribution of solute concentration from the current model ( * 2 ,curr ) and Brownian dynamics ( * 2 ,B ).These quantities are plotted as a function of the axial mean position of solute zone,  * .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.2.29) and from Brownian dynamics simulation (subscript "B").Axial variance computed using Equation 2.29 shows excellent agreement with Brownian dynamics simulations.

Periodic channels
We next apply our model to a periodic channel with sinusoidal radius distribution.Similar to the top plots of Figure 3, the top plot of Figure 4 shows a plot of the channel geometry in coordinates of  * versus  * in a highly exaggerated aspect ratio for clarity.We show results for   0 of 10 where  0 is defined as the radius at  = 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 half of the solute zones are raw results from the Brownian dynamics simulations, while the bottom half 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  * = 1800, and then the subsequent rapid decrease in axial variance caused by the advective effects of the solute expanding from  * = 1800 to just upstream of  * = 3000.Note that the particle clouds from Brownian dynamics are also fairly symmetric about  * .
The middle plot of Figure 4 shows the axial variance of the solute zone as a function of solute average axial location  * .Plotted are the axial variance from the Brownian dynamics simulation and the current model.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 analyze period-averaged axial variance over long times.
Figure 4 also shows the non-dimensional area-averaged flow velocity and the local normalized effective dispersion coefficient,  eff / (c.f.two ordinances on the RHS of the plot).Both  and  * eff are periodic functions, antiphase to the shape of the channel. * 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  * = 2000 to 3000) the effective dispersion coefficient  eff / remains positive and greater than unity.This demonstrates the inability of effective dispersion coefficient alone to describe the axial variance evolution (including here in a periodic channel).The bottom plot of Figure 4 shows the observed and predicted error in dimensionless growth rate of axial variance of solute as a function of  * .Also shown in the figure is the ratio between the square root of axial variance and channel wavelength   /.Consistent with the middle plot of Figure 4, the observed error in axial variance growth rate is small and stochastic, and it matches with the predicted error from Eq. 2.30.The ratio   / 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 Peclet 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 Peclet number to 100.Similar to the top plots of Figure 4, the top plot of Figure 5 shows a plot of the channel geometry in coordinates of  * versus  * 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 half of the solute zones are raw results from the Brownian dynamics simulations, while the bottom half 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 through contractions and expansions, respectively.Our model tends to overestimate the fluctuation in axial variance, especially when the solute transverses through 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.
The middle plot of Figure 5 shows the axial variance of the solute zone as a function of  * .Plotted are the variance from the Brownian dynamics simulation and the current model.There is a good agreement between these two models in trends of development of the variance, and an excellent agreement when  * 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 contraction.
The bottom plot of Figure 5 again shows the observed and predicted error in di- eff are periodic functions and shown as a reference.The bottom plot shows the predicted (Eq.2.30) 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 the axial variance and channel wavelength   /.Variance computed using Equation 2.29 shows excellent agreement with Brownian dynamics simulations.Note 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 Eq. 2.30.This shows the accuracy of our model when   / is small.mensionless growth rate of axial variance of solute as a function of  * .Also shown in the figure is the ratio between the square root of axial variance and channel wavelength   /.The magnitude of observed error grows as the axial variance grows, and the error is most pronounced when the solute transverses through contractions.There is an excellent agreement between the predicted error from Eq. 2.30 when  * is less than 20000, but there is more deviation between the two predictions as   / grows larger.
We further increases the initial Peclet number of the system to 1000 using the same channel of Figure 4 and Figure 5. Similar to the topmost plots of Figure 5, the top of Figure 6 shows a plot of the channel geometry in coordinates of  * versus  * in a highly exaggerated aspect ratio for clarity.Again, the shown are solute zones at non-dimensional times ranging from 5 to 80.The top half of the solute zones are raw results from the Brownian dynamics simulations, while the bottom half are Gaussian distributions with ,B and  * 2 ,curr , each as a function of the non-dimensional axial location along the channel, / 0 .The bottom plot shows the predicted (Eq.2.30) 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   /.Variance computed using Equation 2.29 shows excellent agreement with Brownian dynamics simulations when the channel is expanding, but overestimates the axial variance when the channel is contracting.Note variance averaged along the axial spatial period increases monotonically as expected.The predicted error in axial variance growth rate using Equation 2.30 shows excellent agreement with the observed error, especially when   / is small.As   / increases, there is more deviation between observed and predicted error in axial variance growth rate.
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 on the same order as the wavelength of the channel.
The middle plot of Figure 6 shows the axial variance of the solute zone as a function of  * .Plotted are the variance from the Brownian dynamics simulation and the current model.The agreement between the two is limited to the initial development of the axial variance ( * < 2000).The axial variance growth in Brownian dynamics simulation becomes monotonic after  * > 30000, but the current model is not able to capture this.
The bottom plot of Figure 6 again shows the observed and predicted error in dimensionless growth rate of axial variance of solute as a function of  * .Also shown in the ,B and  * 2 ,curr , each as a function of the non-dimensional axial location along the channel, / 0 .The bottom plot shows the predicted (Eq.2.30) 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   /.Variance computed using Equation 2.29 only shows agreement with Brownian dynamics simulations in the very beginning.Our model fails to predict the monotonic growth of axial variance when the axial variance on the same order as the channel wavelength.The predicted error in axial variance growth rate using Equation 2.30 shows fair agreement with the observed error for  * below 10000.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 is the ratio between the square root of axial variance and channel wavelength   /.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 Eq. 2.30 can only capture the observed error in the initial development of axial variance ( * < 10000), but becomes inaccurate as   / grows above 0.3.This finding suggests that the accuracy of our model depends less on the Peclet number itself but more sensitive to the ratio   /.At large Peclet number, since the axial variance growth rate is much faster (as discussed in Figure 7) and   / quickly approaches and exceeds unity, the current model becomes inaccurate.A zoomed-in comparison of the solute zone subject to Peclet number of 10, 100 and 1000 is in Figure S1.
We next consider the long-term averaged dispersion coefficient in periodic channels.As with previous studies of dispersion (Hoagland & Prud'Homme 1985;Adrover et al. 2019), we will define a new effective dispersion coefficient for the long-term (many period) growth of the axial variance  ∞ eff defined as follows: The non-dimensional version of this quantity will be defined as  ∞ * eff ≡  ∞ eff /.We analyzed long-term dispersion behavior 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   / is greater than about 0.3.We thus computed the long-term growth of the axial variance from our current model as follows: Figure 7 summarizes the results of this study.The three plots in the top row show plots of the radius distribution () normalized by initial radius  0 as a function of the normalized axial coordinate  * .The three functions () are sinusoidal, triangular, and exponential-sine-wave functions of the following forms: where  () is the Heaviside step function.The three plots in the second row of Figure 7 show the effective, long-term dispersion coefficient  ∞ * eff for each case as a function of   0 . ∞ * eff curves are shown for the current model (data point) 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  ∞ * eff = 1 for vanishing   0 , as expected. ∞ * eff increases monotonically with increasing   0 and asymptotes to a  2  0 dependence for large   0 .In the left-hand plot of the bottom row, we also show a comparison of our model results with the work of Adrover et al. (Adrover et al. 2019).Adrover analyzed the long-term dispersion of solutes in a sinusoidal channel and provided an approximate analytical formula for  ∞ * eff as a function of   0 , and this prediction is plotted along with our model in the figure.Adrover's also found that  ∞ * eff tends to unity for small   0 and scales with the second power of   0 for large   0 .We note the excellent agreement among the three predictions.
We conclude the current model can be readily adapted to a wide variety of periodic channel shapes.We further note the similarity among the three  ∞ * eff versus   0 curves for the three cases.This similarity leads us to hypothesize that the long-term solute dispersion in such periodic channels is largely driven by the spatial frequency and the amplitude of the radius oscillation (and   0 ) and may be insensitive to the details of channel shapes.

Arbitrarily-shaped channels
We next demonstrate a novel application of our model to a particular but arbitrary axisymmetric channel shape.Similar to the top plot of Figure 4, the top plot of Figure 8 eff computed using Eq.2.29 and with a Brownian dynamics simulation.The plot on the left also shows a comparison with the expression derived by Adrover et al. (2019) for a sinusoidal channel (Adrover et al. 2019).Note 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.
shows a plot of the channel geometry in coordinates of  * versus  * in a highly exaggerated aspect ratio for clarity.We show results for   0 of 10 where  0 is the radius at  = 0.The channel shows plots of solute zones at non-dimensional times ranging from 25 to 250.The top half of the solute zones are raw results from the Brownian dynamics simulations, while the bottom half 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 through contractions and expansions, respectively.For example, note the sudden decreases in variance due to the expansion just upstream of  * = 500, and then the sudden increases in variance caused by the contraction near  * = 700.
The middle plot of Figure 8 shows the axial variance of the solute zone as a function of nondimensional solute average axial location,  * .Plotted are the axial variance from the Brownian dynamics simulation and the current model.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  * = 500 and  * = 800, and the sudden increases of axial variance just upstream of  * = 700.
The middle plot of Figure 8 also shows the non-dimensional area-averaged flow velocity and the local normalized effective dispersion coefficient,  eff / (c.f.two ordinances on the RHS).Similar to the case in Figure 4, both  and  * eff show opposite trend to the shape of the channel. * eff again increases as the channel constricts and decreases as the channel expands, as we expected.However, even in a region where axial variance is transiently decreasing (e.g. between  * = 400 to 500), the effective dispersion coefficient remains positive and greater than unity.This again highlights the inability of effective dispersion coefficient to describe the variance evolution, here in an arbitrarily-shaped axisymmetric channel.
The bottom plot of Figure 8 shows the scaled rate of change of variance ( 1  d 2  d ) as a function of solute average axial location  * .Plotted are the results computed from Brownian simulations and the current model (Eq.2.29).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 the top plot of the same figure), 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  * = 400 and  * = 700.This shows the capability of our model to predict the axial variance evolution, and also provides validation of using Equation 2.31 to identify the regime of transient negative axial variance growth.

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 Equation 2.35.The shape of channel A monotonically diverges, and the rate of diverging increases downstream.For channel B, there is an increased diverging rate between  * = 300 and 500, coinciding with the region where it is necessary to reduce the variance according to the targeted sinusoidal pattern.Note that the channel shape in Figure 9A has an analytical expression.By setting the d 2  d term on the RHS of Equation 2.35 to 0, we can simplify Equation 2.35 into the following ODE for (): which has an analytical solution as follows: where  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.
Figures 9A and 9B show solutions for dispersion of solute zone in the two engineered channels.The channel in Figure 9A (abbreviated as "channel A") is designed to maintain an approximately constant axial variance of 300 ( 2 *  ( * ) = 300), while the channel in Figure 9B (abbreviated as "channel B") is designed to yield a sinusoidal axial variation of axial variance ( 2 *  ( * ) = 300 + 50 sin(2 * /600)).Similar to the top plot of Figure 3, the top two plots of Figure 9 shows the two engineered channel geometry in coordinates of  * versus  * .We show results for   0 of 10 where  0 is taken as the radius at  = 0.The top half of the solute concentration in the channel are raw data results from Brownian dynamics simulations, while the bottom half of the channels show the area-averaged axial distribution of solute concentration predicted from the current model.
The bottom two plots in Figure 9 are plots of the axial variance of the area-averaged axial distribution of solute concentration from the current model ( * 2 ,curr ), Brownian dynamics ( * 2 ,B ), and the targeted axial variance spatial evolution pattern ( * 2 ,target ).These quantities are plotted as a function of the axial mean position of the solute zone,  * .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 a variance of 300 as we 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 Section 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. 2016).The most common method of detection of analytes is an optical detection technique such as fluorescence, colorimetric, or UV adsorption (Wu & Gu 2011).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 areaaveraged 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.

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 areaaveraged 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 then 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 characteristics wavelength of channel variations.
We further derived an analytical expression which 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 Peclet number, 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 which yields 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 which 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 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 Peclet 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 analyzed the long-term (many period) developments of solute in periodic channels by defining a long-term effective dispersion coefficient.We analyzed three separate periodic channel shapes and showed this long-term behavior 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 which 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.
Lastly, 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 which 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.
Taylor dispersion in arbitrarily shaped axisymmetric channels

D. Supplementary figures
Figure S1 shows plots of the Brownian dynamics simulations (top) and the concentration distribution weighted by local cross-sectional area (bottom) for the channel with a sinusoidal radius distribution (the same channel as in Figure 4-6), with initial Peclet number   0 varying from 10, 100 to 1000 and the mean axial location  the same.As expected, the solute zone subjected to a higher Peclet number spreads faster, and the axial concentration distribution is more skewed.

E. Comparison of effective dispersion coefficient in periodic channels
with solutions by Adrover et al.
Adrover et al. comprehensively analyzed effective dispersion coefficient in a wide range of wavelengths for a periodic channel with a specific radius () =  0 (1+ sin(2/( 0 ))) (?).For our current model, the Taylor dispersion analysis requires that the flow field in the channel be well approximated by lubrication theory.This constraint in turn requires the ratio () of the channel wavelength to channel radius be significantly greater than unity.In this limit, Adrover et al. provided  Here,  2 =  −2 and   =   0  is the Peclet number defined based on wavelength.

Figure 1 .
Figure1.Schematic of an axisymmetric channel with a slowly varying, arbitrary distribution of radius, ().A nominal radius is taken as  0 at  = 0.The slope and the curvature of the cylinder wall are respectively () and (), as shown.  is the characteristic width of a solute zone and  () is the area-averaged axial velocity distribution.The (-direction) length of the channel depicted schematically in the sketch has been compressed relative to its characteristic radius for clarity of presentation.

Figure 2 .
Figure 2. Regimes of transient positive and negative growth of axial variance of solute.(A) Surface of zero variance growth in a space of Peclet number, , and the natural log of the local ratio between variance and squared radius (ln( 2 / 2 )).This surface is computed analytically from Equation2.31.Shown is a surface where the axial variance rate of growth, d 2  d , is approximately 0. (B) The contour plot shows the horizontal cross-section curves of the zero transient axial variance growth surface in a space of Peclet number, , and ln( 2  / 2 ), at different ln( 2  / 2 ) values (labeled on each contour line).

Figure 3 .
Figure 3. Taylor-Aris dispersion in (A) diverging and (B) converging conical channels.The plots in the top row show results from Brownian dynamic simulation (upper half) along with the predicted axial distribution of area-averaged concentration (bottom half).The bottom two plots show solute axial variance as a function of  * .Plotted are axial variance from the current analytical model (subscript "curr", Equation2.29)and from Brownian dynamics simulation (subscript "B").Axial variance computed using Equation2.29 shows excellent agreement with Brownian dynamics simulations.

Figure 4 .
Figure 4. Taylor-Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, ().The top plot shows results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (Equation 2.29).The middle plot shows the distribution of  * 2 ,B ,  * 2 ,curr ,  * (), and  eff /, each as a function of the non-dimensional axial location along the channel, / 0 . * () and  *eff are periodic functions and shown as a reference.The bottom plot shows the predicted (Eq.2.30) 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 the axial variance and channel wavelength   /.Variance computed using Equation2.29 shows excellent agreement with Brownian dynamics simulations.Note 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 Eq. 2.30.This shows the accuracy of our model when   / is small.

Figure 5 .
Figure 5. Taylor-Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, ().The channel is the same as the one in Figure 4, but the initial Peclet number   0 is 100.The top plot shows results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (Eq.2.29).The middle plot shows distribution of  * 2,B and  * 2 ,curr , each as a function of the non-dimensional axial location along the channel, / 0 .The bottom plot shows the predicted (Eq.2.30) 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   /.Variance computed using Equation2.29 shows excellent agreement with Brownian dynamics simulations when the channel is expanding, but overestimates the axial variance when the channel is contracting.Note variance averaged along the axial spatial period increases monotonically as expected.The predicted error in axial variance growth rate using Equation2.30shows excellent agreement with the observed error, especially when   / is small.As   / increases, there is more deviation between observed and predicted error in axial variance growth rate.

Figure 6 .
Figure 6.Taylor-Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, ().The channel is the same as the one in Figure 4, but the initial Peclet number   0 is 1000.The top plot shows results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (Eq.2.29).The middle plot shows distribution of  * 2,B and  * 2 ,curr , each as a function of the non-dimensional axial location along the channel, / 0 .The bottom plot shows the predicted (Eq.2.30) 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   /.Variance computed using Equation2.29only shows agreement with Brownian dynamics simulations in the very beginning.Our model fails to predict the monotonic growth of axial variance when the axial variance on the same order as the channel wavelength.The predicted error in axial variance growth rate using Equation2.30shows fair agreement with the observed error for  * below 10000.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 7 .
Figure 7. Normalized long-term effective dispersion coefficient  ∞ * eff as a function of   0 for three different periodic channels with equal period and similar radius amplitude ( = 0.05,  = 200).All plots show the  ∞ * eff computed using Eq.2.29 and with a Brownian dynamics simulation.The plot on the left also shows a comparison with the expression derived by Adrover et al. (2019) for a sinusoidal channel(Adrover et al. 2019).Note 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.

Figure 8 .
Figure 8. Taylor-Aris dispersion for an example arbitrarily shaped axisymmetric channel.The top plot shows results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (Equation 2.29).The middle plot shows distribution of  * 2 ,B ,  * 2 ,curr ,  * (), and  eff /, each as a function of the non-dimensional axial location along the channel. * () and  eff / are shown as a reference.Variance computed using Equation2.29 shows excellent agreement with Brownian dynamics simulations.The bottom plot shows the same example used for demonstration and benchmarking of Equation2.31.Plotted is a scaled rate of change of axial variance as a function of mean solute position  as computed using Equation2.29 and analytical expression Equation2.31.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 a windows size of Δ * = 2.5, while the light blue curve is the original Brownian simulation results.

Figure 9 .
Figure 9. Engineering the axial variance evolution in Taylor-Aris dispersion.Using Equation 2.35, we designed two channels 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.The top plot shows results from a Brownian dynamic simulation (upper half) and axial solute distribution predicted with the current model (Equation 2.29).The bottom plot shows distribution of  * 2 ,B ,  * 2 ,curr , and  * 2 ,target , each as a function of the non-dimensional axial location along the channel.Axial variance computed using Equation 2.29 and axial variance computed from Brownian dynamic simulation both show excellent agreement with the targeted variance evolution pattern.

Figure S1 .
Figure S1.Taylor-Aris dispersion in a channel with a sinusoidal (periodic) radius distribution, ().The channel is the same as that of Figure 4-6, with initial Peclet numbers   0 of 10 (blue,  = 500 2 0 /(2)), 100 (orange,  = 50 2 0 /(2)) to 1000 (green,  = 5 2 0 /(2)).The top plot shows results from a Brownian dynamic simulation.The bottom plot shows the concentration distribution normalized by local cross-sectional area.For a fixed distance of travel, the solute zones subjected to greater Peclet numbers disperse faster and the solute distribution becomes more skewed.

Table 1 .
Summary of assumptions made in this study.