1. Introduction
Understanding the dispersion of solutes in fluids flowing through channels is crucial for designing devices and optimising processes, such as the transport of nutrients and medications in blood vessels, the dispersion of contaminants in aquifers, industrial chemical reactions and lab-on-a-chip systems (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Kirby Reference Kirby2010; Wu & Gu Reference Wu and Gu2011). Many of these applications involve non-Newtonian fluids, especially those exhibiting viscoelastic behaviour, which are common in chemical engineering, physiological flows (Thurston Reference Thurston1972) and environmental mechanics. Given the wide-ranging significance of these systems, it is essential to understand solute dispersion, particularly in channels with non-circular or varying geometries, for practical applications.
The foundational work of Taylor (Reference Taylor1953) established the theoretical framework for dispersion that describes the long-term behaviour of a solute in a solvent subjected to slow, pressure-driven (Poiseuille) flow through a cylindrical tube. Taylor demonstrated that, at large time scales, axial molecular diffusion becomes negligible compared with advection, while radial diffusion remains significant. This leads to a one-dimensional diffusion equation for the mean concentration. Aris (Reference Aris1956) later generalised this analysis by removing parameter restrictions and introduced the method of moments to track the evolution of mean concentration and variance. Together, these studies form the basis of Taylor–Aris dispersion, which is applicable when radial diffusion time scales are much shorter than advective time scales.
Subsequent research has expanded the analysis of the dispersion phenomenon to diverse flow conditions, including oscillatory, pulsatile and electrokinetically driven flows. For example, Chu et al. (Reference Chu, Garoff, Przybycien, Tilton and Khair2019a ,Reference Chu, Garoff, Tilton and Khair b ) and Barik & Dalal (Reference Barik and Dalal2019) reported that oscillatory flows can yield negative dispersion coefficients two orders of magnitude larger than steady flows with identical amplitudes due to the reversing flow. Similarly, electroosmosis (Muñoz et al. Reference Muñoz, Arcos, Bautista and Méndez2018; Mukherjee et al. Reference Mukherjee, Dhar, DasGupta and Chakraborty2019; Teodoro, Bautista & Méndez Reference Teodoro, Bautista and Méndez2019; Mederos et al. Reference Mederos, Arcos, Bautista and Méndez2020) and diffusioosmosis (Nakad et al. Reference Nakad, Witelski, Domec, Sevanto and Katul2021; Alessio et al. Reference Alessio, Shim, Gupta and Stone2022) have been investigated as mechanisms for solute dispersion.
As mentioned earlier, the dispersion process can occur in flows involving fluids that exhibit non-Newtonian behaviour. A study conducted by Fan & Wang (Reference Fan and Wang1966) analysed the dispersion process and determined effective dispersion coefficients for Bingham plastic and Ellis models based on the classic Taylor–Aris analysis, finding a mathematical model that allows recovery of the Newtonian and the power-law models. Akyildiz (Reference Akyildiz2002) examined viscoelastic effects using the Phan-Thien–Tanner (PTT) model, considering the linearised and exponential functions of the trace of the stress tensor, and demonstrated its direct impact on solute distribution. Other rheological models, such as the Casson fluid (Naragani & Sebastian Reference Naragani and Sebastian2012) and the Carreau–Yasuda fluid (Rana & Murthy Reference Rana and Murthy2016), have also been investigated. Effective transport coefficients were derived and an asymptotic analysis revealed the behaviour of solute dispersion. It was found that the Carreau fluid behaves like a Newtonian fluid at low shear rates, whereas the other fluids exhibit significant differences during the solute dispersion.
However, the geometry of the channel plays a crucial role in influencing dispersion, particularly when the geometries are irregular (Locascio Reference Locascio2004). Pioneering research by Brenner (Reference Brenner1980) established a foundational framework for analysing dispersion in spatially periodic channels. Subsequent studies have adapted this framework to predict long-term dispersion behaviour in various spatially periodic channels and porous media. For example, Hoagland & Prud’Homme (Reference Hoagland and Prud’Homme1985) numerically examined dispersion along the streamwise direction in sinusoidal channels using the method of moments for arbitrary wavelengths. Bolster, Dents & Borgne (Reference Bolster, Dents and Borgne2009) expanded on this by studying dispersion in a sinusoidal axisymmetric channel and comparing their predictions with Brownian dynamics simulations. Additionally, Dorfman (Reference Dorfman2009) included electrokinetic effects in their analysis of thin sinusoidal channels. Research has also been conducted on non-periodic geometries. Smith (Reference Smith1983) explored longitudinal dispersion in channels with varying geometries, focusing on changes in dispersion coefficients across different geometric profiles. Bryden & Brenner (Reference Bryden and Brenner1996) analysed Taylor–Aris dispersion in a diverging conical channel using multiple-time scale analysis to derive an asymptotic expression for the effective dispersion coefficient. Similarly, Stone & Brenner (Reference Stone and Brenner1999) studied dispersion in a different spatially non-periodic system, examining radial dispersion in a flow between large parallel plates that are expanding.
Recent advances in the field include temporal moment analyses (Vikhansky & Ginzburg Reference Vikhansky and Ginzburg2014) and the work by Chang & Santiago (Reference Chang and Santiago2023) that investigates Taylor–Aris dispersion in an arbitrarily shaped axisymmetric channel with a slowly varying radius. These authors derived an expression for solute dynamics based on channel geometry, specifically for Newtonian fluid flow. Building on their findings, we explore Taylor–Aris dispersion in a viscoelastic fluid flowing through axisymmetric channels with gradual changes in radius. We use the sPPT rheological model proposed by Phan-Thien & Tanner (Reference Phan-Thien and Tanner1977). Our objective is to examine how the viscoelasticity of the fluid and various geometries of the channel cross-section affect the spatial and temporal distribution of the solute. We analyse our results both graphically and mathematically. To achieve the above-mentioned goals, we solve the hydrodynamics driven by a pressure gradient along an axisymmetric channel. We apply methodologies developed by Taylor (Reference Taylor1953) and Aris (Reference Aris1956) to derive a mathematical model based on lubrication theory. This allows us to simplify the model by reducing the convection–diffusion equation to a standard diffusion equation. We introduce a change of variables in a new spatial dimension that moves with the fluid and adapts to changes in channel shape. Consequently, we can solve the simplified model numerically using the Crank–Nicolson method to determine the spatial and temporal distribution of the average solute concentration.
The motivation for this study extends beyond fundamental analysis and is deeply rooted in biophysical applications, like in targeted drug delivery. In physiological systems, numerous pathways – such as the circulatory network, the lymphatic vessels and the complex geometries of the interstitial space – can be idealised as channels with axially varying cross-sections. Crucially, many biological fluids (e.g. blood, mucus, synovial fluid) exhibit non-Newtonian, viscoelastic behaviour, which dramatically alters flow profiles and solute transport compared with Newtonian models. Understanding the coupled effect of geometry and fluid rheology on dispersion is essential for predicting the distribution and efficacy of numerous solutes.
To the best of the authors’ knowledge, no studies in the specialised literature have addressed the analysis of Taylor dispersion in channels with varying cross-sections when using viscoelastic fluid flows, such as those described by the sPTT model, which was chosen as a rheological model in this work for its proven capability to capture the key nonlinear rheology of polymers, including shear-thinning and finite extensional viscosity, ensuring physically accurate and numerically stable solutions for this steady-state flow. This makes our work insightful for a variety of applications.
2. Problem formulation
A qualitative schematic of the physical model being studied is presented in figure 1. We analyse the Taylor–Aris dispersion for the flow in an arbitrarily shaped axisymmetric channel with a slowly varying radius
$a(z)$
in which a viscoelastic fluid flows with an area-averaged axial velocity distribution
$U(z)=(1/\pi a^2)\int _0^a u_z 2\pi r \,{\rm d}r$
, with
$u_z$
representing the axial velocity. A nominal radius and area-averaged axial velocity are taken at
$z=0$
as
$a_0$
and
$U_0$
, respectively. The slope of the channel is
$\beta (z)={\rm d}a/{\rm d}z$
and the second derivative is
$\gamma (z)={\rm d}^2a/{\rm d}z^2$
. The origin of the cylindrical coordinate system (
$z, r$
) is located at the left end and the centre of the channel, as shown in the schematic, and we identify two key axial scales: the characteristic axial length scale of radius variation,
$\lambda _z$
, and the characteristic axial dimension of the solute zone,
$\sigma$
. The solute zone is sufficiently long to be considered a study zone, ensuring that
$\lambda _z/\sigma \ll 1$
. We define two parameters,
$\beta _0 = a_0/\lambda _z\ll 1$
and
$\epsilon =a_0/\sigma \ll 1$
, useful for our mathematical modelling. We consider a stationary flow driven by a pressure gradient such that the flow rate is constant.
A schematic diagram depicts an arbitrarily shaped, axisymmetric channel with a slowly varying radius
$a(z)$
, through which a non-Newtonian fluid flows. The area-averaged axial velocity distribution is
$U(z)$
. At the position
$ z=0$
, the nominal radius and area-averaged axial velocity are denoted as
$a_0$
and
$U_0$
, respectively. The slope of the channel is represented by
$\beta (z)$
, as shown in the diagram. The characteristic axial length scale of radius variation is represented by
$\lambda _z$
and the characteristic axial dimension of the solute zone by
$\sigma$
, which is a length long enough to be considered as a study zone for the solute, so that
$\lambda _z/\sigma \ll 1$
. The axial variance of the solute at a position
$(z)$
is given by
$\sigma _z$
.
2.1. Hydrodynamic field
Consider the equations of mass and momentum conservation:
\begin{equation} \ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0 \end{equation}
and
\begin{equation} \rho \frac {\partial \boldsymbol{u}}{\partial t}+\rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}=-\boldsymbol{\nabla }p+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }+\rho \boldsymbol{g}. \end{equation}
The boundary conditions are the no-slip, impermeability and the velocity symmetry at the channel’s centre. It is assumed that the fluid rheology follows the sPTT model (Phan-Thien & Tanner Reference Phan-Thien and Tanner1977), given by
\begin{equation} f(\tau _{\textit{kk}})\boldsymbol{\tau }+\lambda \left \{ \frac {{\rm D}\boldsymbol{\tau }}{{\rm D}t}-\left[(\boldsymbol{\nabla }\boldsymbol{u})^T\boldsymbol{\cdot }\boldsymbol{\tau }+\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right] \right \}=\eta \left(\boldsymbol{\nabla }\boldsymbol{u}+(\boldsymbol{\nabla }\boldsymbol{u})^T\right)\!. \end{equation}
In (2.1)–(2.3), the velocity vector is represented as
$\boldsymbol{u}=(u_z,u_r)$
, where
$\boldsymbol{\nabla}$
denotes the vector operator in cylindrical coordinates. Here,
$\rho$
stands for the fluid density,
$t$
is the time and
$\boldsymbol{g}$
represents the gravitational vector, which will be ignored in this analysis. The stress tensor is denoted by
$\boldsymbol{\tau }$
, while
$\lambda$
refers to the relaxation time of the viscoelastic fluid. The notation
${\rm D}/{\rm D}t$
represents the material derivative and
$\eta$
is the viscosity coefficient. Additionally,
$f(\tau _{\textit{kk}})$
is a function that depends on the trace of the stress tensor, defined as (
$\tau _{\textit{kk}}$
) given by
$f(\tau _{\textit{kk}})=1+\varepsilon \lambda (\tau _{\textit{kk}})/\eta$
(Phan-Thien & Tanner Reference Phan-Thien and Tanner1977). For small deformations, we can use the linearised form of
$f(\tau _{\textit{kk}})$
. It is important to note that the parameter
$\varepsilon$
is inversely proportional to the extensional viscosity of the fluid.
The initial step in addressing this issue involves establishing the hydrodynamic field from (2.1)–(2.3), while using the lubrication approximation (Leal Reference Leal2007). Comprehensive details are available in Appendix A. Thus, the components of the velocity field are represented by the following equations:
\begin{equation} u_z=2U_n\left ( 1-\frac {r^{2}}{a^{2}}\right )+\frac {3}{2}U_{{v\textit{e}}}\left (1-\frac {r^{4}}{a^{4}} \right ) \end{equation}
and
\begin{equation} u_r=2\beta U_n\left (\frac {r}{a}-\frac {r^3}{a^3} \right )+\frac {3}{2}\beta U_{{\textit {ve}}}\left (\frac {r}{a}-\frac {r^5}{a^5} \right )\!, \end{equation}
where
$U_n(z)$
and
$U_{{\textit {ve}}}(z)$
are the Newtonian and viscoelastic components of
$U(z)$
, respectively, such that
$U(z)=U_n(z)+U_{{\textit {ve}}}(z)$
and are given by
\begin{equation} U_{n}(z)=\frac {a^{2}}{8\eta }\left (-\frac {{\rm d}p}{{\rm d}z} \right ) \end{equation}
and
\begin{equation} U_{{\textit {ve}}}(z)= \frac {\varepsilon \lambda ^2 a^4}{24 \eta ^3}\left (-\frac {{\rm d}p}{{\rm d}z} \right )^3. \end{equation}
Equations (2.6) and (2.7) depend on the pressure gradient
$ {\rm d}p/{\rm d}z$
that is unknown along the channel due to variations in its geometry. However, it can be easily calculated by considering the constriction of a constant flow rate, which means that
$ U a^2 = U_0 a_0^2$
holds throughout the channel.
2.2. Taylor–Aris analysis
The concentration field
$ c(z,r,t)$
of a neutral solute transported by a viscoelastic fluid flowing through the channel shown in figure 1 is governed by the convective–diffusion equation:
\begin{equation} \frac {\partial c}{\partial t}+u_r\frac {\partial c}{\partial r}+u_z\frac {\partial c}{\partial z}=D\left [\frac {1}{r}\frac {\partial }{\partial r}\left ( r\frac {\partial c}{\partial r}\right )+\frac {\partial ^2 c}{\partial z^2} \right ]\!. \end{equation}
In this context, let
$D$
denote the molecular diffusion coefficient of the solute within the fluid. For carrying out the Aris–Taylor analysis, we express any variable
$\theta$
as the sum of two components: one that is averaged over the cross-sectional area, denoted as
$\langle \theta \rangle =(2/a(z)^2)\int _0^{a(z)} \theta r\, {\rm d}r$
, and a fluctuating component
$\theta '$
. We thus represent
$\theta (r,z,t)$
as follows:
\begin{equation} \theta (r,z,t) = \langle \theta \rangle (z,t) + \theta '(r,z,t). \end{equation}
By substituting (2.9) for the variables
$c$
,
$u_r$
and
$u_z$
into (2.8), we obtain the following result:
\begin{align} &\frac {\partial \langle c \rangle }{\partial t}+\frac {\partial c'}{\partial t}+\langle u_r \rangle \frac {\partial c' }{\partial r}+ u_r^{\prime}\frac {\partial c'}{\partial r}+ \langle u_z \rangle \frac {\partial \langle c \rangle }{\partial z}+u_z^{\prime}\frac {\partial \langle c \rangle }{\partial z} + \langle u_z \rangle \frac {\partial c' }{\partial z}+ u_z^{\prime}\frac {\partial c'}{\partial z}\nonumber \\&\quad =D\left [\frac {1}{r}\frac {\partial }{\partial r}\left ( r\frac {\partial c'}{\partial r}\right )+\frac {\partial ^2 \langle c\rangle }{\partial z^2}+\frac {\partial ^2 c'}{\partial z^2} \right ]\!. \end{align}
For the velocity components, we have
\begin{equation} \langle u_z \rangle =U_n+U_{{\textit {ve}}}, \quad u_z^{\prime}=U_n\left ( 1-\frac {2r^{2}}{a^{2}}\right )+\frac {1}{2}U_{{\textit {ve}}}\left (1-\frac {3r^{4}}{a^{4}} \right ) \end{equation}
and
\begin{equation} \langle u_r \rangle =\frac {8\beta }{15}U_n+\frac {12\beta }{21}U_{{\textit {ve}}}, \quad u_r^{\prime}=2\beta U_n\left (\frac {r}{a}-\frac {r^3}{a^3}-\frac {4}{15} \right )+\frac {3}{2}\beta U_{{\textit {ve}}}\left (\frac {r}{a}-\frac {r^5}{a^5} -\frac {8}{21}\right )\!. \end{equation}
2.2.1. Average cross-sectional component of the convective–diffusion equation
Averaging (2.10) over the cross-sectional area, we obtain
\begin{align} &\frac {\partial \langle c\rangle }{\partial t}+\langle u_r\rangle \Bigl \langle \frac {\partial c'}{\partial r} \Bigr \rangle +\Bigl \langle u_r^{\prime} \frac {\partial c'}{\partial z} \Bigr \rangle +\langle u_z\rangle \frac {\partial \langle c\rangle }{\partial z}+\langle u_z\rangle \Bigl \langle \frac {\partial c'}{\partial z} \Bigr \rangle +\Bigl \langle u_z^{\prime} \frac {\partial c'}{\partial z} \Bigr \rangle \nonumber \\&\quad =D \left \lbrace \Bigl \langle \frac {1}{r}\frac {\partial }{\partial r}\left (r\frac {\partial c'}{\partial r}\right ) \Bigr \rangle +\frac {\partial ^2 \langle c \rangle }{\partial z^2} +\Bigl \langle \frac {\partial ^2 c'}{\partial z^2} \Bigr \rangle \right \rbrace \!. \end{align}
We simplify each term of (2.13) by applying Leibniz’s rule, integration by parts and the chain rule, as detailed in Appendix B. By substituting these simplifications into (2.13), we obtain the average cross-sectional component of the advection–diffusion equation:
\begin{align} &\frac {\partial \langle c \rangle }{\partial t}+\frac {2\langle u_r \rangle }{a^2}\left (ac'(a)-\int _0^{a}c' \,{\rm d}r \right )+\Bigl \langle u_r^{\prime}\frac {\partial c'}{\partial r} \Bigr \rangle +\langle u_z\rangle \frac {\partial \langle c \rangle }{\partial z} -\frac {2\langle u_z\rangle c'(a) \beta }{a} +\Bigl \langle u_z^{\prime}\frac {\partial c'}{\partial z}\Bigr \rangle \nonumber \\&\quad =D\left \lbrace \frac {2}{a}\frac {\partial c'}{\partial r}\bigg |_{r=a}+\frac {\partial ^2\langle c \rangle }{\partial z^2} -\frac {2\gamma c'(a)}{a}-\frac {4\beta }{a^2}\frac {\partial (rc')}{\partial z}\bigg |_{r=a}\right \rbrace \!. \end{align}
To find the fluctuating component of the previous equation, we can subtract the averaged component (2.14) from (2.8), yielding
\begin{align} &\frac {\partial c'}{\partial t}+\langle u_r\rangle \frac {\partial c'}{\partial r}-\frac {2\langle u_r \rangle }{a^2}\left (a c'(a)-\int _0^{a}c' \,{\rm d}r \right )+u_r^{\prime}\frac {\partial c'}{\partial r}-\Bigl \langle u_r^{\prime}\frac {\partial c'}{\partial r}\Bigr \rangle \nonumber \\&\quad +u_z^{\prime}\frac {\partial \langle c\rangle }{\partial z}+\langle u_z \rangle \frac {\partial c'}{\partial z}+\frac {2\langle u_z\rangle c'(a)\beta }{a}-\Bigl \langle u_z^{\prime}\frac {\partial c'}{\partial z}\Bigr \rangle +u_z^{\prime}\frac {\partial c'}{\partial z} \nonumber \\ &\quad =D\left \lbrace \frac {\partial ^2 c'}{\partial z^2}+\frac {1}{r}\frac {\partial }{\partial r} \left ( r\frac {\partial c'}{\partial r} \right )-\frac {2}{a}\frac {\partial c'}{\partial r}\bigg |_{r=a}+\frac {2\gamma c'(a)}{a}+\frac {4\beta }{a^2}\frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right \rbrace\! . \end{align}
2.2.2. Scale analysis for times
$\sim \sigma /U_0$
To perform a scale analysis, the equations are non-dimensionalised by using the following dimensionless variables with the superscript
$^*$
:
$\langle c \rangle =c_0 \langle c^*\rangle$
,
$c'=c^{\prime}_0 c'^*$
,
$z=\sigma z^*$
,
$r=a_0 r^*$
,
$a(z)=a_0 a^*(z^*)$
,
$\beta (z)=\beta _0 \beta ^*(z^*)$
,
$\gamma (z)=\gamma _0 \gamma ^*(z^*)$
,
$t=(\sigma /U_0)t^*$
,
$u^{\prime}_z=U_0u'^*_z$
,
$\langle u_z \rangle =U_0 \langle u^*_z\rangle$
,
$u^{\prime}_r=\beta _0 U_0u'^*_r$
,
$\langle u_r \rangle =\beta _0 U_0 \langle u^*_r\rangle$
and
$p=(\eta U_0 \sigma /a_0^2)p^*$
. In this context,
$ c_0$
and
$ c^{\prime}_0$
represent the characteristic solute concentrations for the averaged and fluctuating components, respectively, where
$ c^{\prime}_0 \ll c_0$
because
$ c^{\prime}_0$
is a small deviation in concentration from the cross-sectional average. Additionally,
$ \beta _0$
and
$ \gamma _0$
denote the first and second derivatives of
$ a$
at the position
$ z = 0$
.
When analysing this phenomenon, two characteristic time scales can be identified: one associated with molecular diffusion of the solute in the cross-section, which is approximately
$ a_0^2/D$
, and the other related to the advection of the solute, given by
$ \sim \sigma /U_0$
. Typically, we have
$ a_0^2/D \ll \sigma /U_0$
. In this study, we focus on the long time scale, where axial advection dominates over radial diffusion.
Dimensionless numbers will arise in the equations when substituting the dimensionless variables presented here. We identify two Péclet numbers:
$ {\textit{Pe}}_{a0} = U_0 a_0 / D$
and
$ {\textit{Pe}}_{\sigma 0} = U_0 \sigma / D$
, along with the Weisenberg number,
$ Wi = \lambda U_0 / a_0$
.
In the non-dimensionalisation of (2.15), certain terms can be disregarded based on the scale analysis provided in table 1. By focusing solely on the
$\mathcal{O}(1)$
and
$\mathcal{O}(1/{\textit{Pe}}_{a0})$
terms and converting back to the dimensional form, the Taylor dispersion equation retains the dominant advective dispersion term as well as the primary radial diffusion term, leading to
Summary of assumptions and parameters used for scaling analysis made in this study.
\begin{equation} u^{\prime}_z\frac {\partial \langle c\rangle }{\partial z}=D\left [\frac {1}{r}\frac {\partial }{\partial r} \left ( r\frac {\partial c'}{\partial r} \right )-\frac {2}{a}\frac {\partial c'}{\partial r}\bigg |_{r=a} \right ]\!. \end{equation}
2.2.3. Derivation of fluctuating concentration
Equation (2.16) is subject to the condition of impermeability of the wall
$\boldsymbol{\nabla }c \boldsymbol{\cdot }\boldsymbol{n}=0$
, where
\begin{equation} \boldsymbol{\nabla }c \equiv \left (\frac {1}{r}\frac {\partial }{\partial r}(rc), \frac {\partial c}{\partial x} \right )\!. \end{equation}
To obtain the normal vector to the surface, the gradient of the scalar function defines the channel walls, i.e.
$\boldsymbol{n}=\boldsymbol{\nabla }(r-a(z))= ( 1,-\beta )$
. Therefore,
$\boldsymbol{\nabla }c \boldsymbol{\cdot }\boldsymbol{n}=0$
on the surface is given by
\begin{equation} -\frac {\partial c}{\partial z}\bigg |_{r=a} \beta +\frac {1}{a}\frac {\partial }{\partial r}(ca)\bigg |_{r=a}=0, \end{equation}
and written in terms of
$\langle c \rangle$
and
$c'$
, we have
\begin{equation} -\frac {\partial \langle c\rangle }{\partial z}\beta -\frac {\partial c' }{\partial z}\bigg |_{r=a} \beta +\frac {\partial c'}{\partial r}\bigg |_{r=a} =0. \end{equation}
If we make the terms dimensionless, the term
$\partial c'/\partial z|_{r=a}$
can be ignored because it is
$\mathcal{O}(\epsilon ^2)$
. When this is combined with (2.16), the resulting expression is
\begin{equation} D\frac {1}{r}\frac {\partial }{\partial r} \left ( r\frac {\partial c'}{\partial r} \right )=u^{\prime}_z\frac {\partial \langle c\rangle }{\partial z}+\frac {2D\beta }{a}\frac {\partial \langle c\rangle }{\partial z} . \end{equation}
Integrating (2.20), using the expression for
$u^{\prime}_z$
given in (2.11), we have
\begin{equation} r\frac {\partial c'}{\partial r}=\frac {1}{D}\frac {\partial \langle c \rangle }{\partial z}\left [U_n\left (\frac {r^2}{2}-\frac {r^4}{2a^2} \right )+\frac {U_{{\textit {ve}}}}{2}\left ( \frac {r^2}{2}-\frac {r^6}{2a^4} \right ) \right ]+\frac {\beta }{a}\frac {\partial \langle c \rangle }{\partial z}r^2+C_1, \end{equation}
where
$C_1$
represents an integration constant. Considering the symmetry condition at
$r=0$
, we have that
$\partial c'/\partial r=0$
at that point, so
$C_1=0$
. Integrating again, it is possible to obtain
$c'$
:
\begin{equation} c'=\frac {1}{D}\frac {\partial \langle c \rangle }{\partial z}\left [U_n\left (\frac {r^2}{4}-\frac {r^4}{8a^2} \right )+\frac {U_{{\textit {ve}}}}{2}\left ( \frac {r^2}{4}-\frac {r^6}{12a^4} \right ) \right ]+\frac {\beta }{a}\frac {\partial \langle c \rangle }{\partial z}\frac {r^2}{2}+C_2. \end{equation}
The constant
$C_2$
can be found by considering that
$\langle \theta '\rangle =0$
, that is,
$\int _0^a c'r\,{\rm d}r=0$
, so the fluctuating concentration is given by
\begin{align} c'&=\frac {1}{D}\frac {\partial \langle c \rangle }{\partial z}\left [U_n\left (\frac {r^2}{4}-\frac {r^4}{8a^2}-\frac {a^2}{12} \right )+\frac {U_{{\textit {ve}}}}{2}\left ( \frac {r^2}{4}-\frac {r^6}{12a^4} -\frac {5a^2}{48}\right ) \right ] \nonumber \\&\quad +\frac {\beta }{a}\frac {\partial \langle c \rangle }{\partial z}\left (\frac {r^2}{2}-\frac {a^2}{4}\right )\!. \end{align}
2.2.4. Derivation of cross-section-averaged concentration
We revisit (2.14) and substitute (2.23) into all terms involving
$c'$
. After substituting and simplifying the terms, we arrive at
\begin{align} &\frac {\partial \langle c \rangle }{\partial t}+\frac {\partial \langle c \rangle }{\partial z}\left \lbrace \left ( 1-\frac {\beta ^2}{6}+\frac {\gamma a}{12}\right ) U_{n} +\left (1-\frac {\beta ^2}{8} +\frac {\gamma a}{16}\right ) U_{{\textit {ve}}} \right \rbrace \nonumber \\&\quad =D\left \lbrace \frac {\partial \langle c \rangle }{\partial z} \left ( \frac {2\beta }{a}-\frac {3\gamma \beta }{2}+\frac {3\beta ^3}{a}\right ) + \frac {\partial ^2\langle c \rangle }{\partial z^2}\left ( 1-\frac {a\beta }{12D}U_n+\frac {a^2}{48D^2}U^2_n \right .\right . \nonumber \\&\quad\quad \left .\left .-\beta ^2-\frac {a\beta }{16D}U_{{\textit {ve}}}+\frac {a^2}{96D^2}U^2_{{\textit {ve}}}+\frac {7a^2}{240D^2}U_nU_{{\textit {ve}}} \right )\right \rbrace \!. \end{align}
A scale analysis of (2.24) is conducted to identify terms that can be neglected, leading to a final formulation of the Taylor dispersion equation for this study. By using the characteristic scales outlined in § 2.2.2, we categorise the terms into different orders. According to table 1, the terms of
$\mathcal{O}(\beta _0)$
,
$\mathcal{O}(\epsilon )$
,
$\mathcal{O}(\epsilon \beta _0)$
and
$\mathcal{O}(\epsilon ^2)$
can be disregarded, leaving us with the terms
$\mathcal{O}(1)$
and
$\mathcal{O}(1/{\textit{Pe}}_{\sigma 0})$
. The resulting Taylor dispersion equation for a viscoelastic fluid flowing through arbitrarily shaped axisymmetric channels is expressed as follows:
\begin{align} \frac {\partial \langle c \rangle }{\partial t}+\frac {\partial \langle c \rangle }{\partial z}\!\left (\! U_{n} + U_{{\textit {ve}}}-\frac {2D\beta }{a} \!\right ) &=D \frac {\partial ^2\langle c \rangle }{\partial z^2}\!\left (\! 1+\frac {a^2}{48D^2}U^2_n +\frac {a^2}{96D^2}U^2_{{\textit {ve}}}+\frac {7a^2}{240D^2}U_nU_{{\textit {ve}}} \!\kern-0.5pt\right )\!\kern-1pt. \end{align}
Another way to express (2.25) in terms of
$U$
is
\begin{equation} \frac {\partial \langle c \rangle }{\partial t}+\frac {\partial \langle c \rangle }{\partial z}\left ( U-\frac {2D\beta }{a} \right ) =D \frac {\partial ^2\langle c \rangle }{\partial z^2}\left ( 1+\frac {a^2}{48D^2}U^2 -\frac {a^2}{96D^2}U^2_{{\textit {ve}}}-\frac {a^2}{80D^2}U_nU_{{\textit {ve}}} \right )\!. \end{equation}
From (2.25), we can identify the effective Taylor–Aris dispersion coefficient
$D_{{\textit {eff}}}/D$
, given by
Unlike the result obtained by Chang & Santiago (Reference Chang and Santiago2023), (2.27) includes two additional terms: one that accounts for the fluid’s viscoelasticity and another that represents the combined effects of Newtonian and viscoelastic behaviours. Equation (2.26) can also be expressed in a conservative form:
\begin{equation} \frac {\partial \langle c \rangle }{\partial t}+U\frac {\partial \langle c \rangle }{\partial z} =\frac {1}{a^2} \frac {\partial }{\partial z}\left ( \left ( D+\frac {a^2}{48D}U^2 -\frac {a^2}{96D}U^2_{{\textit {ve}}}-\frac {a^2}{80D}U_nU_{{\textit {ve}}} \right ) a^2\frac {\partial \langle c \rangle }{\partial z}\right )\!. \end{equation}
The classical result of Taylor–Aris dispersion in a fully developed laminar flow of a Newtonian fluid through a channel of constant radius can be recovered from (2.26), if
$\beta =0$
and
$U_{{\textit {ve}}}=0$
, i.e.
$Wi=0$
.
3. Solution methodology
When conducting studies on solute dispersion, it is desirable to understand the spatial and temporal distribution of the solute. However, it is not always necessary to completely solve the equation governing the phenomenon; the temporal distribution of the mean and variance can be analysed. Furthermore, other types of simulations, such as Brownian dynamics, can also be performed. In the following sections, the moments of solute distribution are analysed following the methodology developed by Chang & Santiago (Reference Chang and Santiago2023). An adaptation of Brownian distributions is performed based on their work. Finally, the complete solute distribution given in (2.26) is solved for certain limiting cases.
3.1. Time evolution of mean and variance
Following the methodology developed by Chang & Santiago (Reference Chang and Santiago2023), which is based on analysing the first and second moments for the mean and variance, at times sufficiently large compared with the times associated with radial diffusion, a mathematical model is formulated for the temporal distribution of the mean and variance. The average concentration, also known as the zeroth moment of solute distribution, is determined as
\begin{equation} \bar {c}=\int _{-\infty }^\infty a^2(z')\langle c \rangle (z',t) \,{\rm d}z'=\textrm {constant}, \end{equation}
with
$z'$
as a dummy variable of integration. The axial average of a variable is given as
\begin{equation} \bar {(\boldsymbol{\cdot })}=\frac {\int _{-\infty }^\infty (\boldsymbol{\cdot }) a^2(z')\langle c \rangle (z',t)\, {\rm d}z'}{\bar {c}}. \end{equation}
To derive the first moment or mean, we multiply (2.26) by
$z$
and apply the axial average described in (3.2),
\begin{align} &\int _{-\infty }^{\infty }za^2\frac {\partial \langle c \rangle }{\partial t}\,{\rm d}z+\int _{-\infty }^{\infty }za^2\frac {\partial \langle c \rangle }{\partial z}\left (U-\frac {2D\beta }{a} \right )\,{\rm d}z \nonumber \\&\quad =D \int _{-\infty }^{\infty }za^2\frac {\partial ^2\langle c \rangle }{\partial z^2}\left ( 1+\frac {a^2}{48D^2}U^2 -\frac {a^2}{96D^2}U^2_{{\textit {ve}}}-\frac {a^2}{80D^2}U_nU_{{\textit {ve}}} \right ){\rm d}z. \end{align}
Once the integrals have been evaluated and divided by
$\bar {c}$
, the dynamics of the axial mean position
$\bar {z}(t)$
is defined by the equation
\begin{equation} \frac {{\rm d}\bar {z}}{{\rm d}t}=U_0a_0^2\overline {\left (\frac {1}{a^2}\right )}+2D\overline {\left (\frac {\beta }{a}\right )}. \end{equation}
The temporal distribution of the mean position depends on the channel geometry, unlike the classic case of a constant-radius channel, where it depends only on the average velocity. Furthermore, viscoelasticity has no effect on the dynamics of the mean position. We consider a Taylor expansion of a function
$f(z)$
around the axial mean position
$\bar {z}$
to be able to perform an order analysis to simplify the model
\begin{equation} f(z)\approx f(\bar {z})+(z-\bar {z})f'(\bar {z})+\frac {1}{2}(z-\bar {z})^2f''(\bar {z})+\boldsymbol{\cdot }\boldsymbol{\cdot }\boldsymbol{\cdot }, \end{equation}
where
$f'$
and
$f''$
represent the first and second derivatives of
$f$
, respectively. Using the Taylor series around the mean position in (3.4) and performing an analysis of orders of magnitude, the variance
$\sigma _z^2=(z-\bar {z})^2$
is small with respect to the axial scale
$\sigma$
, that is,
$\sigma _z/\sigma \ll 1$
. Defining a new variable
$\widetilde {f}=f(\bar {z})$
simplifies the notation of the following equations:
\begin{equation} \frac {{\rm d}\bar {z}}{{\rm d}t}=U(\bar {z})+2D\widetilde {\left (\frac {\beta }{a}\right )} +O\left (U_0\left (\frac {\sigma _z}{\sigma }\right )^2,\frac {D}{\sigma } \left (\frac {\sigma _z}{\sigma }\right )^2\right )\!. \end{equation}
Similarly, to derive the second moment or variance of (2.26), we multiply (2.26) by
$(z-\bar {z})^2$
, apply the axial average integration and divide by
$\bar {c}$
, resulting in
\begin{align} \frac {{\rm d}\sigma _z}{{\rm d}t}&=2U_0a_0^2\left (\overline { \left (\frac {z}{a^2}\right )}-\bar {z}\overline {\left (\frac {1}{a^2} \right )}\right )+2D+4D\overline {\frac {\beta }{a}(z-\bar {z})}+\frac {U_0^2a_0^4}{24D} \overline {\left (\frac {1}{a^2} \right )}\nonumber \\&\quad -\frac {1}{48D^2}\overline {a^2U^2_{{\textit {ve}}}}-\frac {1}{40D^2}\overline {a^2U_nU_{{\textit {ve}}}}. \end{align}
Equation (3.7) can be simplified using a Taylor series around the axial mean position, where these series for the forms
$zf(z)$
and
$(z-\bar {z})f(z)$
are given as
\begin{align}& zf(z)\approx \bar {z}f(\bar {z})+(z-\bar {z})\left ( \bar {z}f'(\bar {z})+f(\bar {z})\right )+(z-\bar {z})^2 \left ( \frac {1}{2}\bar {z}f''(\bar {z})+f'(\bar {z})\right )+\cdot\cdot\cdot, \end{align}
\begin{align}&\qquad (z-\bar {z})f(z)\approx (z-\bar {z})f(\bar {z})+(z-\bar {z})^2f'(\bar {z})+\frac {1}{2}(z-\bar {z})^3f''(\bar {z})+\cdot\cdot\cdot. \\[9pt] \nonumber \end{align}
Using (3.5), (3.8) and (3.9) in (3.7), and applying an order of magnitude analysis, we obtain
\begin{align} \frac {{\rm d}\sigma _z^2}{{\rm d}t}&=2D+\frac {U_0^2a_0^4}{24D}\widetilde {\left (\frac {1}{a^2} \right )}-4U_0a_0^2\sigma _z^2\widetilde {\left (\frac {\beta }{a^3}\right )} -\frac {1}{48D}\widetilde {a^2U_{{\textit {ve}}}^2}-\frac {1}{40D}\widetilde {a^2U_nU_{{\textit {ve}}}}\nonumber \\&\quad +\textit{DO}\left ( \left (\frac {\sigma _z}{\sigma }\right )^2, {\textit{Pe}}_a^2 \left (\frac {\sigma _z}{\sigma }\right )^2, {\textit{Pe}}_a \frac {\sigma _z^3}{a_0\sigma ^2}Sk \right )\!, \end{align}
with
$Sk$
as the axial skewness of the solute, defined as
$Sk=\overline {(z-\bar {z})x^3}/\sigma _z^3$
. Unlike the mean axial position, the evolution of the variance defined in (3.10) does depend on the viscoelasticity of the fluid.
3.1.1. Boundary between the regimes of transient positive and negative growth of axial variance
Changes in geometry can produce changes in velocity and thus in advective dispersion, reflected in a variation in the axial variance of the solute. Chang & Santiago (Reference Chang and Santiago2023) explored the decrease of the axial variance of the solute in an expanding channel, where they explored the qualitative characterisation by positive
$\beta$
values,
${\textit{Pe}}_a$
values large enough for advective dispersion to dominate over molecular diffusion and a sufficiently large value of the variance with respect to the cross-section
$\sigma ^2_z/a^2$
. In this work, an analysis is carried out with the purpose of including the effect of the viscoelasticity of the fluid starting from (3.10). Therefore, it is necessary to define a new local Péclet number with two new components dependent on the Newtonian and viscoelastic velocity components, respectively,
\begin{equation} \widetilde {\textit{Pe}}_a=\frac {\widetilde{aU}\ }{D}=\widetilde {\textit{Pe}}_{a-n}+\widetilde {\textit{Pe}}_{a-{\textit {ve}}} \end{equation}
with
$\widetilde {\textit{Pe}}_{a-n}=\widetilde {aU_n}/D$
and
$\widetilde {\textit{Pe}}_{a-{\textit {ve}}}=\widetilde {aU_{{\textit {ve}}}}/D$
. Rearranging terms in (3.10), we have
\begin{equation} \frac {1}{D}\frac {{\rm d}\sigma _z^2}{{\rm d}t}=2\left (1+\frac {\widetilde {\textit{Pe}}^2_a}{48}-\frac {\widetilde {\textit{Pe}}^2_{a-{\textit {ve}}}}{96}-\frac {\widetilde {\textit{Pe}}_{a-n}\widetilde {\textit{Pe}}_{a-{\textit {ve}}}}{80} \right )-4\widetilde {\beta }\widetilde {\textit{Pe}}_a\left (\frac {\sigma _z}{\widetilde {a}} \right )^2. \end{equation}
Starting from (2.6) and (2.7), it is possible to write
$U$
only as a function of
$U_n$
, that is, write
$U_{ve}$
as a function of
$U_n$
, as follows:
\begin{equation} U=U_n+\frac {64\varepsilon \lambda ^2}{3a^2}U_n^3. \end{equation}
Likewise, it is also possible to express
$\widetilde {\textit{Pe}}_{a-{\textit {ve}}}$
as a function of
$\widetilde {\textit{Pe}}_{a-n}$
:
\begin{equation} \widetilde {\textit{Pe}}_{a-{\textit {ve}}}=\frac {64\varepsilon \widetilde {\textit{Wi}}^2}{3\widetilde {\textit{Pe}}^2_a}\widetilde {\textit{Pe}}^3_{a-n}, \end{equation}
where
$\widetilde {\textit{Wi}}=\lambda \widetilde {U}/\widetilde {a}$
is the local Weissenberg number. The first term on the right-hand side of (3.12) refers to the axial effective dispersion, due to molecular diffusion and the advective effect, which will always be positive, since despite the fact that the viscoelastic components are negative,
$\widetilde {{\textit{Pe}}_a}\gt \widetilde {\textit{Pe}}_{a-n},\widetilde {\textit{Pe}}_{a-{\textit {ve}}}$
always holds. However, the second term on the right-hand side is a factor due to geometry. For negative
$\widetilde {\beta }$
values, the growth of the variance will always be positive; however, for positive
$\widetilde {\beta }$
values, there may be limits where the growth is temporarily negative. If
$\widetilde {\beta }\gt 0$
and this term is large enough, it can counteract the variance due to effective dispersion, temporarily causing negative growth. The limit for
$\widetilde {\beta }\gt 0$
is analysed later, where these two effects counteract each other and zero growth occurs, that is,
${\rm d}\sigma ^2_z/{\rm d}t=0$
, being the border between positive and negative growth, as follows:
\begin{equation} \left (\frac {\sigma _z}{\widetilde {a}} \right )^2=\frac {1}{4\widetilde {\beta }}\left (\frac {2}{ \widetilde {\textit{Pe}}_a}+\frac {\widetilde {\textit{Pe}}_a}{24}-\frac {\widetilde {\textit{Pe}}^2_{a-{\textit {ve}}}}{48 \widetilde {\textit{Pe}}_a}-\frac {\widetilde {\textit{Pe}}_{a-n}\widetilde {\textit{Pe}}_{a-{\textit {ve}}}}{40 \widetilde {\textit{Pe}}_a} \right )\!. \end{equation}
The expression given in (3.15) shows a relationship between the ratio of the variance and the local channel radius
$\sigma ^2_z/\widetilde {a}^2$
, the local Péclet number
$\widetilde {\textit{Pe}}_a$
, the local slope of the radius
$\widetilde {\beta }$
and the local Weissenberg number
$\widetilde {\textit{Wi}}$
implicitly in
$\widetilde {Pa}_{a-{\textit {ve}}}$
for the zero-growth boundary. This is key to designing channels that control solute dispersion. It is useful to find the Péclet number at which the minimum variance occurs. To do this, we differentiate (3.15) with respect to
$\widetilde {\textit{Pe}}_a$
and find the minimum value. From (3.11) and (3.14), we find that
${\rm d}\widetilde {\textit{Pe}\ }_{a-n}/{\rm d}\widetilde {\textit{Pe}}_{a}=\widetilde {\textit{Pe}}_{a-n}/\widetilde {\textit{Pe}}_{a}$
and
${\rm d}\widetilde {\textit{Pe}}_{a-{\textit {ve}}}/{\rm d}\widetilde {\textit{Pe}}_{a}=\widetilde {\textit{Pe}}_{a-{\textit {ve}}}/\widetilde {\textit{Pe}}_{a}$
. The expression for the Péclet number when the minimum variance occurs is given by
\begin{equation} \widetilde {\textit{Pe}}_a=\sqrt {48+\frac {1}{2}\widetilde {\textit{Pe}}^2_{a-{\textit {ve}}}+\frac {3}{5}\widetilde {\textit{Pe}}_{a-n}\widetilde {\textit{Pe}}_{a-{\textit {ve}}}}. \end{equation}
Note that the value at which the minimum variance occurs is not a fixed value, unlike the Newtonian case which corresponds to
$\sqrt {48}$
, but rather depends on the viscoelasticity. To find the relationship between
$\widetilde {\textit{Pe}}_a$
and
$\widetilde {\textit{Wi}}$
, (3.16) must be solved together with (3.11) and (3.14).
3.1.2. Approximations useful for channel design
Following the methodology used by Chang & Santiago (Reference Chang and Santiago2023), it is possible to further simplify (3.12) by implementing the approximation
\begin{equation} \frac {\rm d}{{\rm d}t}\approx \frac {{\rm d}\bar {z}}{{\rm d}t}\frac {\rm d}{{\rm d}z}. \end{equation}
Using this approximation in (3.10), writing
$\beta$
as its definition, using (3.6) and (3.13), we obtain an ODE (ordinary differential equation) for
$a$
for a desired
$\sigma _z^2$
distribution:
\begin{equation} \frac {{\rm d}a}{{\rm d}z}\approx \frac {-\dfrac {{\rm d}\sigma _z^2}{{\rm d}z}+\dfrac {U_0a_0^2}{24D}+\dfrac {2D}{U_0a_0^2}a^2-\dfrac {256\varepsilon ^2\lambda ^4}{27U_0a_0^2D}U_n^6-a^2\dfrac {8\varepsilon \lambda ^2}{15U_0a_0^2D}U_n^4}{\left (\dfrac {2D}{U_0a_0^2}\dfrac {{\rm d}\sigma _z^2}{{\rm d}z} \right )a+\dfrac {4\sigma _z^2}{a}}, \end{equation}
where
$U_n$
, according to (3.13), can be obtained as the root of
\begin{equation} a^2U_n+\frac {64\varepsilon \lambda ^2}{3}U_n^3=U_0a_0^2. \end{equation}
It is possible to design channels that can result in a specific pattern of variance evolution in space when solved numerically by (3.18).
3.2. Brownian dynamics simulations
To validate our expressions on the evolution of the variance and the position of the average, we used Brownian dynamics simulations. To do this, we used an adaptation of the code developed by Chang and available on Github, modifying the velocity profiles given in (2.4) and (2.5) and the differential equation of the variance given in (3.10), and the data output for analysis and comparison with respect to the simulation available online here. Each simulation consisted of tracking 5000 point particles, with an initial axial variance of
$\sigma ^2_{z0}=10a_0^2$
. It is developed based on the direct Euler method. For each case, five simulations were performed using different random seed initial conditions, and the reported values correspond to their average. The numerical evaluation of the axially averaged quantities
$\overline {f(x)}$
is given by
$\overline {f(x)} = \varSigma _{i=1}^{N_{{\textit {particle}}}}f(x_i)/N_{{\textit {particle}}}$
, where
$x_i$
denotes the axial position of the ith particle.
3.3. Spatial and temporal solute distribution
The equations required to determine
$ \langle c \rangle$
in their dimensionless form are derived from (2.6), (2.7), (2.25) and (2.27):
\begin{align}& U_{n}^*(z)=\frac {a^{*2}}{8}\left (-\frac {{\rm d}p^*}{{\rm d}z^*} \right )\!, \quad U_{ve}^*(z)= \frac {\varepsilon {\textit{Wi}}^2 a^{*4}}{24}\left (-\frac {{\rm d}p^*}{{\rm d}z^*} \right )^3, \end{align}
\begin{align}& \frac {\partial \langle c \rangle ^*}{\partial t^*}+\frac {\partial \langle c \rangle ^*}{\partial z^*}\left ( U_{n}^* + U_{{\textit {ve}}}^*-\frac {1}{{\textit{Pe}}_{\sigma 0}} \frac {2\beta ^* }{a^*} \right ) \nonumber \\&\quad = \frac {1}{{\textit{Pe}}_{\sigma 0}}\frac {\partial ^2\langle c \rangle ^*}{\partial z^{*2}}\left ( 1+{\textit{Pe}}_{a0}^2a^{*2}\left (\frac {U^{*2}_n}{48} +\frac {U^{*2}_{{\textit {ve}}}}{96}+\frac {7U_n^*U_{{\textit {ve}}}^*}{240} \right )\right ) \end{align}
and
\begin{equation} D^*_{{\textit {eff}}}= 1+{\textit{Pe}}_a^2a^{*2} \left (\frac {U^{*2}}{48} -\frac {U^{*2}_{{\textit {ve}}}}{96}-\frac {U_n^*U_{{\textit {ve}}}^*}{80} \right )\!. \end{equation}
Equation (3.21) is subject to the boundary and the initial conditions:
\begin{equation} \langle c \rangle ^*(z^*\to \pm \infty ,t^*)=0 \end{equation}
and
\begin{equation} \langle c \rangle ^*(z^*,0)=\delta (z), \end{equation}
where
$\delta (z)$
represents the Dirac delta function. While it is acknowledged that the assumption of a cross-sectionally uniform concentration profile, inherent in the derivation of (3.21), is not valid immediately after injection (
$t^*=0$
) for this initial condition (Chatwin Reference Chatwin1970), the fundamental solution to this equation provides an exact description of long-term dispersion. This is valid once transverse diffusion has had time to homogenise the concentration profile, a process that occurs over a short time scale
$t_c \sim a^2 / D$
. For times
$t \gg a^2 / D$
, which are the focus of this study, the solution to (3.21) accurately captures the spatio-temporal evolution of the cross-sectional average solute concentration. In this section, we will derive an analytical solution to (3.21) for a channel with a constant radius. We will employ a change of variable that simplifies the model, allowing it to be solved numerically for any arbitrary geometry.
3.3.1. Constant radius channel case
The simplest geometry that the channel can have is a uniform cross-section. To simplify the problem, we introduce a new variable,
$\xi (z, t) = z^* - t^*$
, which represents a spatial variable moving at the fluid’s velocity. Additionally, we define
$\tau (t^*) = t^*$
. Therefore, we can express (3.21) in terms of these new variables as follows:
\begin{equation} \frac {\partial \langle c \rangle ^*}{\partial \tau } = D_{{\textit {eff}}-c}\frac {\partial ^2\langle c \rangle ^*}{\partial \xi ^{2}}, \end{equation}
with
$D^*_{{\textit {eff}}-c}=D^*_{{\textit {eff}}}/{\textit{Pe}}_{\sigma 0}$
, and subject to
\begin{equation} \langle c \rangle ^*(\xi \to \pm \infty ,\tau )=0, \quad \langle c \rangle ^*(\xi ,0)=\delta (\xi ). \end{equation}
Using the Fourier transform, one can derive a well-known analytical solution for (3.25) as follows:
\begin{equation} \langle c\rangle ^*=\frac {1}{\sqrt {4\pi D^*_{{\textit {eff}}-c}\tau }}\exp {\left (-\frac {\xi ^2}{4D_{{\textit {eff}}-c}\tau }\right )}. \end{equation}
The solution (3.27) corresponds to a Gaussian profile, analogous to that found by Taylor (Reference Taylor1953) in his pioneering work. He established the validity of this solution for
${\textit{Pe}}_a \gt 6.9$
based on experimental results.
3.3.2. Simplifying the mathematical model
In § 3.3.1, a change of variable was introduced to simplify the mathematical model, enabling the possibility of finding an analytical solution to the (3.21). However, since the channel radius was constant, the velocities – and, consequently, the effective diffusion coefficient – remain unchanged in the flow direction. This assumption is not valid for arbitrary geometries. Motivated by this limitation, we seek variables that will enable the mathematical model to be simplified within a general framework applicable to any function
$a(z)$
.
A change of variables is sought that allows us to simplify (3.21). To do this, two variables are proposed,
$\tau (t)$
and
$\zeta (z,t)$
, whose form is not yet known, but writing (3.21) in terms of these two unknown variables yields
\begin{equation} \frac {{\rm d} \tau }{{\rm d} t^*} \frac {\partial \langle c \rangle ^*}{\partial \tau } +\frac {\partial \zeta }{\partial t^*} \frac {\partial \langle c \rangle ^*}{\partial \zeta }+\frac {\partial \zeta }{\partial z^*} \frac {\partial \langle c \rangle ^*}{\partial \zeta }\left ( U^*_{n} + U^*_{{\textit {ve}}}-\frac {2\beta ^* }{{\textit{Pe}}_\sigma a^*} \right ) =\frac {D^*_{{\textit {eff}}}}{{\textit{Pe}}_\sigma } \left (\frac {\partial \zeta }{\partial z^*} \right )^2\frac {\partial ^2 \langle c \rangle ^*}{\partial \zeta ^2 }. \end{equation}
We conduct a search for the variable
$\zeta$
in a way that simplifies (3.28) for a numerical solution by eliminating the second and third terms of (3.28). The standard form of a parabolic partial differential equation is expressed as
$c_t = Dc_{zz}$
. To transform (3.28) into this standard form, we will analyse the second and third terms to identify the variable
$ \zeta$
that allows these two terms to vanish. Consequently, we will derive a differential equation for
$ \zeta$
in the following form:
\begin{equation} \frac {\partial \zeta }{\partial t^*} +\frac {\partial \zeta }{\partial z^*} \left ( U^*_{n} + U^*_{ve}-\frac {2\beta ^* }{{\textit{Pe}}_{\sigma 0} a^*} \right )=0. \end{equation}
A solution for (3.29), assuming a constant flow rate,
$(U_n^*(z)+U^*_{{\textit {ve}}}(z^*))a^{*2}=1$
, is given by
\begin{equation} \zeta =\int \frac { a^{*2}}{1-2\beta ^* a^*/{\textit{Pe}}_{\sigma 0}}\,{\rm d}z^* -t^* . \end{equation}
Considering in the previous equation that
${\textit{Pe}}_{a0}\gt 6.9$
and that
${\textit{Pe}}_{\sigma 0}=U_0\sigma /D= (U_0 a_0/D ) ( \sigma /a_0 )={\textit{Pe}}_{a0}/\epsilon$
, and since
$\epsilon \ll 1$
, therefore,
${\textit{Pe}}_\sigma \gg 1$
, so
$2\beta ^* a^*/{\textit{Pe}}_{\sigma 0} \ll 1$
, then the proposed variable can be written, in a first approximation, as
\begin{equation} \zeta =\int a^{*2} \,{\rm d}z^* -t^* . \end{equation}
The variable
$\zeta$
is derived mathematically, but can also be interpreted physically, similarly to the discussion in § 3.3.1. In this context,
$\zeta$
is a new variable that moves with the fluid velocity and is adjusted to counteract the effects of changes in the channel. Specifically, in (3.31) for the constant ratio case
$(a^*=1$
,
$\beta ^*=0$
), the argument of the function aligns with the definition presented in § 3.3.1, where
$\zeta = \xi = z^* - t^*$
.
Considering
$\tau =t$
and the proposal derived in (3.31) and (3.21), we obtain
\begin{equation} \frac {\partial \langle c \rangle ^*}{\partial \tau } =D^*_{{\textit {eff}}-a} \frac {\partial ^2 \langle c \rangle ^*}{\partial \zeta ^2 }, \end{equation}
where
\begin{equation} D^*_{{\textit {eff}}-a}=\frac {D^*_{{\textit {eff}}}}{{\textit{Pe}}_{\sigma 0}} a^{*4}, \end{equation}
subject to
\begin{equation} \langle c \rangle ^*(\zeta \to \pm \infty ,\tau )=0, \quad \langle c \rangle ^*(\xi ,0)=\delta (\zeta ). \end{equation}
The effective diffusion coefficient,
$D^*_{{\textit {eff}}}$
, defined in (3.22), and the channel radio
$a^*$
depend on their spatial position,
$z^*$
. To solve (3.32), we must express the entire equation in terms of the variables
$\zeta$
and
$\tau$
. This means that we need to convert
$D^*_{{\textit {eff}}}$
and
$a^*$
into functions of these variables:
$D^*_{{\textit {eff}}}(z) \to D^*_{{\textit {eff}}}(\zeta , \tau )$
and
$a^*(z^*) \to a^*(\zeta , \tau )$
. The relationship between
$z^*$
and the variables
$\zeta$
and
$\tau$
is derived from (3.31) using numerical methods, specifically the Newton–Raphson method.
Mathematically, (3.32) is an equation with variable coefficients, which requires the use of numerical methods for its solution. The Crank–Nicolson method is employed here as it provides numerical stability and second-order accuracy. This method is obtained by averaging the implicit and explicit approximations of the diffusion equation. The implementation of this method is done using MATLAB software, and the results will be presented in the following section.
In certain applications, understanding the temporal position
$ \bar {z}^*$
of the mean concentration (which corresponds to the peak concentration for symmetrical distributions) is often more relevant than knowing the entire solute distribution. The solutions to an equation such as (3.32) tend to develop around
$ \zeta = 0$
in a nearly symmetrical manner (with zero skewness) due to the conditions of this study that involve slow variations in the channel shape. However, this does not imply that the concentration peak remains stationary in real space
$( z^* )$
; therefore, to analyse the movement of the solute, it is essential to refer back to the original independent variables. By assuming that the mean concentration is positioned at
$\zeta = 0$
, we can determine the spatial and temporal evolution of the mean concentration by solving the following equation, based on (3.31):
\begin{equation} \int a^*(\bar {z}^*)^2\,{\rm d}\bar {z}^* -t^*=0. \end{equation}
4. Results and discussion
Using classical lubrication theory, the velocity distribution of a non-Newtonian fluid through an axisymmetric channel of arbitrary shape given in (2.4) and (2.5) was found. These velocities depend on the radial (
$r$
) and axial (
$z$
) positions. They are written in terms of two average velocities,
$U_n$
and
$U_{ve}$
, which correspond to the Newtonian and viscoelastic components of the velocity averaged over the cross-section. Subsequently, starting from the advection–diffusion equation, a Taylor analysis was performed, allowing us to obtain (3.21), which describes the dynamics of the average concentration in an arbitrarily shaped channel through which a viscoelastic fluid flows. Using (3.21), the first moment for the average position and the second moment for the variance were analysed, allowing to obtain a system of ODEs that describe them, in addition to making a proposal to numerically solve the temporal and spatial distribution of the average concentration under the condition of
${\textit{Pe}}_{\sigma 0}\gg 1$
. Brownian simulations were performed to validate the results of this work. The dimensionless parameter that takes relevance in hydrodynamics is the Weissenberg number, involved in
$U_{{\textit {ve}}}$
. However, for the concentration distribution, the dimensionless parameter that becomes relevant is the Péclet number, calculated from both
$a_0$
(
${\textit{Pe}}_{a0}$
) and
$\sigma$
(
${\textit{Pe}}_{\sigma 0}$
). Finally, another aspect that governs the entire phenomenon is precisely the channel geometry. We consider the channel geometry to be given in the form
$a/a_0=1+f(z)$
, where
$f(z)$
is an arbitrary function that describes the variation (which is small, of order
$\beta _0$
) of the channel’s radius along the channel.
To estimate the dimensionless parameters used for the calculations, a suitable combination of values for the physical parameters reported in the literature was selected (Mukherjee et al. Reference Mukherjee, Dhar, DasGupta and Chakraborty2019). Towards that, the corresponding estimates of various dimensional parameters have been taken as
$a_0\sim 10^{-5}$
m,
$U_0\sim 10^{-3}$
m s
$^{-1}$
,
$D\sim 10^{-9}$
m
$^2$
s
$^{-1}$
,
$\sigma \sim 10^{-2}$
m,
$\beta _0 \sim 10^{-3}$
m,
$\lambda \sim 10^{-3}$
s and the extensibility parameter of the sPTT model is assumed as
$0\leq \varepsilon \leq 1$
.
Velocity profile for different values of the parameter
$\varepsilon {\textit{Wi}}^2$
, considering
${\textit{Pe}}_{a0}=10$
in a divergent linear channel with
$\beta =0.01$
at (a)
$z/a_0=0$
and
$(b)$
$z/a_0=10$
.
4.1. Fluid velocity
The velocity profiles, as described by (2.4), (2.5) and (3.20), are graphically illustrated in figure 2. These profiles facilitate an analysis of how fluid viscoelasticity and channel geometry influence the hydrodynamic behaviour of the flow.
In the case of a Newtonian fluid (where
$\varepsilon {\textit{Wi}}^2=0$
), the velocity profile displays the classic parabolic shape typical of Poiseuille flow. However, as the viscoelastic parameter
$\varepsilon {\textit{Wi}}^2$
increases, the profile becomes noticeably flatter. This change occurs due to the introduction of elastic stresses that redistribute the shear stress, thereby decreasing the apparent viscosity through the ‘shear thinning’ effect described by the function
$f(\tau _{\textit{kk}})$
in the PTT model, as noted by Oliveira & Pinho (Reference Oliveira and Pinho1999). At higher values of
$\varepsilon {\textit{Wi}}^2$
, the velocity profile increasingly resembles a plug-like flow. Given that the flow rate remains constant compared with the Newtonian case, the velocity at the centre of the channel diminishes as
$\varepsilon {\textit{Wi}}^2$
increases. This behaviour underscores the significant impact of viscoelasticity on flow dynamics when compared with purely viscous (Newtonian) fluids.
The effect of channel geometry is illustrated in figure 2. Both figures 2(a) and 2(b) display the velocity profile of the fluid flowing through a diverging linear channel described by the equation
$ a/a_0 = 1 + 0.01z/a_0$
, at positions
$ z/a_0 = 0$
and
$ z/a_0 = 10$
, respectively. The parameters used in these graphs are
$ {\textit{Pe}}_{a0} = 10$
,
$ \beta = 0.01$
and various values of
$ \varepsilon {\textit{Wi}}^2$
. At
$ z = 1$
, where the channel has a larger radius, the fluid velocities are lower compared with those at the beginning of the channel at
$ z = 10$
, which can be attributed to the constant flow rate. Thus, it is evident from the graphs that the radius size directly influences the velocity profile, resulting in lower deformation rates and stresses within the fluid.
4.2. Effective dispersion coefficient
The effective dispersion coefficient
$D_{{\textit {eff}}}$
quantifies how easily a solute disperses in a medium, integrating influences from both advective transport and molecular diffusion. As derived in (3.22),
$D_{{\textit {eff}}}$
is shaped by three crucial factors: the Péclet number
${\textit{Pe}}_{a0}$
, the term
$\varepsilon {\textit{Wi}}^2$
(which is implicitly included in
$U_{ve}$
) and the geometry characterised by
$a$
. Figure 3
$(a)$
illustrates the impact of viscoelasticity on
$D_{{\textit {eff}}}/D$
. It is evident that as the parameter
$\varepsilon {\textit{Wi}}^2$
increases, the value of
$D_{{\textit {eff}}}/D$
decreases throughout the channel. According to (2.27), the presence of a viscoelastic component reduces
$D_{{\textit {eff}}}/D$
compared with the Newtonian case. The presence of viscoelasticity dampens the hydrodynamic dispersion due to a flattening of the velocity profile, as shown in figure 2. Specifically, viscoelasticity reduces velocity gradients near the centre of the channel, resulting in decreased radial diffusion and, consequently, lower dispersion rates. This observation aligns with findings from similar studies conducted by Akyildiz (Reference Akyildiz2002).
Figure 3
$(b)$
illustrates the impact of channel geometry on the effective dispersion coefficient,
$D_{{\textit {eff}}}/D$
. This behaviour is closely related to the channel shape, as previously noted by Chang & Santiago (Reference Chang and Santiago2023). For sinusoidal functions represented by the solid black and dashed blue lines,
$D_{{\textit {eff}}}$
reaches its lowest value when the channel is at its widest point and vice versa. In the case of a diverging channel, represented by the dash-dotted red line, we observe that
$D_{{\textit {eff}}}$
decreases as the channel widens. Conversely, the green dotted line represents a converging channel, showing that
$D_{{\textit {eff}}}$
increases as the channel narrows. This behaviour can be attributed to the constant flow rate condition: as the cross-sectional area increases, the fluid velocity tends to decrease, leading to a reduction in
$D_{{\textit {eff}}}$
.
Results of the effective diffusion coefficient considering
${\textit{Pe}}_{a0}=10$
$(a)$
for a diverging linear channel with
$\beta =0.01$
,
$\epsilon =0.001$
and various values of
$\varepsilon {\textit{Wi}}^2$
and
$(b)$
for different geometries, considering
$\varepsilon {\textit{Wi}}^2=10$
.
4.3. Boundary between the regimes of transient positive and negative growth of axial variance
According to (3.15), there is a given combination of the local Péclet number values
$\widetilde {\textit{Pe}}_a$
and local axial variation of the channel radius
$\tilde {\beta }$
that produce zero transient growth. Viscoelasticity is also found implicitly by
$\varepsilon \widetilde {\textit{Wi}}^2$
, which is within
$\widetilde {\textit{Pe}}_{a-{\textit {ve}}}$
in (3.15). This happens because the term related to the effective axial dispersion is counteracted by the variance produced by the channel geometry in (3.12). It is worth noting that this balance between these two terms only occurs for channels with a positive
$\beta$
value; otherwise, the transient growth of the variance will always be positive. Physically, a divergent channel could transiently have a negative growth due to the widening of the channel, but there is a point where this negative growth is counteracted by the effective dispersion (molecular and advective), and that is the case studied here. Figure 4 shows the regimes of zero transient growth of the axial variance obtained from (3.15) expressed by means of the variance growth
$(\ln (\sigma ^2_z/a^2))$
for different values of the local Péclet number
${\textit{Pe}}_a$
and the axial local variation of the radius
$\beta$
. Here, the tilde is not used to represent local values, since the difference is made by the subscript 0 for the values that are calculated from the reference values (
$U_0$
and
$a_0$
).
Zero transient growth regimes of the solute axial variance obtained from (3.15), expressed by means of the variance growth surface
$(\ln (\sigma ^2_z/a^2))$
for different values of the Péclet number
${\textit{Pe}}_a$
and the axial variation of the radius
$\beta$
for the cases of
$(a)$
a Newtonian fluid (i.e.
$\varepsilon {\textit{Wi}}^2=0$
) and
$(b)$
a viscoelastic fluid with
$\varepsilon {\textit{Wi}}^2=10$
. Contour lines are shown corresponding to different values of
$\ln (\sigma ^2_z/a^2)$
indicated on each contour line, for the cases of
$(c)$
a Newtonian fluid and
$(d)$
a viscoelastic fluid with
$\varepsilon {\textit{Wi}}^2=10$
.
Figure 4
$(a,\!b)$
show the zero transient growth surface, for a Newtonian case with
$\varepsilon {\textit{Wi}}^2=0$
in panel (a) and for a viscoelastic fluid with
$\varepsilon {\textit{Wi}}^2=10$
in panel (b), being that they follow a graphically similar behaviour, but with different values that can be better appreciated with contour lines. The space above the surface represents a negative transitory growth of the variance since it would mean that the term related to the geometric effect in (3.12) is the dominant one and, in contrast, the space below the surface represents a positive transitory growth of the variance.
However, figure 4
$(c)$
shows the contour lines for different values of
$(\ln (\sigma ^2_z/a^2))$
for the case of a Newtonian fluid, while figure 4
$(d)$
shows it for the case of a viscoelastic fluid with
$\varepsilon {\textit{Wi}}^2=10$
. Unlike figure 4
$(a,\!b)$
, the viscoelastic effect is evident, with these curves having different behaviour. An important point that can be observed in these figures is that to have a zero variance transient growth, for the case of viscoelastic fluid, smaller axial variations of the channel radius are required (lower values of
$\beta$
) compared with a Newtonian fluid, for a given Péclet number. Furthermore, according to what was found in (3.16), the Péclet number at which the smallest variance occurs does not correspond to the fixed value of
$\sqrt {48}$
at which it occurs in the Newtonian case, but rather occurs at a larger number that depends on the viscoelasticity of the fluid.
4.4. Axial variance analysis in different geometries
We analyse how the variance behaves across channels with different geometry, comparing the solutions obtained from the system of ODEs described in (3.6) and (3.10), with the corresponding Brownian simulations. The viscoelastic effect on the distribution of variance is compared with the effect that Deff has on the variance, and finally the validity of the model for large Péclet numbers is analysed.
Figure 5 shows the axial variance
$\sigma _z^2/a^2$
and the mean position
$\bar {z}/a_0$
of the solute along the linear channels, comparing different values of
$\varepsilon {\textit{Wi}}^2$
. The solid lines were determined from (3.6) and (3.10), while the dotted lines were determined from Brownian simulations. Figure 5
$(a)$
corresponds to a diverging channel whose radius is given by
$a/a_0=1+0.001z/a_0$
and figure 5
$(b)$
corresponds to a convergent channel whose radius is given by
$a/a_0=1-0.001z/a_0$
. In the case of the convergent channel, we can observe how the variance grows in greater proportion with respect to its mean position, compared with the divergent channel, which does so more slowly. In both cases, viscoelasticity reduces the variance compared with the Newtonian case, and this corresponds to the results shown in figure 3
$(a)$
, where it is observed that viscoelasticity reduces the value of
$D_{{\textit {eff}}}/D$
. The Brownian simulations and our model agree very well, which allows us to have confidence in the variance distribution found for the cases analysed.
Axial variance of the solute along the channel determined from (3.6) and (3.10) represented by the solid lines, for different values of
$\varepsilon {\textit{Wi}}^2$
, and with dotted lines for the Brownian simulation, with
$\epsilon =0.0001$
.
$(a)$
Diverging channel whose radius is given by
$a/a_0=1+0.001z/a_0$
and
$(b)$
convergent channel whose radius is given by
$a/a_0=1-0.001z/a_0$
, both with
${\textit{Pe}}_{a0}=10$
.
Axial variance of the solute along a sinusoidal periodic channel whose radius is given by
$a/a_0=1+0.2\sin {z/400a_0}$
determined from (3.6) and (3.10), with
$\varepsilon {\textit{Wi}}^2=10$
and
$\epsilon =0.001$
.
$(a)$
Axial variance of the solute (left vertical axis) represented by the black solid line and with dotted lines for the Brownian simulation with
${\textit{Pe}}_{a0}=10$
, and the effective dispersion coefficient
$D_{{\textit {eff}}}/D$
from (3.22) (left vertical axis) along the same channel represented by the blue solid line.
$(b)$
Axial variance represented with solid lines for different Péclet numbers and with dotted lines for its value from the Brownian simulation.
Figure 6 shows the axial variance
$\sigma _z^2/a^2$
and the mean position
$\bar {z}/a_0$
of the solute along a sinusoidal channel whose radius is given by
$a/a_0=1+0.2\sin {z/400a_0}$
, comparing the evolution of the variance with
$D_{{\textit {eff}}}/D$
, in addition to analysing the agreement of the studied model and the Brownian simulations (dotted lines) for different Péclet numbers. Figure 6
$(a)$
shows that the shape of
$D_{{\textit {eff}}}/D$
is not sufficient to fully predict how the variance will evolve, since although it has very similar distribution shapes, in different regions, they behave very differently, due to the effects of the channel geometry. Mathematically, this difference can be observed from (3.12), where writing the variance as a function of
$D_{{\textit {eff}}}$
gives
${\rm d}\sigma _z^2/{\rm d}t=2D_{{\textit {eff}}}-4D\widetilde {\beta }\widetilde {\textit{Pe}}_a (\sigma _z/\widetilde {a} )^2$
, where the second term on the right-hand side marks a difference between
$D_{{\textit {eff}}}$
and
$\sigma ^2_z$
, which does not exist for channels with constant shape. That is why for
$\beta =0$
, it is enough to analyse
$D_{{\textit {eff}}}$
to understand how
$\sigma _z^2$
behaves. Also with this geometry, the Brownian simulations have very good agreement with the studied model. Figure 6
$(b)$
shows the variance that develops in the same channel, but for different
${\textit{Pe}}_{a0}$
, where for values of
${\textit{Pe}}_{a0}=10$
and
${\textit{Pe}}_{a0}=100$
, we find a good agreement between this model and the Brownian simulations, but for the case of
${\textit{Pe}}_{a0}=1000$
, there are differences that may be important.
4.5. Engineering channels design for given axial variance patterns
Using the analytical approximation given in (3.18), it is possible to obtain a numerical solution to design a geometry that produces a desired variance distribution pattern. Figure 7
$(a,\!b)$
shows two cases where, starting from a target variance distribution, the geometry that produces it is obtained solving numerically (3.18) and, in turn, from it, the variance is calculated based on (3.6) and (3.10) and also from Browian simulations, with
$\varepsilon {\textit{Wi}}^2=10$
and
${\textit{Pe}}_{a0}=10$
.
In figure 7
$(a)$
, a constant variance distribution is the target with a value of
$\sigma _z^2/a^2=300$
. We can observe, with black solid lines, the geometry that results from solving numerically (3.18) for this case, in addition to three other lines that correspond to the target variance, the variance calculated from this geometry found and the variance calculated with Brownian simulations also for the geometry found. Additionally, figure 7
$(b)$
shows a similar analysis, but for a target variance given by
$\sigma _z^2/a^2=300+50\sin {2\pi z/600a_0}$
. In both cases, we can observe a fairly acceptable match between the target, calculated and Brownian variances.
For scientific and engineering applications involving microfluidics, the ability to design a channel that achieves the desired solute variance for solute dispersion in viscoelastic fluids can represent a significant advance.
Engineering the evolution of the axial variance from (3.18). Two-channel designs are shown, where the channel radius variation is shown with black solid lines (left vertical axis) and the variance distribution calculated from (3.18) is shown with blue solid lines. The variance determined from Brownian simulations is shown with red dotted lines and the design target variance is shown with green dashed lines (right vertical axis). The design for panel
$(a)$
consists of obtaining a constant variance
$\sigma _z^2/a^2=300$
and for panel
$(b)$
, one to maintain a sinusoidally distributed variance given by
$\sigma _z^2/a^2=300+50\sin {2\pi z/600a_0}$
, both with
$\varepsilon {\textit{Wi}}^2=10$
and
${\textit{Pe}}_{a0}=10$
.
A comparison made between the analytical solution for the case of a constant radius, as given in (3.27), and the numerical solution represented by (3.32). This comparison considers the parameters
$ a/a_0 = 1$
,
$ {\textit{Pe}}_{a0} = 10$
,
$\epsilon =0.001$
and
$ \varepsilon {\textit{Wi}}^2 = 10$
, evaluated at time points
$ t/(a_0^2/2D) = 2$
,
$ 20$
,
$ 40$
and
$ 60$
. Both solutions are derived from our work.
4.6. Average concentration distribution
The dispersion dynamics in channels reveals fundamentally different behaviours between constant-radius and axially varying geometries. In a uniform cylindrical channel, the variance of the solute increases linearly with time, following the relationship
$2D_{\textit{eff}} t$
as described by classical Taylor dispersion theory. However, this straightforward relationship does not hold in channels with axial geometric variations. To better understand how the average concentration changes over time, one can analyse the evolution of the mean concentration and its variance, as we did in the previous section, or by finding a numerical solution, as is done in this one. The equation that is solved numerically is the one proposed in (3.32), by using the variables
$\tau$
and
$\zeta$
proposed for the limit
${\textit{Pe}}_{a0}/\epsilon \gg 1$
. To validate our model, it is solved for the case of a constant radius channel,
$\beta =0$
, and compared with the analytical solution given in (3.27), which is a well-known classical concentration distribution solution. This comparison can be observed in figure 8, showing a perfect agreement between both solutions. Furthermore, the results for different geometries will be compared qualitatively with the variance
$\sigma _z^2/a^2$
obtained from the solutions of (3.6) and (3.10) in § 4.4, and quantitatively with the mean position
$\bar {z}/a_0$
obtained from these solutions in § 4.7.
Numerical results of the dimensionless average concentration distribution for a channel of different geometry (
$a/a_0$
), with
${\textit{Pe}}_a=10$
and
$\epsilon =0.001$
for various times, where the continuous lines represent the dispersion in a Newtonian fluid (
$\varepsilon {\textit{Wi}}^2=0$
) and the dotted lines represent it in a viscoelastic medium (
$\varepsilon {\textit{Wi}}^2=10$
).
Figure 9 displays the numerical results of the dimensionless average concentration distribution for a channel with different geometries, represented by distinct functions of
$a$
, with
${\textit{Pe}}_{a0}=10$
at times
$t/(a_0^2/2D)=2$
,
$20$
,
$60$
and
$100$
. The graphs in this figure include two types of lines: the continuous line indicates dispersion in a Newtonian fluid (
$\varepsilon {\textit{Wi}}^2=0$
), while the dotted line represents the behaviour in a viscoelastic fluid (
$\varepsilon {\textit{Wi}}^2=10$
). From the graphs, we can observe that the position of the mean concentration during the dispersion of a viscoelastic fluid aligns with that of a Newtonian fluid. From (3.6), viscoelasticity is not a parameter that influences the average velocity position, since the flow rate is the same, meaning the average velocity, which primarily drives the axial movement of the solute, is consistent for Newtonian and viscoelastic fluids. However, the viscoelastic effects significantly influence the variance of the solute distribution, observed graphically in figure 9, and mathematically in (3.10). In a viscoelastic fluid, the separation of the solute relative to the mean (variance) occurs more slowly compared with a Newtonian fluid. This aligns with the discussion in § 4.2, which explains how the viscoelastic effect reduces
$D_{\textit{eff}}$
, thereby impeding the ease with which the solute can disperse in the fluid.
Concerning the impact of geometry on the average concentration distribution, figure 9
$(a)$
illustrates the average concentration across a cylindrical channel with a radius of
$a/a_0=1$
. Figure 9
$(b,\!c)$
depicts the concentration distribution for linear convergent and divergent channels, respectively, while figure 9
$(d)$
represents a sinusoidal channel. For channels with a constant radius, the dimensionless position of the mean concentration is directly proportional with respect to dimensionless time. In contrast, for convergent channels, this position advances more rapidly than in constant-radius channels. This acceleration is attributable to the increase in effective diffusivity (
$D_{{\textit {eff}}}$
) in convergent channels, as discussed in § 4.2. Conversely, divergent channels decrease
$D_{{\textit {eff}}}$
, resulting in a slowdown in dispersion. For the sinusoidal channel, the mean position can either advance or delay compared with where it would be in a cylindrical channel. These results about the variance are qualitatively consistent with the analysis of variance performed in § 4.4, where it was discussed that the viscoelasticity of the fluid decreased the axial variance of the solute and that, in divergent channels, the variance evolved more slowly compared with how it does in convergent channels.
4.7. Mean concentration spatial and temporal evolution
The numerical results presented in § 4.6 illustrate how the mean concentration evolves as it moves along the channel. This movement is primarily influenced by geometric effects. Accurately determining the position of the mean concentration at any given moment is essential for various practical applications. Importantly, the average concentration distributions depicted in figure 9, for the parameters analysed, are symmetric around their mean position, resulting in approximately zero skewness. This symmetry is a consequence of the smooth and gradual changes in geometry. These findings support the model presented in (3.35), which offers a simple yet effective analytical expression for predicting the evolution of mean concentration.
Time evolution of the mean position (
$\bar {z}$
) of the solute distribution with
$\epsilon =0.001$
.
$(a)$
Numerical results of the proposed evolution of the mean concentration position given in (3.35) for different geometries.
$(b)$
Comparison between the temporal evolution of the mean position given by (3.6) represented by solid lines and that given by (3.35) represented by dotted lines, in a channel whose radius is given by
$a/a_0=1+0.2\sin {z/400a_0}$
for different Péclet numbers.
Figure 10
$(a)$
illustrates the solutions to (3.35) for various channel geometries. The plot depicts the temporal evolution of the mean concentration position, denoted as
$\bar {z}/a_0$
, where the slope of each curve corresponds to its propagation speed. This graph effectively highlights the geometric effects discussed throughout this work. The solid black line represents the mean concentration evolution in a uniform cylindrical channel, which aligns with previously reported results. The other curves demonstrate how variations in geometry impact the mean concentration position over time. Specifically, the blue dashed line represents a linearly divergent channel, while the red dash-dotted line corresponds to a linearly convergent channel. As expected, the mean concentration propagates faster in convergent geometries and slower in divergent ones. Meanwhile, the green dotted line illustrates a sinusoidal channel, where the mean concentration oscillates around the baseline of the cylindrical channel, alternating between regions with higher and lower propagation speeds. Understanding how
$\bar {z}/a_0$
depends on geometry, as well as its propagation velocity, can significantly enhance processes such as species mixing or separation in microfluidic systems.
Figure 10
$(b)$
shows a comparison between the temporal evolution of the mean position given by (3.6) represented by solid lines and that given by (3.35) represented by dotted lines, in a channel whose radius is given by
$a/a_0=1+0.2\sin {z/400a_0}$
for different Péclet numbers. The considerations made to obtain (3.35) were very low skewness values and
${\textit{Pe}}_{a0}/\epsilon \gg 1$
. The limiting case to analyse the validity of (3.35) occurs for low values of
${\textit{Pe}}_{a0}$
that may put at risk that
${\textit{Pe}}_{a0}/\epsilon \gg 1$
is fulfilled; however, the validity of Taylor’s analysis occurs for values of
${\textit{Pe}}_{a_0}\gt 6.9$
according with experimental data (Taylor Reference Taylor1953), so the lowest case analysed is
${\textit{Pe}}_{a0}=10$
. We can observe a good approximation between both solutions, even for values of
${\textit{Pe}}_{a0}=10$
, allowing us to validate our proposal given in (3.35), which allows us to easily obtain the temporal position of the average concentration along a channel of variable geometry.
5. Conclusions
The combined effects of viscoelasticity and channel geometry significantly alter the velocity distribution of the fluid flow. Viscoelastic effects lead to a flatter velocity profile and shear thinning, which reduces the velocity magnitude in the centre of the channel while maintaining a constant flow rate. In contrast, the geometry of the channel can either decrease the velocity if the channel expands or increase it if the channel narrows. These observations are consistent with theoretical expectations and offer valuable insight into the complex relationship between fluid rheology and conduit shape in non-Newtonian flows.
The analysis demonstrates that the effective dispersion coefficient in viscoelastic flows is primarily influenced by two competing factors: fluid viscoelasticity and channel geometry. As the viscoelastic parameter
$\varepsilon {\textit{Wi}}^2$
increases, there is a significant reduction in the effective dispersion coefficient (
$D_{{\textit {eff}}}$
) throughout the channel. This reduction is primarily attributed to the flattening of the velocity profile, which approaches a plug-like flow. This change reduces both the velocity gradients and the effects of radial diffusion. The geometric characteristics of the channel further modify this behaviour. In diverging channels, the effective dispersion coefficient decreases as the cross-section widens, primarily due to slower fluid velocities. Conversely, converging channels exhibit the opposite trend, where the effective dispersion coefficient increases. This pattern can also be applied to channels that vary in size, such as sinusoidal channels.
We obtained a system of ordinary differential equations that allow us to find the average position and axial variance of a solute in a viscoelastic fluid for a given geometry that changes smoothly axially. To have a zero variance transient growth, for the case of viscoelastic fluid, smaller axial variations of the channel radius are required (lower values of
$\beta$
) compared with a Newtonian fluid, for a given Péclet number, and the value of
${\textit{Pe}}_{a}$
at which the smallest variance occurs is different and greater than
$\sqrt {48}$
which occurs in a Newtonian fluid. The variance was analysed for different geometries, finding that it decreases due to the viscoelastic effect and that it has a behaviour that differs from the direct behaviour of
$D_{{\textit {eff}}}/D$
, and that the solution matches the Brownian simulations successfully, except for values of
${\textit{Pe}}_{a0}=1000$
where our model differs from the Brownian simulations. The analytical approach to design an engineering channel for a target variance has been shown to have good results compared with Brownian simulations, being a powerful tool capable of designing channel geometry from a desired axial variance distribution of the solute.
Solutions for the distribution of average concentration have been obtained based on a proposal based on
${\textit{Pe}}_{a0}/\epsilon \gg 1$
values, which turned out to be successful in predicting the spatial and temporal distribution of average concentration. This solution demonstrates that while viscoelasticity reduces solute dispersion by decreasing the effective diffusion coefficient (
$D_{{\textit {eff}}}$
) and slowing the growth of variance, it maintains the transport velocity of the mean concentration, as flow rates are consistent between the Newtonian and viscoelastic cases, coinciding with the analysis of variance and position carried out in the work. Channel geometry is a crucial parameter for controlling dispersion; converging channels enhance spreading by increasing
$D_{{\textit {eff}}}$
whereas diverging geometries suppress it. Sinusoidal variations result in sinusoidal dispersion patterns. The combined use of transformed coordinates (
$\zeta ,\tau$
) and numerical methods validates a practical framework for analysing dispersion in axially varying geometries.
Finally, a proposal was made to determine the position of the mean concentration at any given moment, based on the consideration of skewness
$\ll 1$
. This study highlights the significant role of channel geometry in determining the evolution of mean concentration, while viscoelastic effects are found to be negligible. The model proposed in this work provides a simple yet effective tool for predicting concentration dynamics across various geometries, including uniform, convergent, divergent and sinusoidal channels, as well as any other geometries with slowly varying changes, which was successfully validated with the solution of a system of ODEs for the mean position and axial variance.
In summary, these findings offer essential insights for microfluidic design, suggesting that a combination of geometric and rheological adjustments can create tailored dispersion profiles. This can be beneficial for applications ranging from chemical processing, such as industrial polymer processing, to biomedical devices where controlled mixing and transport are critical. Future research could build on this knowledge by exploring more complex rheological models and examining transient, oscillatory or pulsating flow conditions to establish broader applicability for these dispersion mechanisms.
Acknowledgements
C.T. acknowledges the support of the Secretaria de Ciencias Humanidades Tecnología e Innovación program for Postdoctoral fellowship at the Escuela Superior de Ingeniería Mecánica y Eléctrica Azcapotzalco from Instituto Politécnico Nacional de México. O.B. expresses their gratitude to the Secretaría de Investigación y Posgrado at the Instituto Politécnico Nacional in México for Research Grants 20250739, as well as to the Red de Energía and the Red de Nanociencia y Micro-nanotecnología.
Funding
This work was supported by the Estancias Posdoctorales por México Consejo Nacional de Humanidades, Ciencia y Tecnología through research grants with application numbers 4923376 and CVU 895339.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the hydrodynamic field
In its scalar form, for the axisymmetric case, (2.1)–(2.3) can be written as
\begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \frac {1}{r}\frac {\partial }{\partial r}(ru_r)+\frac {\partial u_z}{\partial z}=0, \end{align}
\begin{align}&\qquad\qquad \rho \left ( \frac {\partial u_r}{\partial t}+u_r\frac {\partial u_r}{\partial r}+u_z\frac {\partial u_r}{\partial z} \right )=-\frac {\partial p}{\partial r}+\frac {1}{r}\frac {\partial }{\partial r}(r\tau _{rr})+\frac {\partial \tau _{zr}}{\partial z}, \end{align}
\begin{align}&\qquad\qquad \rho \left ( \frac {\partial u_z}{\partial t}+u_r\frac {\partial u_z}{\partial r}+u_z\frac {\partial u_z}{\partial z} \right )=-\frac {\partial p}{\partial z}+\frac {1}{r}\frac {\partial }{\partial r}(r\tau _{rz})+\frac {\partial \tau _{zz}}{\partial z}, \end{align}
\begin{align}& f(\tau _{\textit{kk}})\tau _{rr}+\lambda \left \{ \frac {\partial \tau _{rr}}{\partial t}+ u_r \frac {\partial \tau _{rr}}{\partial r}+u_z \frac {\partial \tau _{rr}}{\partial z} -2\left ( \tau _{rr}\frac {\partial u_r}{\partial r}+\tau _{rz}\frac {\partial u_r}{\partial z} \right )\right \}=2\eta \frac {\partial u_r}{\partial r}, \end{align}
\begin{align}& f(\tau _{\textit{kk}})\tau _{zz}+\lambda \left \{ \frac {\partial \tau _{zz}}{\partial t}+ u_r \frac {\partial \tau _{zz}}{\partial r}+u_z \frac {\partial \tau _{zz}}{\partial z} -2\left ( \tau _{zz}\frac {\partial u_z}{\partial z}+\tau _{rz}\frac {\partial u_z}{\partial r} \right )\right \}=2\eta \frac {\partial u_z}{\partial z} \\[2pt] \nonumber \end{align}
and
\begin{align} &f(\tau _{\textit{kk}})\tau _{zr}+\lambda \left \{ \frac {\partial \tau _{zr}}{\partial t}+ u_r \frac {\partial \tau _{rz}}{\partial r}+u_z \frac {\partial \tau _{rz}}{\partial z}-\tau _{rr}\frac {\partial u_z}{\partial r}-\tau _{zz}\frac {\partial u_r}{\partial z} -\tau _{rz}\left ( \frac {\partial u_r}{\partial r}+\frac {\partial u_z}{\partial z} \right )\right \}\nonumber\\&\quad =\eta \left (\frac {\partial u_z}{\partial r}+\frac {\partial u_r}{\partial z}\right )\!, \end{align}
with
$f(\tau _{\textit{kk}})=1+\varepsilon \lambda (\tau _{rr}+\tau _{zz})/\eta$
and are subject to
\begin{equation} u_z(r=0,z) \text{, } u_r(r=0) \text{ and } \tau _{rz}(r=0,z) \text{ are finite,}\hspace {1cm} u_z(r=a,z)=u_r(r=a,z)=0. \end{equation}
Using scale analysis with the characteristic scales outlined in § 4.2, we neglect the terms of order
$\epsilon$
,
$\epsilon ^2$
,
$\epsilon Re$
,
$\epsilon ^2 Re$
,
$\epsilon Wi$
and
$\epsilon ^2 Wi$
. The equations that follow are presented in dimensionless form, but for simplicity, we omit the superscript
$^*$
. As a result, (A1)–(A6) can be written as
\begin{align}&\quad \frac {1}{r}\frac {\partial }{\partial r}(ru_r)+\frac {\partial u_z}{\partial z}=0, \end{align}
\begin{align}&\qquad\qquad \frac {\partial p}{\partial r}=0, \end{align}
\begin{align}&\qquad\, \frac {\partial p}{\partial z}=\frac {1}{r}\frac {\partial }{\partial r}(r\tau _{rz}), \end{align}
\begin{align}&\qquad\quad f(\tau _{\textit{kk}})\tau _{rr}=0, \end{align}
\begin{align}& f(\tau _{\textit{kk}})\tau _{zz}-2Wi \tau _{rz}\frac {\partial u_z}{\partial r}=0 \end{align}
and
\begin{equation} f(\tau _{\textit{kk}})\tau _{zr}-Wi \tau _{rr}\frac {\partial u_z}{\partial r}= \frac {\partial u_z}{\partial r}. \end{equation}
Integrating (A10) and using the boundary condition given in (A7), we obtain
\begin{equation} \tau _{rz}=\frac {r}{2}\frac {{\rm d}p}{{\rm d}z}, \end{equation}
combining (A12) and (A13), we obtain
\begin{equation} \tau _{zz}=2Wi\tau _{rz}^{2}, \end{equation}
and finally replacing (A14) and (A15) in (A13), we have
\begin{equation} \left (1+2\varepsilon {\textit{Wi}}^2 \left (\frac {r}{2}\frac {{\rm d}p}{{\rm d}z} \right )^2\right )\left (\frac {r}{2}\frac {{\rm d}p}{{\rm d}z}\right )=\frac {\partial u_z}{\partial r}. \end{equation}
Solving (A16) and applying the no-slip condition given in (A7), we obtain
\begin{equation} u_z=\left (-\frac {a^{2}}{4}\frac {{\rm d}p}{{\rm d}z}\right )\left ( 1-\frac {r^{2}}{a^{2}}\right )+\frac {\varepsilon {\textit{Wi}}^2 a^{4}}{16}\left (-\frac {{\rm d}p}{{\rm d}z} \right )^3\left (1-\frac {r^{4}}{a^{4}} \right )\!. \end{equation}
If we calculate the average velocity from the volumetric flow rate of the fluid, we can find that
\begin{equation} U(z)=\frac {1}{\pi a^2 }\int _0^{a}u_z2\pi r \,{\rm d}r=\left (-\frac {a^{2}}{8}\frac {{\rm d}p}{{\rm d}z} \right )+ \frac {\varepsilon {\textit{Wi}}^2 a^{4}}{24}\left (-\frac {{\rm d}p}{{\rm d}z} \right )^3=U_n(z)+U_{ve}(z). \end{equation}
Writing (A17) in terms of these two averages, we have
\begin{equation} u_z=2U_n\left ( 1-\frac {r^{2}}{a^{2}}\right )+\frac {3}{2}U_{ve}\left (1-\frac {r^{4}}{a^{4}} \right )\!. \end{equation}
However, it is necessary to know
$u_r$
induced by changes in the microchannel geometry. To do this, it is necessary to resort to the conservation of mass given in (A8) and the conditions given in (A7). If we differentiate (A17) with respect to
$z$
, we get
\begin{equation} \frac {\partial u_z}{\partial z}=2\frac {{\rm d} U_n}{{\rm d} z}\left (1-\frac {r^2}{a^2} \right )+\frac {4r^2\beta U_n}{a^3}+\frac {3}{2}\frac {{\rm d} U_{ve}}{{\rm d} z}\left (1-\frac {r^4}{a^4} \right )+\frac {6r^4\beta U_{ve}}{a^5}. \end{equation}
Here,
${\rm d}U_n/{\rm d}z$
and
${\rm d}U_{ve}/{\rm d}z$
can be obtained by separating the velocity into its Newtonian part and its viscoelastic part
$u_z=u_{zn}+u_{zve}$
and
$u_r=u_{rn}+u_{rve}$
as follows:
\begin{align} \frac {{\rm d} U_n}{{\rm d} z}=\frac {d}{{\rm d} z}\left (\frac {2}{a^2}\int _0^a u_{zn}r\,{\rm d}r \right )=\frac {-4\beta }{a^3}\int _0^a u_{zn}r\,{\rm d}r+\frac {2}{a^2}\frac {d }{{\rm d} z}\int _0^a u_{zn}r\,{\rm d}r, \end{align}
of which the last derivative can be obtained using Leibniz’s rule and the no-slip condition given in (A7):
\begin{align} \frac {d }{{\rm d} z}\int _0^a u_{zn}r\,{\rm d}r=\int _0^a \frac {\partial (u_{zn}r) }{\partial z}\,{\rm d}r+\underbrace {(u_{zn}r)|_a \frac {{\rm d}a}{{\rm d}z}}_{=0}- \underbrace {(u_{zn}r)|_0 (0)}_{=0}. \end{align}
By continuity,
$\int _0^a \partial (u_{zn}r) / \partial z\, {\rm d}r=-(u_{rn}r)|_0^a=0$
, so all terms of the integral are zero, so
\begin{equation} \frac {{\rm d} U_n}{{\rm d} z}=\frac {-4\beta }{a^3}\int _0^a u_{zn}r\,{\rm d}r=\frac {-2\beta }{a}\left (\frac {2}{a^2}\int _0^au_{zn}r\,{\rm d}r \right )=\frac {-2\beta }{a}U_n. \end{equation}
Similarly,
\begin{equation} \frac {{\rm d} \bar {U}_{ve}}{{\rm d} z}=\frac {-2\beta }{a}U_{ve}. \end{equation}
By continuity and using (A20), (A23) and (A24), we have
\begin{equation} u_r=-\frac {1}{r}\int \frac {\partial u_z}{\partial z} r \,{\rm d}r =-\frac {1}{r}\int \left ( \frac {4\beta \bar {U}_n}{a}\left (\frac {2r^2}{a^2}-1 \right )+\frac {3\beta \bar {U}_{ve}}{a}\left (\frac {3r^4}{a^4}-1 \right )\right ) r \,{\rm d}r. \end{equation}
Finally, by solving (A25) and applying the condition at the centre of the channel given in (A7), we obtain
\begin{equation} u_r=2\beta U_n\left (\frac {r}{a}-\frac {r^3}{a^3} \right )+\frac {3}{2}\beta U_{ve}\left (\frac {r}{a}-\frac {r^5}{a^5} \right )\!. \end{equation}
Appendix B. Simplification of terms in the average of the convective-diffusion equation
Using Leibniz’s rule and the chain rule, some terms of (2.13) are developed as follows:
-
(i)
$\langle {\partial c'}/{\partial r} \rangle$
(B1)separating the integrals and evaluating the first one results in
\begin{align} \frac {2}{a^2}\int _0^{a} \frac {\partial c'}{\partial r}r\,{\rm d}r=\frac {2}{a^2}\int _0^{a} \left (\frac {\partial }{\partial r}(rc')-c' \right )\,{\rm d}r, \end{align}
(B2)
\begin{align} \frac {2}{a^2} \left \lbrace \int _0^{a} \frac {\partial (rc')}{\partial r}\,{\rm d}r- \int _0^{a}c'\, {\rm d}r \right \rbrace =\frac {2}{a^2} \left \lbrace ac'(a)- \int _0^{a}c' \,{\rm d}r \right \rbrace ; \end{align}
-
(ii)
$\langle {\partial c'}/{\partial r} \rangle$
(B3)by Leibniz’s rule,
\begin{align} \frac {2}{a^2}\int _0^{a} \frac {\partial c'}{\partial z}r\,{\rm d}r=\frac {2}{a^2}\int _0^{a} \frac {\partial (rc')}{\partial z} \,{\rm d}r, \end{align}
(B4)therefore,
\begin{align} \frac {d}{{\rm d}z}\int _0^{a}rc' \,{\rm d}r=\frac {2}{a^2}\int _0^{a} \frac {\partial (rc')}{\partial z} \,{\rm d}r+(rc')|_{r=a}\beta , \end{align}
(B5)
\begin{align} \frac {2}{a^2}\int _0^{a} \frac {\partial (rc')}{\partial z} \,{\rm d}r=\frac {2}{a^2} \left [\frac {d}{{\rm d}z}\int _0^{a}rc'\, {\rm d}r- ac'(a)\beta \right ] =\frac {\rm d}{{\rm d}z}\langle c' \rangle -\frac {2c'(a)\beta }{a}=-\frac {2c'(a)\beta }{a}; \end{align}
-
(iii)
$\langle ({1}/{r})({\partial }/{\partial r}) ( r( {\partial c'}/{\partial r} )) \rangle$
(B6)
\begin{align} \frac {2}{a^2}\int _0^{a} \frac {1}{r}\frac {\partial }{\partial r} \left ( r\frac {\partial c'}{\partial r} \right )r\,{\rm d}r=\frac {2}{a^2}\int _0^{a}\frac {\partial }{\partial r} \left ( r\frac {\partial c'}{\partial r} \right )\,{\rm d}r =\frac {2}{a^2}\left [a\frac {\partial c'}{\partial r}\bigg |_{r=a}-0 \right ]=\frac {2}{a}\frac {\partial c'}{\partial r}\bigg |_{r=a}; \end{align}
-
(iv)
$\langle {\partial ^2 c'}/{\partial z^2} \rangle$
(B7)according to Leibniz’s rule, we have that
\begin{align} \frac {2}{a^2}\int _0^{a}\frac {\partial ^2 c'}{\partial z^2}r\,{\rm d}r=\frac {2}{a^2}\int _0^{a}\frac {\partial ^2 (rc')}{\partial z^2}\,{\rm d}r, \end{align}
(B8)Using Leibniz’s rule again, it turns out
\begin{align} \frac {2}{a^2}\int _0^{a}\frac {\partial ^2 (rc')}{\partial z^2}\,{\rm d}r=\frac {2}{a^2}\left [ \frac {\rm d}{{\rm d}z}\int _0^{a}\frac {\partial (rc')}{\partial z}\,{\rm d}r-\beta \frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right ]\!. \end{align}
(B9)the indicated term is equal to 0, assuming that
\begin{align} \frac {2}{a^2}\left [ \frac {\rm d}{{\rm d}z}\left [\frac {\rm d}{{\rm d}z}\underbrace {\int _0^{a}c'r\,{\rm d}r}_{=0}-(rc')|_{r=a}\beta \right ]-\beta \frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right ]\!, \end{align}
$\langle c'\rangle =0$
. Continuing with the simplification, we have(B10)considering that
\begin{align} &-\frac {2}{a^2}\left [ \frac {\rm d}{{\rm d}z}\left [(rc')|_{r_=a}\beta \right ]+\beta \frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right ] \nonumber \\&\quad =-\frac {2}{a^2}\left [\frac {\rm d}{{\rm d}z}\left [(rc')|_{r_=a}\right ]\beta +(rc')|_{r_=a}\gamma +\beta \frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right ]\!, \end{align}
${\rm d}/{\rm d}z [(rc')|_{r_=a} ]=\partial (rc')/\partial z|_{r=a}$
, finally, it can be written as(B11)
\begin{align} -\frac {2}{a^2}\left [(rc')|_{r_=a}\gamma +2\beta \frac {\partial (rc')}{\partial z}\bigg |_{r=a} \right ] =-\frac {2\gamma c'(a)}{a}-\frac {4\beta }{a^2}\frac {\partial (rc')}{\partial z}\bigg |_{r=a}. \end{align}
