2.1. Governing equations

Consider a vertical wavy channel of average gap size $h_o$ filled with a Newtonian fluid of density $\rho$ and kinematic viscosity $\nu$ (for CSF, $\rho \simeq 10^{3}$ kg m$^{-3}$ and $\nu \simeq 0.7 \times 10^{-6}$ m$^{2}$ s$^{-1}$). The channel, open at both ends, is bounded by a flat surface and a wavy wall of wave length $\lambda \gg h_o$, so that the resulting channel width is $h=h_o[1+\beta \cos (2 {\rm \pi}x^{*}/\lambda )]$, where $x^{*}$ is the longitudinal distance measured from the upper end and $\beta < 1$ is the relative amplitude of the wall undulation. The total channel length is $n \lambda$, with $n$ representing a general integer number, so that the channel comprises $n$ identical cells. The flow is described using cartesian coordinates $(x^{*},y^{*})$, with $y^{*}$ measured from the flat surface, and corresponding velocity components $\boldsymbol{v}^{*}=(u^{*},v^{*})$. The Navier–Stokes equations describing the planar unsteady flow are written in the Boussinesq approximation

(2.1)\begin{gather} \boldsymbol{\nabla}^{*} \boldsymbol{\cdot} \boldsymbol{v}^{*} =0, \end{gather}

(2.2)\begin{gather}\frac{\partial \boldsymbol{v}^{*}}{\partial t^{*}}+\boldsymbol{v}^{*} \boldsymbol{\cdot} \boldsymbol{\nabla}^{*} \boldsymbol{v}^{*} ={-}\frac{1}{\rho} \boldsymbol{\nabla}^{*} p^{*}+ \nu \boldsymbol{\nabla}^{*2} \boldsymbol{v}^{*} - \frac{\rho-\rho_s}{\rho} g c \boldsymbol{e}_x, \end{gather}

(2.3)\begin{gather}\frac{\partial c}{\partial t^{*}}+\boldsymbol{v}^{*} \boldsymbol{\cdot} \boldsymbol{\nabla}^{*} c=\kappa \boldsymbol{\nabla}^{*2} c, \end{gather}

where $p^{*}$ is the sum of the pressure difference from the upper end and the constant-density hydrostatic component $-\rho g x^{*}$, $c$ is the solute volume concentration, $\kappa$ is the solute diffusivity, $\boldsymbol {\nabla }^{*}=(\partial /\partial x^{*},\partial /\partial y^{*})$ and $\boldsymbol{e}_x$ is the unit vector aligned with the gravitational acceleration.

A pressure difference $n \Delta p \cos (\omega t^{*})$ oscillating harmonically in time is prescribed between the upper and lower ends of the canal, driving a periodic fluid motion with angular frequency $\omega$. The resulting slender flow is characterised by longitudinal velocities of order $u_c=\Delta p/(\rho \omega \lambda )$, as follows from a balance between the local acceleration $\partial u^{*}/\partial t^{*} \sim u_c \omega$ and the pressure gradient $\rho ^{-1} \partial p^{*}/\partial x^{*} \sim \Delta p/(\lambda \rho )$, and much smaller transverse velocities of order $v_c=(h_o/\lambda )u_c \ll u_c$, as follows from the continuity balance $\partial u^{*}/\partial x^{*} \sim \partial v^{*}/\partial y^{*}$.

2.2. Controlling parameters

The analysis assumes that the viscous time across the channel $h_o^{2}/\nu$ is comparable to the characteristic oscillation time $\omega ^{-1}$, resulting in Womersley numbers

(2.4)\begin{equation} \alpha=\left(\frac{\omega h_o^{2}}{\nu}\right)^{1/2}, \end{equation}

of order unity. The limit $\alpha \sim 1$ is instrumental in analysing cardiac-driven CSF flow ($\omega =2{\rm \pi}$ s$^{-1}$) in the spinal canal, for which typical values of $\alpha$ are in the range $3 \lesssim \alpha \lesssim 6$, as can be seen by evaluating the above expression with $h_o \simeq 1\unicode{x2013}2$ mm and $\nu =0.7 \times 10^{-6}$ m$^{2}$ s$^{-1}$. In the lumbar region, the typical drug-delivery site in ITDD procedures, the average CSF speeds are of order $u_c \sim 1$ cm s$^{-1}$, so that the associated stroke lengths $u_c/\omega$ are much smaller than the characteristic longitudinal distance $\lambda \simeq 2\unicode{x2013}4$ cm. Their ratio

(2.5)\begin{equation} \varepsilon = \frac{u_c/\omega}{\lambda}, \end{equation}

of order $\varepsilon \simeq 0.05$ for spinal CSF flow, defines the small parameter employed in the following asymptotic description. As shown earlier (Hall Reference Hall1974; Grotberg Reference Grotberg1984), the solution at leading order is determined by a balance between the pressure gradient, the local acceleration and the viscous forces, with the convective acceleration introducing small corrections of relative magnitude $\varepsilon$. Although the leading-order motion is harmonic, the velocity corrections include a non-zero steady-streaming component.

The familiar periodic channel flow described previously is altered by gravitational forces when a solute of density $\rho _s \ne \rho$ is introduced in the channel. The extent of the resulting buoyancy-induced motion can be measured by the associated Richardson number

(2.6)\begin{equation} {{\textit{Ri}}}=\frac{\rho-\rho_s}{\rho} \frac{g \lambda}{u_c^{2}}, \end{equation}

which compares the order of magnitude of the effective gravitational acceleration $g (\rho -\rho _s)/\rho$ with that of the convective acceleration $\boldsymbol {v}^{*} \boldsymbol {\cdot} \boldsymbol {\nabla }^{*} \boldsymbol {v}^{*} \sim u_c^{2}/\lambda$. Our analysis addresses the limit ${{\textit {Ri}}} \sim 1$, which is relevant for drug dispersion in ITDD procedures, as can be seen by evaluating (2.6) with $\lambda \simeq 2$ cm and $u_c \sim 1$ cm s$^{-1}$ for density differences in the range $10^{-3} \lesssim |\rho -\rho _s|/\rho \lesssim 10^{-2}$.

Also motivated by ITDD applications, we consider solutes with diffusivities $\kappa$ much smaller than the kinematic viscosity, that always being the case of diffusion in liquid phase. As $\nu /\kappa \sim 1000$ and $\varepsilon \simeq 0.05$ in ITDD applications, the following analysis of solute dispersion will specifically address the distinguished limit $\kappa /\nu \sim \varepsilon ^{2}$, with solute diffusion correspondingly characterised by the reduced Schmidt number

(2.7)\begin{equation} \sigma=\varepsilon^{2} \frac{\nu}{\kappa}, \end{equation}

assumed to be of order unity.

2.3. Non-dimensional formulation

We address the motion that follows from the deposition of the solute inside an intermediate cell along the channel. The problem is non-dimensionalised using the scales identified previously to give the dimensionless variables

(2.8a–f)\begin{equation} t=\omega t^{*}, \quad x=\frac{x^{*}}{\lambda}, \quad y=\frac{y^{*}}{h_o}, \quad u=\frac{u^{*}}{u_c}, \quad v=\frac{v^{*}}{v_c}, \quad p=\frac{p^{*}}{\Delta p} \end{equation}

and associated conservation equations

(2.9)\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, \end{gather}

(2.10)\begin{gather}\frac{\partial u}{\partial t}+\varepsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} u={-}\frac{\partial p}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} u}{\partial y^{2}}- \varepsilon \,{{\textit{Ri}}} \,c, \end{gather}

(2.11)\begin{gather}\frac{\partial p}{\partial y}=0, \end{gather}

(2.12)\begin{gather}\frac{\partial c}{\partial t}+\varepsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} c = \frac{\varepsilon^{2}}{\alpha^{2} \sigma} \frac{\partial^{2} c}{\partial y^{2}}, \end{gather}

where $\boldsymbol {v}=(u,v)$ and $\boldsymbol {\nabla }=(\partial /\partial x,\partial /\partial y)$. In writing (2.9)–(2.12) from (2.1)–(2.3) we have used the slender-flow approximation resulting from the limit $h_o \ll \lambda$. Thus, the terms representing longitudinal diffusion of momentum and mass have been neglected in (2.10) and (2.12), because they are a factor $(h_o/\lambda )^{2}$ smaller than those associated with transverse diffusion. At the same level of approximation, the transverse component of the momentum equation takes the reduced form (2.11). The velocity and concentration must satisfy the boundary conditions

(2.13)\begin{equation} u=v=\frac{\partial c}{\partial y}=0 \quad {\rm at} \ \left\{ \begin{array}{@{}l} y=0 \\ y=H=1+\beta \cos(2 {\rm \pi}x) \end{array} \right., \end{equation}

corresponding to non-permeable no-slip surfaces, whereas the reduced pressure $p(x,t)$, independent of $y$, is identically zero at $x=0$ and takes the value $p=n \cos t$ at the lower end $x=n$.

In the absence of buoyancy (i.e. for ${{\textit {Ri}}}=0$), the solution for the velocity is periodic in time, including a steady component of order $\varepsilon$, and also periodic in space, so that the velocity distribution found in each cell is identical. On the other hand, for ${{\textit {Ri}}} \ne 0$ the motion is coupled to the solute transport, albeit weakly, with the result that the velocity necessarily evolves in time following the dispersion of the solute, which is driven partly by the steady streaming motion, with characteristic velocities $\varepsilon u_c$. It can be anticipated that the characteristic time for the slow evolution is that associated with the dispersion of the solute inside the deposition cell $\lambda /(\varepsilon u_c)=\varepsilon ^{-2} \omega ^{-1}$, much larger than the characteristic oscillation time $\omega ^{-1}$. These considerations suggest the introduction of a second time variable $\tau = \varepsilon ^{2} t$ for describing the slow evolution, additional to the fast time-scale variable $t$ describing the oscillatory motion. In the two-time-scale formalism, the time derivatives in (2.10) and (2.12) are replaced with $\partial /\partial t+\varepsilon ^{2} \partial /\partial \tau$ and the different variables are expressed in terms of the power expansions

(2.14)\begin{gather} u=u_0(x,y,t,\tau)+\varepsilon u_1(x,y,t,\tau)+\cdots, \end{gather}

(2.15)\begin{gather}v=v_0(x,y,t,\tau)+\varepsilon v_1(x,y,t,\tau)+\cdots, \end{gather}

(2.16)\begin{gather}p=p_0(x,t,\tau)+\varepsilon p_1(x,t,\tau)+\cdots, \end{gather}

(2.17)\begin{gather}c=c_0(x,y,t,\tau)+\varepsilon c_1(x,y,t,\tau)+\cdots, \end{gather}

with all functions assumed to be $2{\rm \pi}$ periodic in the fast time scale $t$. The asymptotic procedure leads to a hierarchy of problems that can be solved sequentially, as shown in the following.

3.1. Leading-order oscillatory flow

Convective acceleration and buoyancy are negligible at leading order, so that the velocity $\boldsymbol{v}_0=(u_0,v_0)$ satisfies a linear problem that can be solved in terms of the reduced variables

(3.1a–c)\begin{equation} u_0={\rm Re}(\mathrm{i} \mathrm{e}^{\mathrm{i} t} U),\quad v_0={\rm Re} (\mathrm{i} \mathrm{e}^{\mathrm{i} t} V),\quad p_0={\rm Re} (\mathrm{e}^{\mathrm{i} t} P), \end{equation}

where the complex functions $U(x,y)$, $V(x,y)$, and $P(x)$ satisfy

(3.2a,b)\begin{equation} \frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0 \quad {\rm and} \quad -U ={-}\frac{{\rm d} P}{{\rm d} x}+\frac{\mathrm{i}}{\alpha^{2}} \frac{\partial^{2} U}{\partial y^{2}}. \end{equation}

The second equation above can be integrated with boundary conditions $U=0$ at $y=0,H$ to give

(3.3)\begin{equation} U=\frac{{\rm d} P}{{\rm d} x} G(x,y), \end{equation}

where

(3.4a,b)\begin{equation} G=1-\frac{\cosh[\varLambda(2y/H-1)]}{\cosh \varLambda} \quad {\rm and} \quad \varLambda=\frac{\alpha}{2} \frac{1+\mathrm{i}}{\sqrt{2}} H(x). \end{equation}

The result can be used to integrate the first equation in (3.2a,b) subject to $V=0$ at $y=0$, yielding

(3.5)\begin{equation} V={-}\frac{\partial}{\partial x} \left(\int_0^{y} U \,{\rm d} \hat{y} \right)={-}\frac{\partial}{\partial x} \left(\frac{{\rm d} P}{{\rm d} x} \int_0^{y} G\, {\rm d}\hat{y} \right), \end{equation}

where

(3.6)\begin{equation} \int_0^{y} G \,{\rm d}\hat{y}=y-\frac{H}{2 \varLambda} \left\{\frac{\sinh[\varLambda(2y/H-1)]+\sinh \varLambda}{\cosh \varLambda} \right\}, \end{equation}

with $\hat {y}$ representing a dummy integration variable. Note that both velocity components $U$ and $V$ are spatially periodic in $x$ through the function $H=1+\beta \cos (2 {\rm \pi}x)$. The determination of the longitudinal pressure gradient ${\rm d}P/{\rm d} x$ that completes the solution begins by using the condition $V=0$ at $y=H$ in the first equation of (3.5) to give

(3.7)\begin{equation} \frac{{\rm d}}{{\rm d} x} \left(\int_0^{H} U \,{\rm d} y \right)=0, \end{equation}

indicating that the reduced flow rate

(3.8)\begin{equation} Q=\int_0^{H} U \,{\rm d} y=\frac{{\rm d} P}{{\rm d} x} \int_0^{H} G\, {\rm d} y \end{equation}

is constant. Further progress requires use of the conditions $P(0)=P(n)-n=0$, consistent with the boundary values $p(0,t)=0$ and $p(n,t)=n \cos t$ stated below (2.13). Using (3.6) to evaluate the integral $\int _0^{H} G \,{\rm d} y$ leads to the equation

(3.9)\begin{equation} Q=\frac{{\rm d} P}{{\rm d} x} H (1-\varLambda^{{-}1} \tanh \varLambda), \end{equation}

which can be integrated subject to $P(0)=0$ to yield the pressure distribution

(3.10)\begin{equation} P=Q \int_0^{x} \frac{{\rm d} \hat{x}}{H (1-\varLambda^{{-}1} \tanh \varLambda)}. \end{equation}

Using now the condition $P(n)=n$ provides

(3.11)\begin{equation} Q=n \left[\int_0^{n} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)} \right]^{{-}1}. \end{equation}

Owing to the spatial periodicity of $H$ and $\varLambda$, it follows that

(3.12)\begin{equation} \int_0^{n} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}=n \int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}, \end{equation}

thereby finally yielding

(3.13)\begin{equation} Q=\left[\int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}\right]^{{-}1} \end{equation}

and, from (3.9),

(3.14)\begin{equation} \frac{{\rm d} P}{{\rm d} x}=\left[H (1-\varLambda^{{-}1} \tanh \varLambda) \int_0^{1} \frac{{\rm d} x}{H (1-\varLambda^{{-}1} \tanh \varLambda)}\right]^{{-}1}, \end{equation}

independent of $n$. It is worth pointing out that, because at this order the velocity is harmonic, the associated time-averaged values $\langle u_0 \rangle$ and $\langle v_0 \rangle$ with $\langle {\cdot } \rangle = ({1}/{2 {\rm \pi}}) \int _t^{t+2 {\rm \pi}}\,{\rm d}t$ are identically zero, so that the steady bulk motion of the fluid occurs through the velocity corrections at the following order.

3.2. First-order corrections

Collecting terms of order $\varepsilon$ in (2.9) and (2.10) yields

(3.15)\begin{gather} \frac{\partial u_1}{\partial x}+\frac{\partial v_1}{\partial y}=0, \end{gather}

(3.16)\begin{gather}\frac{\partial u_1}{\partial t}+ \frac{\partial}{\partial x}(u_0^{2})+ \frac{\partial}{\partial y}(u_0 v_0)={-}\frac{\partial p_1}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} u_1}{\partial y^{2}}- {{\textit{Ri}}} \, c_0, \end{gather}

to be integrated with boundary conditions $u_1=v_1=0$ at $y=0,H$ and $p_1=0$ at $x=0,n$. There is interest in computing the corresponding time-averaged velocity correction $\langle \boldsymbol{v}_1 \rangle =(\langle u_1 \rangle,\langle v_1 \rangle )$. Taking the time average of (3.15) and (3.16) provides

(3.17a,b)\begin{equation} \frac{\partial \langle u_1 \rangle}{\partial x}+\frac{\partial \langle v_1 \rangle}{\partial y}=0 \quad {\rm and} \quad F(x,y)={-}\frac{\partial \langle p_1 \rangle}{\partial x}+\frac{1}{\alpha^{2}} \frac{\partial^{2} \langle u_1 \rangle}{\partial y^{2}}- {{\textit{Ri}}} \, c_0, \end{equation}

where the known function $F=\partial \langle u_0^{2} \rangle /\partial x+\partial \langle u_0 v_0 \rangle /\partial y$ can be expressed in terms of the complex velocities $U$ and $V$ defined above in the form

(3.18)\begin{equation} F=\frac{1}{2} {\rm Re}\left[\frac{\partial}{\partial x}(U \bar{U})+\frac{\partial}{\partial y}(V \bar{U}) \right], \end{equation}

a result following from the identity $\langle {\rm Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} f_1) \, {\rm Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} f_2) \rangle = {\rm Re}(f_1 \bar {f}_2)/2$, which applies to any generic time-independent complex functions $f_1$ and $f_2$, with the bar denoting complex conjugates. In writing (3.17a,b), we have anticipated that, at leading order, the solute concentration is independent of the fast time scale $t$, as follows from (2.12) when $\varepsilon \ll 1$, so that its time-averaged value $\langle c_0 \rangle$ reduces simply to $\langle c_0 \rangle =c_0$. Also of interest is that, because of the symmetry of $H(x)$, the periodic function $F$ defined in (3.18) is antisymmetric with respect to $x=1/2$, so that $F(x,y)=-F(1-x,y)$.

As can be concluded from (3.16), the velocity corrections arise partly owing to the convective acceleration and partly owing to the buoyancy force. In computing the corresponding time average, it is convenient to compute both contributions separately by introducing $\langle \boldsymbol{v}_1 \rangle =\boldsymbol{v}_{SS}+\boldsymbol{v}_{B}$ and $\langle p_1 \rangle =p_{SS}+p_{B}$, where $\boldsymbol{v}_{SS}(x,y)=(u_{SS},v_{SS})$ and $p_{SS}(x)$ describe the familiar steady-streaming associated with the nonlinear convective terms, which is independent of time and periodic in $x$, and $\boldsymbol{v}_{B}(x,y,\tau )=(u_{B},v_{B})$ and $p_{B}(x,\tau )$ describe the buoyancy-induced corrections, which evolve in the long time scale $\tau$ as the solute spreads in the channel.

3.3. Steady streaming

The solution procedure needed to compute the velocity corrections parallels that followed earlier at leading order. Thus, the longitudinal component of the steady-streaming velocity

(3.19)\begin{equation} \frac{u_{SS}}{\alpha^{2}}={-}\frac{{\rm d} p_{SS}}{{\rm d} x} \frac{1}{2} (H-y)y+y\int_0^{y} F \, {\rm d} \hat{y} -\int_0^{y} F \hat{y} \, {\rm d} \hat{y}-y \int_0^{H} F\left(1-\frac{y}{H} \right) {\rm d} {y} \end{equation}

is determined by integrating the momentum equation (3.16) written for ${{\textit {Ri}}}=0$ subject to the boundary conditions $u_{SS}=0$ at $y=0,H$. The result can be substituted into (3.15) to give

(3.20)\begin{align} \frac{v_{SS}}{\alpha^{2}}&=\frac{\partial}{\partial x} \left[\frac{{\rm d} p_{SS}}{{\rm d} x} \frac{y^{2}}{2} \left(\frac{H}{2}-\frac{y}{3}\right)+\frac{y^{2}}{2}\int_y^{H} F \left(1-\frac{\hat{y}}{H} \right) {\rm d} \hat{y} \right.\nonumber\\ &\quad \left.+\,y\left(1-\frac{y}{2H}\right)\int_0^{y} F \hat{y} \, {\rm d} \hat{y}-\frac{1}{2} \int_0^{y} F \hat{y}^{2} \, {\rm d} \hat{y}\right] \end{align}

upon integrating with $v_{SS}=0$ at $y=0$. To determine the unknown pressure gradient ${\rm d} p_{SS}/{\rm d} x$ we begin by using $v_{SS}=0$ at $y=H$ in the above equation to give

(3.21)\begin{equation} \frac{\rm d}{{\rm d} x} \left(\int_0^{H} \frac{u_{SS}}{\alpha^{2}} \,{\rm d} y\right)=\frac{\rm d}{{\rm d} x} \left(\frac{{\rm d} p_{SS}}\,{{\rm d} x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} F y (H-y) \, {\rm d} y \right)=0, \end{equation}

which indicates that the flow rate

(3.22)\begin{equation} Q_{SS}=\int_0^{H} u_{SS} \,{\rm d} y=\alpha^{2} \left[ \frac{{\rm d} p_{SS}}{{\rm d} x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} F y (H-y) \, {\rm d} y\right] \end{equation}

is constant. Its value can be determined by integrating a second time (3.22) subject to $p_{SS}(n)=p_{SS}(0)=0$ to yield

(3.23)\begin{equation} Q_{SS}= \frac{\displaystyle\alpha^{2} \int_0^{1} H^{{-}3} \left[\int_0^{H} F y (H-y) \,{\rm d} y \right] {\rm d} x}{2 \displaystyle\int_0^{1} H^{{-}3} \,{\rm d} x}, \end{equation}

when account is taken of the spatial periodicity of $H$ and $F$. Since $H(x)$ is symmetric about $x=1/2$ while $F(x,y)$ is antisymmetric, the double integral in the numerator of the above equation is identically zero, so that

(3.24)\begin{equation} Q_{SS}=\int_0^{H} u_{SS} \,{\rm d} y= 0, \end{equation}

in agreement with previous findings regarding steady streaming in tubes (Hall Reference Hall1974) and channels (Lo Jacono et al. Reference Lo Jacono, Plouraboué and Bergeon2005; Guibert et al. Reference Guibert, Plouraboué and Bergeon2010). Using the condition $Q_{SS}=0$ in (3.22) finally yields

(3.25)\begin{equation} \frac{{\rm d} p_{SS}}{{\rm d} x}={-}6 H^{{-}3} \int_0^{H} F y (H-y) \, {\rm d} y, \end{equation}

for the pressure gradient, thereby completing the solution.

3.4. Buoyancy-induced velocity

The corresponding solution for the buoyancy-induced velocity can be obtained by simply replacing $F(x,y)$ with ${{\textit {Ri}}} \, c_0(x,y,\tau )$ in (3.19) and (3.20), yielding

(3.26)\begin{equation} \frac{u_{B}}{\alpha^{2} {{\textit{Ri}}}}={-}\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{1}{2} (H-y)y+y\int_0^{y} c_0 \, {\rm d} \hat{y} -\int_0^{y} c_0 \hat{y} \, {\rm d} \hat{y}-y \int_0^{H} c_0\left(1-\frac{y}{H} \right) {\rm d} {y} \end{equation}

and

(3.27)\begin{align} \frac{v_{B}}{\alpha^{2} {{\textit{Ri}}}}&=\frac{\partial}{\partial x} \left[\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{y^{2}}{2} \left(\frac{H}{2}-\frac{y}{3}\right)+\frac{y^{2}}{2}\int_y^{H} c_0 \left(1-\frac{\hat{y}}{H} \right) {\rm d} \hat{y} \right.\nonumber\\ &\quad \left.+\,y\left(1-\frac{y}{2H}\right)\int_0^{y} c_0 \hat{y} \, {\rm d} \hat{y}-\frac{1}{2} \int_0^{y} c_0 \hat{y}^{2} \, {\rm d} \hat{y}\right]. \end{align}

Using the condition $v_{B}=0$ at $y=H$ in the above equation and integrating once gives

(3.28)\begin{equation} \frac{Q_{B}}{\alpha^{2} {{\textit{Ri}}}}=\frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x} \frac{H^{3}}{12}+\frac{1}{2} \int_0^{H} c_0 y (H-y) \, {\rm d} y \end{equation}

for the buoyancy-induced flow rate $Q_{B}=\int _0^{H} u_{B} \,{\rm d} y$. Integrating a second time with $p_{B}=0$ at $x=0,n$ to give

(3.29)\begin{equation} Q_{B}(\tau) =\frac{\displaystyle\alpha^{2} {{\textit{Ri}}} \int_0^{n} H^{{-}3} \left[\int_0^{H} c_0 y (H-y) \,{\rm d} y \right] {\rm d} x}{\displaystyle 2 n \int_0^{1} H^{{-}3}\, {\rm d} x}, \end{equation}

finally determines the pressure gradient

(3.30)\begin{equation} \frac{1}{{{\textit{Ri}}}} \frac{\partial p_{B}}{\partial x}=\frac{6}{H^{3}} \left\{\frac{\displaystyle\int_0^{n} H^{{-}3} \left[\int_0^{H} c_0 y (H-y) \,{\rm d} y \right] {\rm d} x}{\displaystyle n \int_0^{1} H^{{-}3} \,{\rm d} x}-\int_0^{H} c_0 y (H-y) \, {\rm d} y \right\}, \end{equation}

which can be used in (3.26) and (3.27) to complete the determination of the buoyancy-induced velocity. Note that, because $c_0$ is not spatially periodic, the solution carries a dependence on the channel length $n$ through the pressure gradient $\partial p_{B}/\partial x$.

5.1. Buoyancy-free flow

As mentioned previously, in the absence of buoyancy, the flow induced by the imposed pressure gradient is periodic in time and space. The steady Lagrangian motion for $\varepsilon \ll 1$ is given in this case by the sum of the steady-streaming velocity $\boldsymbol{v}_{SS}$ and the Stokes-drift velocity $\boldsymbol{v}_{SD}$. These two contributions as well as their sum are shown in figure 2 for $\beta =0.4$ and two different values of $\alpha$. As the flow in each cell is identical, it suffices to show the solution for $0\le x \le 1$, symmetric with respect to the centre line $x=0.5$. For each value of $\alpha$, streamlines are plotted using a fixed increment $\delta \psi$ of the streamfunction $\psi$, defined in the usual way (e.g. $\partial \psi /\partial y=u_{SS}$ and $\partial \psi /\partial x=-v_{SS}$ for steady streaming) with $\psi =0$ along the wall, so that the interline spacing provides a measure of the local velocity. To further quantify the motion, colour contours are used to represent the associated vorticity $\varOmega =(h_o/\lambda )^{2} \partial v/\partial x-\partial u/\partial y$, which reduces to $\varOmega =-\partial u/\partial y$ in the slender flow approximation.

Figure 2. Streamlines and colour contours of vorticity corresponding to the steady-streaming velocity $\boldsymbol{v}_{SS}$, Stokes-drift velocity $\boldsymbol{v}_{SD}$ and steady mean Lagrangian velocity $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ in a canal with $\beta =0.4$ for (a) $\alpha =4$ and (b) $\alpha =16$. Results of DNS of a non-buoyant flow (i.e. ${{\textit {Ri}}}=0$) with $h_o/\lambda =1/20$ and $\varepsilon =0.02$ are also shown, including the rescaled time-averaged Eulerian velocity field $\langle \boldsymbol{v} \rangle /\varepsilon$ (first column) and the rescaled Lagrangian velocity $\boldsymbol{v}_{L}/\varepsilon$ (fifth column), the latter determined by following tracer particles, as explained in the text. To facilitate the comparisons, fixed constant streamline spacings $\delta \psi =0.003$ and $\delta \psi =0.03$ are used for the upper and lower plots, respectively.

The spatially periodic, time-independent, steady-streaming velocity computed with (3.19) and (3.20) supplemented with (3.25) is shown in the second column of figure 2. The results are qualitatively similar to those presented in Guibert et al. (Reference Guibert, Plouraboué and Bergeon2010). For $\alpha =4$, the flow structure of each half cell exhibits two counter-rotating vortices, whereas for $\alpha =16$ the flow develops an additional, much weaker vortex, located near the section with largest width. As expected, the vorticity, having peak values of order unity for $\alpha =4$, increases with increasing flow frequency as a result of augmented wall production to reach peak values exceeding $\varOmega =40$ for $\alpha =16$.

The steady-streaming results are compared with time-averaged velocity fields obtained in DNS with $\varepsilon =0.02$. In the comparison, the time-averaged DNS velocity is expressed in the rescaled form $\langle \boldsymbol{v} \rangle /\varepsilon \sim 1$, consistent with the scaling employed in defining $\boldsymbol{v}_{SS}$. The two functions $\boldsymbol{v}_{SS}$ and $\langle \boldsymbol{v} \rangle /\varepsilon$ are seen to be almost identical, thereby giving additional confidence in the mathematical development. For instance, the peak values of the stream function and vorticity corresponding to the time-averaged DNS velocity $\langle \boldsymbol{v} \rangle /\varepsilon$ are $\psi _{PEAK}=\pm (0.0115, 0.1680)$ and $\varOmega _{PEAK}=\pm (1.4465, 40.786)$ for $\alpha =(4,16)$, whereas the corresponding values for the steady-streaming motion are $\psi _{PEAK}=\pm (0.0115,0.1699)$ and $\varOmega _{PEAK}=\pm (1.4474,40.787)$. The small relative differences remain below about 1 %, as is consistent with the order of the asymptotic development.

The third column in figure 2 displays the Stokes-drift velocity field evaluated with (4.6a,b). As it is clear from a quantitative comparison with the corresponding steady-streaming results, both bulk-flow velocities have comparable magnitude for $\alpha =4$, whereas for $\alpha =16$ the Stokes drift provides a much smaller relative contribution to the Lagrangian drift. The dominant role of steady streaming in flows at high Womersley numbers is found also away from the wall in oscillating flow over circular cylinders (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954; Raney, Corelli & Westervelt Reference Raney, Corelli and Westervelt1954), for example. The mean Lagrangian velocity $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ corresponding to the asymptotic limit $\varepsilon \ll 1$ compares favourably with the velocity $\boldsymbol{v}_{L}/\varepsilon$ obtained numerically by following tracer particles in the DNS computation for $\varepsilon =0.02$, shown in the last column of figure 2, although the relative errors are somewhat larger than those of the Eulerian velocity. For instance, the peak values of the stream function and vorticity corresponding to $\boldsymbol{v}_{L}/\varepsilon$ are $\psi _{PEAK}=\pm (0.0235, 0.1674)$ and $\varOmega _{PEAK}=\pm (2.0454, 45.9497)$ for $\alpha =(4,16)$, while the corresponding values for $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ are $\psi _{PEAK}=\pm (0.0235, 0.1491)$ and $\varOmega _{PEAK}=\pm (1.9906, 40.787)$.

5.2. Buoyancy-free solute dispersion

The reduced transport (4.8) resulting from the two-time-scale asymptotic analysis indicates that the solute relies on the Lagrangian drift for longitudinal dispersion. As a consequence, the existence of the closed recirculating vortices displayed in figure 2 implies that in oscillatory buoyancy-free channel flow a solute released in a given cell would be unable to reach their neighbouring cells, thereby precluding its progression along the canal. To illustrate this important feature of the flow, we show in figure 3 the temporal evolution of a bolus of solute with reduced Schmidt number $\sigma =1$ released at the initial instant of time in the central cell of a three-cell canal. The initial concentration is given by the truncated Gaussian distribution

(5.1)\begin{equation} c_i = \min\left\{1,\frac{3}{2}\exp\left[{-}16\left(\frac{x-x_0}{\delta}\right)^{2}\right]\right\}, \end{equation}

which represents a band of solute with characteristic width $\delta$ centred at $x_o$ having a saturated core flanked by thin layers across which the concentration decays to zero. Results obtained by integration of (4.8) for $x_0=1.75$ and $\delta =0.2$ are compared in figure 3 at different instants of time in the interval $0\le \tau \le 8$ with DNS computations. Note that, with $\tau =\varepsilon ^{2} t$, for the value $\varepsilon =0.02$ used in the DNS, this interval of time corresponds to $0\le t \le 20\,000$ (i.e. about $20\,000/2{\rm \pi} \simeq 3200$ oscillatory cycles).

Figure 3. Snapshots of solute concentration for ${{\textit {Ri}}}=0$, $\beta =0.2$, $\alpha =4$ and $\sigma =1$ as obtained at five different instants of time from the reduced model and from DNS for $h_o/\lambda =1/20$ and $\varepsilon =0.02$. In addition to colour contours of local concentration, the figure shows distributions of solute per unit channel length $\int _0^{H} c \,{\rm d} y$ for the model (solid curve) and for the DNS (dot-dashed curves), with the initial distribution $\int _0^{H} c_i \,{\rm d} y$ shown as a dotted curve. For reference, the figure shows streamlines with constant spacing $\delta \psi =0.003$ for the Lagrangian mean drift, which is characterised using the asymptotic prediction $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}$ for the model results and the value of $\boldsymbol{v}_{L}/\varepsilon$ determined numerically for the DNS results.

As seen in figure 3 the model accurately reproduces the DNS results. To facilitate the quantitative comparisons, in addition to colour contours showing the solute concentration, the figure includes side plots for the amount of solute per unit channel length $C=\int _0^{H} c \, {\rm d} y$ at different instants of time, with the initial distribution $C_i=H(x) c_i(x)$ included for reference as a dotted curve. The model predictions lie very close to the DNS results, in that the normalised value of the integrated departure $(\int _0^{n} |C_{DN}-C_{MODEL}|\,{\rm d} x)/(\int _0^{n} C_i \,{\rm d} x)$, which provides a metric for the accuracy of the model, remains below $0.003$ over the entire range of times considered in the figure.

For the buoyancy-free conditions considered in figure 3, the steady Lagrangian motion is seen to stir the solute about the deposition location, uniformising its concentration within the recirculating cell. The effect of longitudinal diffusion, present in the DNS results, is found to be rather limited, in that, even at the latest instant of time considered, the presence of the solute in the adjacent cells is negligibly small. This tendency of the solute to remain trapped inside Lagrangian vortices has potential implications concerning the drug-dispersion rate in ITDD procedures. Although the Lagrangian flow in the spinal canal does not exhibit the spatial periodicity of the canonical configuration investigated here, closed recirculating vortices, associated with the changes in the eccentricity of the spinal cord along the canal, have been found to characterise the CSF bulk motion (Coenen et al. Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019). Typically, there are three main vortices, extending along the cervical, thoracic and lumbar regions. As ITDD injection occurs in the lumbar region, the buoyancy-free results in figure 3 seem to indicate that, when the density of the drug matches exactly the CSF density, the drug is bound to linger in the lumbar vortex near the injection site without reaching the thoracic region. This could be advantageous in applications involving pain medication, which is meant to be delivered to the spinal cord, but not in applications involving anticancer drugs targeting brain tumors, for example. As shown in the following, buoyancy-induced motion has the potential to drastically change the associated transport rate, in accordance with clinical observations (Mitchell et al. Reference Mitchell, Bowler, Scott and Edström1988; Povey et al. Reference Povey, Jacobsen and Westergaard-Nielsen1989; Richardson et al. Reference Richardson, Thakur, Abramowicz and Wissler1996; Veering et al. Reference Veering, Immink-Speet, Burm, Stienstra and van Kleef2001).

5.3. Slowly varying buoyancy-induced motion

As reasoned previously, buoyancy forces, acting on solutes with density $\rho _s \ne \rho$, alter the steady Lagrangian drift by adding an additional component that varies slowly in the long time scale $\tau$ following the solute dispersion, so that the flow and the solute transport are intimately coupled, as described by (4.8) supplemented with (3.26), (3.27) and (3.30). The corresponding behaviour is characterised in figure 4 for a light solute with ${{\textit {Ri}}}=1$ spreading upwards. Note that, because of the problem symmetry, results corresponding to a heavy solute with ${{\textit {Ri}}}=-1$ can be generated by simply reversing the direction of the gravity vector, i.e. by rotating the figure $180^{\circ }$.

Figure 4. Snapshots of solute concentration for ${{\textit {Ri}}}=1$, $\beta =0.2$, $\alpha =4$ and $\sigma =1$ as obtained at five different instants of time from the reduced model and from DNS for $h_o/\lambda =1/20$ and $\varepsilon =0.02$. In addition to colour contours of local concentration, the figure shows distributions of solute per unit channel length $\int _0^{H} c \,{\rm d} y$ for the model (solid curve) and for the DNS (dot-dashed curves), with the initial distribution $\int _0^{H} c_i \,{\rm d} y$ shown as a dotted curve. The plots include streamlines with constant spacing $\delta \psi =0.003$ for the varying Lagrangian mean drift, which is evaluated with use of $\boldsymbol{v}_{SS}+\boldsymbol{v}_{SD}+\boldsymbol {v}_{B}$ (model results) and from the displacement of the tracer particles (DNS results).

As in the buoyancy-free flow depicted in figure 3, the solution in figure 4 includes Lagrangian streamlines, colour contours of solute concentration, and streamwise distributions of integrated solute concentration $C=\int _0^{H} c \, {\rm d} y$ along the canal. Buoyancy has a dramatic effect on the dispersion of the solute, as is apparent by comparing the results in both figures. Gravitational forces acting on the light solute induce a longitudinal pressure gradient that modifies drastically the resulting Lagrangian drift, as can be seen by comparing the streamlines in figure 3 with those in figure 4. The pattern of symmetric recirculating cells with unconnected streamlines existing for ${{\textit {Ri}}}=0$ is replaced for ${{\textit {Ri}}}=1$ by a more complicated streamline pattern featuring a net upward flow rate $Q_{B}(\tau )$ (see the solid curves in figure 5, to be discussed later). As can be seen, although the flow rate and the associated streamline pattern vary slowly in time, the observed changes are not very pronounced. Also of interest is that the quantitative agreement between the streamlines predicted by the model and the DNS results is again remarkable, thereby giving additional confidence in our development.

Figure 5. Influence of (a) the contraction ratio $\beta$, (b) reduced Schmidt number $\sigma$, (c) Richardson number ${{\textit {Ri}}}$ and (d) Womersley number $\alpha$ on the temporal evolution of the buoyancy-induced flow rate $Q_{B}$. The values of the parameters in each case are: (a) $\alpha =4$, $\sigma =Ri=1$; (b) $\alpha =4$, $Ri=1$ and $\beta =0.2$; (c) $\alpha =4$, $\sigma =1$ and $\beta =0.2$; (d) $\sigma =Ri=1$ and $\beta =0.2$.

The changes in the Lagrangian motion have a dramatic reflection in the solute dispersion. As shown in figure 4, the solute is transported upwards following the Lagrangian streamlines connecting the cells, enabling its upward progression. The bolus of solute is distorted by the recirculating flow as it travels upwards, driven by the buoyancy-induced draft. The variation with time of the concentration distribution predicted with the model is in excellent agreement with the DNS results. The model is shown to predict not only the mean location of the bolus but also its shape and elongation. The relative departures, measured by $(\int _0^{n} |C_{DNS}-C_{MODEL}|\, {\rm d} x)/(\int _0^{n} C_i \,{\rm d} x)$, remain below $0.009$ over the entire range of times shown in the figure. In assessing the potential benefits of the time-averaged formulation, it is important to emphasise that, although the integration of the integro-differential equation (4.8) over times $\tau \sim 1$ can be completed in a few hours using a laptop computer, generating the DNS results shown in figure 4, spanning over 3000 cardiac cycles, required 10 days in a computational cluster using a total of 72 cores.

5.4. Parametric dependence of the results

The case shown in figure 4, corresponding to $\beta =0.2$, $\alpha =4$, $\sigma =Ri=1$ is used in figures 5 and 6 as a basis to investigate the influence of the different parameters on the dispersion of a buoyant solute. To that end, results are generated with use of (4.8) by modifying one of the four controlling parameters at a time, while keeping the other three fixed at the values selected earlier. Figure 5 shows the variation with time of the buoyancy-induced flow rate $Q_{B}$, whereas figure 6 shows instantaneous solute-concentration distribution and associated Lagrangian streamlines at a fixed time, namely, $\tau =(6,6,6,2)$ for figures 6(a)–6(d), with corresponding results for the base case at these times shown in two of the subpanels of figure 4.

Figure 6. Influence of (a) the contraction ratio $\beta$, (b) reduced Schmidt number $\sigma$, (c) Richardson number ${{\textit {Ri}}}$ and (d) Womersley number $\alpha$ on the solute-concentration distribution and associated Lagrangian streamlines. The snapshots are taken at $\tau =6$ for (a–c) and at $\tau =2$ for (d). The values of the parameters in each case are: (a) $\alpha =4$, $\sigma =Ri=1$; (b) $\alpha =4$, $Ri=1$ and $\beta =0.2$; (c) $\alpha =4$, $\sigma =1$ and $\beta =0.2$; (d) $\sigma =Ri=1$ and $\beta =0.2$. The streamline spacing for the mean Lagrangian velocity is $\delta \psi =0.003$ for (a–c) and $\delta \psi =0.005$ for (d).

We begin by discussing the effect of the channel geometry. As shown in figure 6(a), increasing the undulation of the channel from $\beta =0.1$ to $\beta =0.4$ tends to increase the magnitude of the buoyancy-induced longitudinal velocity in the region of minimum cross-sectional area, where streamlines are closely spaced together for larger $\beta$, but these changes result in only moderately small variations of the flow rate $Q_{B}$, as shown in figure 5(a). As a result, the bolus of solute becomes more elongated as $\beta$ increases, but advances upward at approximately the same rate, so that the maximum concentration occupies approximately the same location at $\tau =6$, as shown in figure 6(a).

The effect of the Schmidt number $\sigma$, entering in the formulation only through the factor affecting the transverse diffusion rate on the right-hand side of (4.8), is investigated in figures 5(b) and 6(b). The changes observed in streamline pattern and flow rate when changing the Schmidt number from $\sigma =0.25$ to $\sigma =8$ are not very significant, so that the differences in solute evolution in figure 6(b) are attributable to the direct effect of transverse diffusion (or lack thereof). The snapshot corresponding to $\sigma =8$ displays the behaviour expected at large Schmidt numbers, for which fluid particles maintain a nearly constant concentration in their slow Lagrangian evolution, as described by the limiting form of (4.8) for $\sigma \gg 1$. In the opposite limit $\sigma \ll 1$, solute diffusion leads to a rapid uniformisation of the concentration, as can be seen by integrating $\partial ^{2} c/\partial y^{2}=0$ (i.e. the reduced form of (4.8) when $\sigma \ll 1$) subject to the non-permeability condition $\partial c/\partial y=0$ at $y=0,H$ to give $c=c(x,\tau )$. As a result, the bolus remains relatively compact as it moves along the channel with a velocity determined by the flow rate. The reduced transport equation governing the transport of solutes with $\sigma \ll 1$ can be derived from (4.9) by noting that in this limit the integrated solute concentration $C=\int _0^{H} c \, {\rm d} y$ becomes $C=H(x) c(x,\tau )$ while the solute flux $\phi (x,\tau )=\int _0^{H} (u_{SD}+u_{SS}+u_{B})c \, {\rm d} y$ reduces to $\phi =Q_{B} c$. Substituting these simplified expressions into (4.9) and using (3.29) to evaluate $Q_{B}$ finally yields

(5.2)\begin{equation} \frac{\partial c}{\partial \tau}+\frac{\alpha^{2} {{\textit{Ri}}} \displaystyle\int_0^{n} c \,{\rm d} x}{12 n H \displaystyle\int_0^{1} H^{{-}3} \,{\rm d} x} \frac{\partial c}{\partial x}=0 \end{equation}

as the limiting form of (4.8) for $\sigma \ll 1$.

As is to be expected from the velocity expressions derived in § 3.4, the Richardson number and the Womersley number have a pronounced effect on the mean Lagrangian motion. As shown in figures 5(c) and 5(d), the flow rate exhibits dependences on $\textit{Ri}$ and $\alpha$ that are approximately linear and approximately quadratic, respectively, consistent with (3.29). These dependences have a reflection on the evolution of the solute bolus shown in figures 6(c) and 6(d). With limited updraft for ${{\textit {Ri}}}=0.25$, the bolus is seen to spread about the injection location, without significant upward progression at the instant of time $\tau =6$ considered in the figure. An increase in ${{\textit {Ri}}}$ promotes the displacement of the bolus, but its longitudinal extent remains approximately equal in all three cases. By way of contrast, an increase in $\alpha$ increases the upward displacement and also enhances bolus distortion. The reason for the latter is that larger values of $\alpha$ hinder transverse diffusion, as can be inferred from (4.8), with the result that fluid particles travel following the Lagrangian recirculating paths with a nearly constant concentration, rapidly deforming the compact concentration distribution of the initial bolus.