2. Formulation

Let us consider the flow configuration depicted in figure 1, emerging when an incompressible fluid stream with harmonic velocity $U_\infty \cos (\omega t')$ flows past an infinite array of equally spaced identical cylinders aligned with the unperturbed flow. The semi-distance $L$ between the centres of contiguous cylinders is assumed to be comparable to the cylinders radius $a$, their ratio defining the geometrical parameter $\ell =L/a \geqslant 1$. As previously anticipated, the two controlling flow parameters are the dimensionless stroke length $\varepsilon$, defined in (1.1), and the Womersley number $M$, defined in (1.2). DNS corresponding to order-unity values of the three parameters $\ell$, $M$ and $\varepsilon$ are to be presented below along with results corresponding to the small-stroke-length asymptotic limit $\varepsilon \ll 1$, when the velocity displays a harmonic temporal dependence at leading order, while the first-order corrections, of order $\varepsilon U_\infty$, contain a steady contribution.

The problem is scaled with use of $a$, $\omega ^{-1}$, $U_\infty$ and $\rho \omega a U_\infty$ as characteristic values of length, time, velocity and spatial pressure difference, with $\rho$ denoting the fluid density. Correspondingly, the unperturbed flow velocity of the external oscillating stream becomes $u_\infty =\cos {t}$ with $t=\omega t'$. Since the resulting velocity $\boldsymbol {v}$ is periodic in the streamwise direction, the solution can be described by considering the flow about an individual cylinder, with the origin of the coordinate system placed at the cylinder centre. The description employs Cartesian coordinates $\boldsymbol {x}=(x,y)$ and Cartesian velocity components $\boldsymbol {v}=(u,v)$, with $x$ aligned in the direction of the unperturbed flow and $r=(x^2+y^2)^{1/2}$ denoting the distance to the cylinder centre, as indicated in figure 1. Since, in the regime $\varepsilon \lesssim 1$ and $M\sim 1$ investigated below, the flow can be anticipated to remain symmetric about the $y=0$ plane, in the computations it suffices to consider the integration domain extending for $x^2+y^2>1$ with $y>0$ and $-\ell < x < \ell$. The velocity must satisfy the continuity and momentum equations

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

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{v}}{\partial t}+\varepsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} =-\boldsymbol{\nabla} p +\frac{1}{M^2} {\nabla}^2 \boldsymbol{v}, \end{gather}$$

subject to the non-slip condition

(2.3)\begin{equation} \boldsymbol{v} = 0 \quad {\rm at} \ r=1, \end{equation}

the far-field condition

(2.4)\begin{equation} \boldsymbol{v} = (\cos t, 0) \quad {\rm as} \ y \rightarrow \infty \quad {\rm for} -\ell \leqslant x \leqslant \ell, \end{equation}

the centreline symmetry condition

(2.5)\begin{equation} \frac{\partial u}{\partial y}=v=0 \quad {\rm at} \ y=0 \quad {\rm for} \ 1 \leqslant |x| \leqslant \ell \end{equation}

and the condition of $2\ell$ spatial periodicity in the $x$ direction. The free-stream velocity condition (2.4) is consistent with a far-field pressure distribution approaching $p=x \sin t$ as $y \rightarrow \infty$.

The above problem was integrated numerically using the immersed boundary method with the projection approach given by Taira & Colonius (Reference Taira and Colonius2007) in a Cartesian non-uniform staggered mesh extending up to $y=30$. The value of the associated grid spacing $\varDelta$, smaller near the cylinder surface, was reduced for increasing values of the Womersley number as needed to resolve the associated near-wall Stokes layer with sufficient accuracy, yielding, for instance, $\varDelta =0.04$ for $M=1$ and $\varDelta =0.01$ for $M=16$. The spatial width of the cylinder nodes employed in the implementation of the immersed boundary method was selected to be equal to the smallest spacing of the Cartesian mesh. The time step $\delta t$ was correspondingly adjusted to give a Courant number $\delta t/\varDelta$ below 0.25. Following standard practice (see e.g. Alaminos-Quesada Reference Alaminos-Quesada2021), the implementation of the far-field condition (2.4) was facilitated in the simulations by rewriting the problem in terms of the axial velocity perturbation $u'=u-\cos t$, which satisfies $u'=-\cos t$ at $r=1$ and $u' \rightarrow 0$ as $y \rightarrow \infty$ along with the symmetry and periodicity conditions stated above. As explained in Appendix A, the numerical method was validated through comparisons with previously reported results corresponding to a single cylinder.

3. The limit of small stroke lengths

Following standard practice, the flow description in the limit $\varepsilon \ll 1$ utilizes expansions for the different flow variables in powers of $\varepsilon$, i.e. $\boldsymbol {v}=\boldsymbol {v}_0+\varepsilon \boldsymbol {v}_1 +\cdots$ and $p=p_0+\varepsilon p_1 +\cdots$. As seen below, the leading-order solution has a zero time average, i.e. $\langle \boldsymbol {v}_0 \rangle =0$, with $\langle {\cdot } \rangle = ({1}/{2 {\rm \pi}} )\int _t^{t+2 {\rm \pi}} {\cdot } \, {\rm d}t$, whereas the first-order correction $\boldsymbol {v}_1$, accounting for the effects of convective acceleration, includes a non-zero steady-streaming component $\boldsymbol {v}_{SS}=\langle \boldsymbol {v}_1 \rangle$.

3.1. Leading-order oscillatory flow

At leading order in the limit $\varepsilon \ll 1$, convective acceleration does not enter in the momentum balance equation (2.2). The resulting linear problem can be conveniently solved by introducing $\boldsymbol {v}_0={\rm Re} (\mathrm {e}^{\mathrm {i} t} \boldsymbol {V})$ and $p_0={\rm Re} (\mathrm {e}^{\mathrm {i} t} P)$ with $\boldsymbol {V}(x,y)=(U,V)$ and $P(x,y)$ representing complex functions satisfying

(3.1a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{V}=0, \quad \mathrm{i} \boldsymbol{V}=-\boldsymbol{\nabla} P +\frac{1}{M^2} {\nabla}^2 \boldsymbol{V}, \end{equation}

with boundary conditions

(3.2) \begin{equation} \left. \begin{array}{lll@{}} \boldsymbol{V}= 0 & {\rm at} \ r=1, & \\ \boldsymbol{V} = (1, 0) & {\rm as} \ y \rightarrow \infty & {\rm for} \ -\ell \leqslant x \leqslant \ell, \\ \partial U/\partial y=V=0 & {\rm at} \ y=0 & {\rm for} \ 1 \leqslant |x| \leqslant \ell, \end{array} \right\} \end{equation}

as follows from (2.1)–(2.5), along with the condition of $2\ell$ spatial periodicity in the $x$ direction.

Except for the limiting case $\ell \gg 1$, which reduces to that of flow over a single cylinder (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954; Chong et al. Reference Chong, Kelly, Smith and Eldredge2013), no analytic solution is available, and the above problem must be solved numerically. To that aim, (3.1a,b) were written in weak form and implemented in the finite element solver COMSOL Multiphysics. Solutions were computed on an unstructured triangular mesh that extended laterally to $y = 30$. Mesh elements were clustered near the cylinder surface, the typical element size ranging from 0.01 at that surface to 0.2 near the far-field boundary. It was checked that further increases in lateral domain extension as well as in mesh refinement did not alter the results.

For a general value of $M$, the resulting complex velocity $\boldsymbol {V}(x,y)$ has real and imaginary parts. Note, however, that, in the inviscid limit $M \gg 1$, the solution contains an imaginary part only in the thin Stokes layer of thickness $1/M$ that develops on the cylinder surface, outside of which the flow is irrotational, such that $\boldsymbol {V}(x,y)=\boldsymbol {\nabla } \varPhi$. The associated velocity potential $\varPhi$, a real function, satisfies ${\nabla }^2 \varPhi =0$ subject to ideal-flow boundary conditions stemming from (3.2), including, for instance, the no-penetration condition $\partial \varPhi /\partial r=0$ at $r=1$. The problem was considered recently by Crowdy (Reference Crowdy2016), who provided a quasi-analytical solution for the corresponding complex potential. For illustrative purposes, the streamlines of the potential flow corresponding to the specific case $\ell =2$ are included in the schematic of figure 1.

3.2. Steady streaming

The steady-streaming velocity $\boldsymbol {v}_{SS}=\langle \boldsymbol {v}_1 \rangle =(u_{SS},v_{SS})$ is determined from the problem that arises at the following order. Collecting terms of order $\varepsilon$ in (2.1) and (2.2) and taking the time average leads to

(3.3a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}_{{SS}}=0, \quad \frac{1}{2}\,{\rm Re}(\boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{V}^*) =-\boldsymbol{\nabla} \langle \,p_1 \rangle +\frac{1}{M^2} {\nabla}^2 \boldsymbol{v}_{{SS}}, \end{equation}

after writing $\langle \boldsymbol {v}_0 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_0\rangle =\frac {1}{2}\ \textrm {Re}(\boldsymbol {V} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {V}^* )$, which follows from the identity

(3.4)\begin{equation} \langle {\rm Re}(\mathrm{e}^{\mathrm{i} t} A) \, {\rm Re}(\mathrm{e}^{\mathrm{i} t} B) \rangle = {\rm Re}(A B^*)/2, \end{equation}

applying to any generic time-independent complex functions $A$ and $B$, with the asterisk $^*$ denoting complex conjugates. The resulting recirculating cells, symmetric about the $x=0$ plane, can be correspondingly obtained by integrating (3.3a,b) in the first quadrant subject to the boundary conditions

(3.5) \begin{equation} \left. \begin{array}{lll@{}} \boldsymbol{v}_{{{SS}}}= 0 & {\rm at} \ r=1, & \\ \boldsymbol{v}_{{{SS}}} \rightarrow 0 & {\rm as} \ y \rightarrow \infty & {\rm for} \ 0 \leqslant x \leqslant \ell,\\ {\partial u_{{SS}}}/{\partial y}=v_{{{SS}}}=0 & {\rm at} \ y=0 & {\rm for} 1 \leqslant x \leqslant \ell, \end{array} \right\} \end{equation}

consistent with (2.3)–(2.5), and the condition of $2\ell$ spatial periodicity in the $x$ direction. At this order, the steady-streaming pressure $\langle \,p_1 \rangle$ vanishes in the far field, as is consistent with the velocity condition $\boldsymbol {v}_{SS} \rightarrow 0$ as $y \rightarrow \infty$.

Equations (3.3a,b) were integrated using the numerical method employed earlier for the leading-order problem. Representative results are shown in figure 2 for four values of the inter-cylinder spacing $\ell$, including as extreme cases the configuration with touching cylinders ($\ell =1$) and the familiar single-cylinder case, recovered in the present array configuration when $\ell =\infty$. Because of the condition of flow periodicity and the symmetry of the cylinder array, the vertical lines $x=0, 1 \leqslant y < \infty$ and $x=\ell, 0 \leqslant y < \infty$ are streamlines of the steady-streaming flow. Only the first quadrant is shown in figure 2, since the flow structure is identical in all four quadrants. Streamlines are plotted using a fixed increment $\delta \psi$ of the streamfunction $\psi _{SS}$ computed from $\partial \psi _{SS}/\partial y=u_{SS}$ and $\partial \psi _{SS}/\partial x=-v_{SS}$, with $\psi _{SS}=0$ on the domain boundary. The spacing is $\delta \psi =0.005$ for $\ell =1.5$ and $\ell =3.0$, with a smaller spacing $\delta \psi =0.002$ used for $\ell =1$, as needed to represent the associated weak motion, and a larger spacing $\delta \psi =0.01$ for $\ell =\infty$, in accordance with the associated vigorous motion. In addition to streamlines, colour contours are used to represent the vorticity $\varOmega =\partial v/\partial x-\partial u/\partial y$, with the level indicated in the colour bar on the far right.

Figure 2. Streamlines and colour contours of vorticity $\varOmega$ corresponding to the steady-streaming motion with different inter-cylinder distance $\ell$ for $M=2$ (a) and $M=16$ (b). Streamlines are represented using a constant spacing $\delta \psi$, with $\delta \psi =0.002$ for $\ell =1$, $\delta \psi =0.005$ for $\ell =1.5$ and $3$, and $\delta \psi =0.01$ for $\ell =\infty$. Corresponding vorticity levels are indicated in the colour bar on the right.

As seen in figure 2, the streaming structure arising for finite values of $\ell$ is qualitatively similar to that of a single cylinder (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954). For $M=2$ the flow displays one vortex in each quadrant, with the clockwise circulation (negative vorticity) exhibited by the vortex in the first quadrant corresponding to fluid approaching the cylinder along the oscillation axis $y=0$. This vortex is known to progressively approach the cylinder wall on increasing $M$ (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954) and, for the case $M=16$ shown in figure 2(b), is seen to be embedded in the high-vorticity Stokes layer that develops near the cylinder surface. A second vortex with opposite circulation, clearly visible in the results for $M=16$, appears outside in configurations with $M$ exceeding a critical value $M_c$. For the case of a single cylinder, the value $M_c \simeq 6.08$ can be determined from the exact solution (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954) as the value of $M$ for which the streamfunction $\psi _{SS}$ vanishes in the far field. From our numerical computations, it was seen that the value of $M_c$ is somewhat larger for the cylinder array (e.g. $M_c \simeq 7$ for $\ell =2$).

The presence of the neighbouring cylinders has a noticeable effect on the shape of the resulting vortices, as can be seen by comparing the results for $\ell =(1,1.5,3)$ with the canonical case of a single cylinder ($\ell =\infty$) shown in the last column of figure 2. For $M=2$ the core of the vortex, which for $\ell =\infty$ is located along the ${\rm \pi} /4$ ray, is displaced towards the vertical axis $x=0$ on decreasing the inter-cylinder spacing, producing vortices that are much more slender, with the case $\ell =1$ displaying the largest distortion. For $M=16$, the outer vortex, which for the single cylinder exhibits open streamlines with no vortex core, displays for $\ell \ne \infty$ a well-defined core surrounded by closed streamlines. This qualitative change, also observed in the flow about an oscillating cylinder when enclosed by a concentric cylindrical surface (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954), is attributable to the effect of confinement, which also produces a drastic reduction in the magnitude of the streaming motion. The extent of the reduction can be quantified by comparing the peak value of the streamfunction, given by $\psi _{{SS},{peak}}=-0.1602$ for $M=2$ and $\psi _{{SS},{peak}}=(-0.0493/0.243)$ (inner/outer vortex) for $M=16$ in the case of the isolated cylinder ($\ell =\infty$) and $\psi _{{SS},{peak}}=-0.0041$ for $M=2$ and by $\psi _{{SS},{peak}}=(-0.0438/0.0022)$ (inner/outer vortex) for $M=16$ in the case of an array of touching cylinders ($\ell =1$).

3.3. Mean Eulerian velocity for finite stroke lengths

The steady-streaming velocity $\boldsymbol {v}_{SS}=\langle \boldsymbol {v}_1 \rangle$ provides the leading-order description for the mean Eulerian velocity $\langle \boldsymbol {v} \rangle =\varepsilon \boldsymbol {v}_{SS}$ in the asymptotic limit $\varepsilon \ll 1$. In principle, the description can be improved by computing higher-order terms in the asymptotic expansion for $\langle \boldsymbol {v} \rangle = \varepsilon \langle \boldsymbol {v}_1 \rangle + \varepsilon ^2 \langle \boldsymbol {v}_2 \rangle +\varepsilon ^3 \langle \boldsymbol {v}_3 \rangle +\cdots$. The development must begin by computing the unsteady component of the first-order velocity correction $\boldsymbol {v}_1$, which can be shown to be of the form $\boldsymbol {v}_1-\langle \boldsymbol {v}_1 \rangle =\textrm {Re} (\mathrm {e}^{2 \mathrm {i} t} \boldsymbol {V}_1 )$, where $\boldsymbol {V}_1(x,r)$ is a complex function, the expression of which was obtained by Chong et al. (Reference Chong, Kelly, Smith and Eldredge2013) for the case of a single isolated cylinder. The equations that determine $\langle \boldsymbol {v}_2 \rangle$, analogous to (3.3a,b), with the convective term in the momentum equation replaced by $\langle \boldsymbol {v}_0 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_1 \rangle +\langle \boldsymbol {v}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_0\rangle$, are to be integrated with the homogeneous boundary conditions stated in (3.5), with $\langle \boldsymbol {v}_2 \rangle$ replacing $\boldsymbol {v}_{SS}$. Since $\boldsymbol {v}_0=\textrm {Re} (\mathrm {e}^{\mathrm {i} t} \boldsymbol {V})$ and $\boldsymbol {v}_1=\langle \boldsymbol {v}_1 \rangle +\textrm {Re} (\mathrm {e}^{2 \mathrm {i} t} \boldsymbol {V}_1 )$, it follows that $\langle \boldsymbol {v}_0 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_1 \rangle +\langle \boldsymbol {v}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_0\rangle =0$, with the consequence that integration of the steady-streaming problem that arises at order $\varepsilon ^2$ yields $\langle \boldsymbol {v}_2 \rangle =0$. Therefore, the corrections to the mean Eulerian velocity would enter only at the following order, i.e. $\langle \boldsymbol {v} \rangle = \varepsilon \langle \boldsymbol {v}_1 \rangle +\varepsilon ^3 \langle \boldsymbol {v}_3 \rangle +\cdots$, indicating that the leading-order expression $\langle \boldsymbol {v} \rangle = \varepsilon \boldsymbol {v}_{SS}= \varepsilon \langle \boldsymbol {v}_1 \rangle$ computed here contains small relative errors of order $\varepsilon ^2$.

The accuracy of the asymptotic description $\langle \boldsymbol {v} \rangle = \varepsilon \boldsymbol {v}_{SS}$ was tested through comparisons with the mean Eulerian velocity $\langle \boldsymbol {v} \rangle =({1}/{2 {\rm \pi}} )\int _t^{t+2 {\rm \pi}} \boldsymbol {v} \, \textrm {d}t$ determined in direct integrations of the complete problem (2.1)–(2.5). Selected numerical results corresponding to $\ell =2$ and $M=2$ are shown in figure 3($b$–$e$) for $\varepsilon =(0.1,0.5,1.0,2.0)$. Since the time-averaged velocity can be anticipated to be of order $\varepsilon$, as suggested by the asymptotic analysis for $\varepsilon \ll 1$, the rescaled velocity $\langle \boldsymbol {v} \rangle /\varepsilon$ is used in computing the streamlines and vorticity contours shown in figure 3. The results are to be compared with those of the steady-streaming velocity $\boldsymbol {v}_{SS}$, shown in figure 3($a$). Close agreement is found between the DNS results for $\varepsilon =0.1$ and the $\varepsilon \ll 1$ predictions, with the associated velocity fields being nearly identical, as seen in figure 3. A quantitative measure of the existing differences, of the order of 1 % for $\varepsilon =0.1$, consistent with the relative errors of order $\varepsilon ^2$ anticipated in the discussion of the preceding paragraph, is provided by the peak values of the corresponding streamfunctions at the vortex centre, given by $\psi _{SS}=-0.0419$ for $\varepsilon \ll 1$ and $\langle \psi \rangle /\varepsilon =-0.0416$ for $\varepsilon =0.1$. It is remarkable that, although larger differences are found as the oscillation amplitude becomes comparable to the cylinder radius, the $\varepsilon \ll 1$ description remains reasonably accurate even for $\varepsilon =0.5$, for which $\langle \psi \rangle /\varepsilon =-0.0390$ at the vortex centre. For completeness, a figure showing the spatial distribution of $|\psi _{SS}-\langle \psi \rangle /\varepsilon |$ is included in Appendix B.

Figure 3. Streamlines and colour contours of vorticity $\varOmega$ for $\ell =2$ and $M=2$. Besides results corresponding to the steady-streaming velocity $\boldsymbol {v}_{SS}$, shown in panel ($a$), results are given for the rescaled time-averaged Eulerian velocity $\langle \boldsymbol {v} \rangle /\varepsilon$ determined in the DNS computations for $\varepsilon =(0.1,0.5,1.0,2.0)$ in panels ($b$)–($e$). Streamlines are represented using a constant spacing $\delta \psi =0.005$. Corresponding vorticity levels are indicated in the colour bar on the right.

3.4. Stokes drift

As pointed out by Raney et al. (Reference Raney, Corelli and Westervelt1954) when addressing oscillating flow over a cylinder, the Lagrangian mean motion of the fluid particles comes partly from the Eulerian mean motion (i.e. $\langle \boldsymbol {v} \rangle =\varepsilon \boldsymbol {v}_{SS}$) and partly from the so-called Stokes drift (Stokes Reference Stokes1847), a purely kinematic effect arising in non-uniform oscillating flows. As a result, streamlines visualized in experiments by tracing the motion of dyed fluid do not coincide in general with those determined from the steady-streaming velocity (Raney et al. Reference Raney, Corelli and Westervelt1954; Larrieu, Hinch & Charru Reference Larrieu, Hinch and Charru2009; Chong et al. Reference Chong, Kelly, Smith and Eldredge2013). Since the velocity of the Lagrangian mean motion $\boldsymbol {v}_{SS}+\boldsymbol {v}_{SD}$, where

(3.6)\begin{equation} \boldsymbol{v}_{{SD}}=\left\langle \int \boldsymbol{v}_0 \,{\rm d}t \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_0 \right\rangle \end{equation}

represents the contribution of the Stokes drift, determines the convective transport of solutes, there is interest in quantifying numerically $\boldsymbol {v}_{SD}$ for the cylinder array, thereby complementing the analytical results developed previously for the single cylinder (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954; Raney et al. Reference Raney, Corelli and Westervelt1954; Chong et al. Reference Chong, Kelly, Smith and Eldredge2013).

The expression (3.6) for the Stokes-drift velocity, which can be systematically derived using a two-time-scale analysis, as shown in Appendix C, can be written in the form

(3.7)\begin{equation} \boldsymbol{v}_{{SD}}=\tfrac{1}{2}\, {\rm Im}(\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{V}^*), \end{equation}

by using $\boldsymbol {v}_0=\textrm {Re} (\mathrm {e}^{\mathrm {i} t} \boldsymbol {V})$ along with the identity $\langle \textrm {Re}(\mathrm {i} \mathrm {e}^{\mathrm {i} t} A) \, \textrm {Re}(\mathrm {e}^{\mathrm {i} t} B) \rangle =-\textrm {Im}(A B^*)/2$. It is of interest that the real part of the complex function $\tfrac {1}{2} \boldsymbol {V} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {V}^*$ determines the steady streaming, as revealed by (3.3a,b), whereas its imaginary part is the Stokes-drift velocity (3.7). Note that, as mentioned before, for large values of $M$ viscous forces are confined to a thin Stokes layer, outside of which the flow is potential and the function $\boldsymbol {V}$ is real, so that the associated Stokes drift can be expected to vanish for $M \gg 1$, as follows from (3.7).

3.5. Evaluation of the Lagrangian mean velocity

The expression (3.7) was used to evaluate the Stokes-drift velocity $\boldsymbol {v}_{SD}$ for a cylinder array with $\ell =2$, with associated streamlines and vorticity contours given in the middle column of figure 4. The first two columns of figure 4 show the corresponding steady-streaming velocity $\boldsymbol {v}_{SS}$ (second column from the left) along with the rescaled time-averaged Eulerian velocity $\langle \boldsymbol {v} \rangle /\varepsilon$ determined in DNS computations with $\varepsilon =0.1$ (leftmost column), the two sets of results being nearly indistinguishable. Besides the two Womersley numbers $M=2$ and $M=16$ considered earlier in the computations of figure 2, figure 4 includes results for $M=1$, a case for which the Stokes drift is stronger than the steady-streaming motion. To facilitate comparisons, in plotting the streamlines for each value of $M$, the spacing of the Stokes-drift streamfunction $\psi _{SD}$ is that used for the corresponding steady-streaming plot.

Figure 4. Streamlines and colour contours of vorticity $\varOmega$ corresponding to the steady-streaming velocity $\boldsymbol {v}_{SS}$, Stokes-drift velocity $\boldsymbol {v}_{SD}$ and steady mean Lagrangian velocity $\boldsymbol {v}_{L}=\boldsymbol {v}_{SS}+\boldsymbol {v}_{SD}$ for $\ell =2$ and $M=1$ (a), $M=2$ (b) and $M=16$ (c). Corresponding DNS results for $\varepsilon =0.1$ 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). For each value of $M$, streamlines are represented using a constant spacing $\delta \psi =0.002$ ($M=1$) and $\delta \psi =0.005$ ($M=2$ and $M=16$), with the corresponding vorticity levels indicated in the colour bar on the right.

As can be seen, the Stokes-drift results for $M=1$ display a primary clockwise-rotating vortex occupying most of the quadrant, along with a much weaker counter-rotating vortex of negligibly small circulation near the oscillation axis $y=0$. For this value of $M$, this primary vortex is significantly stronger than the corresponding steady-streaming vortex. This can be verified by comparing the magnitude $|\psi _{peak}|$ of the peak values of the associated streamfunctions at the vortex centre. Since $\psi$ is defined to be zero on the cylinder surface, the value of $|\psi _{peak}|$, whose variation with $M$ is represented in figure 5, gives a measure of the volume flow rate driven by the recirculating vortex motion. As can be seen, for $M=1$ the peak value of $\psi _{SD}$ is significantly larger than that of $\psi _{SS}$, with the result that the Lagrangian velocity $\boldsymbol {v}_{SS}+\boldsymbol {v}_{SD}$ is largely determined by its Stokes-drift component, as reflected in the shape of the corresponding Lagrangian vortex, shown in the fourth column of figure 4(a).

Figure 5. The variation with $M$ of the magnitude $|\psi _{peak}|$ of the local peak values of the streamfunction $\psi _{SS}$ (dashed curves), $\psi _{SD}$ (dotted curves) and $\psi _{SS}+\psi _{SD}$ (solid curves) at the centre of the outer ($o$) and inner ($i$) vortices for the $\ell =2$ configuration.

The Stokes-drift motion develops an additional vortex, external to the primary vortex, when the Womersley number is increased to values exceeding a critical value (e.g. $M\simeq 1.5$ for $\ell =2$). As seen in the plots of peak streamfunction in figure 5, this external Stokes-drift vortex, clearly visible in figure 4(b), increases in strength for increasing $M$ to prevail over the inner vortex for $M\gtrsim 2.5$. Figure 5 also reveals that, for the cases $M=2$ and $M=16$ of figures 4(b) and 4(c), the Stokes drift is significantly weaker than the steady streaming, so that the Lagrangian motion is largely determined by the latter.

Figure 5 also shows the peak value of the streamfunction $\psi _{SS}+\psi _{SD}$ associated with the Lagrangian motion. Regarding the resulting curve, it is of interest that, since the inner and outer vortices have opposite circulation, leading to peak values of the streamfunction with different sign, there is an intermediate range of values of $M$ for which the strength of the Lagrangian vortex is smaller than that of the steady-streaming vortex. The comparison of the different curves in figure 5 reveals that the Stokes drift prevails for sufficiently small values of the Womersley number $M \ll 1$, for which $\psi _{SS} \ll \psi _{SD}$, whereas in the opposite limit $M \gg 1$ the Stokes-drift motion fades away, as anticipated above, below (3.7), so that $\psi _{SS} \gg \psi _{SD}$. The trends identified in figure 5 therefore confirm that the Stokes drift can be neglected only if $M \gg 1$, whereas for $M \lesssim 1$ it must be necessarily accounted for when seeking an accurate description of the Lagrangian motion, in agreement with previous findings (Raney et al. Reference Raney, Corelli and Westervelt1954; Chong et al. Reference Chong, Kelly, Smith and Eldredge2013).

To validate the asymptotic prediction $\boldsymbol {v}_{SS}+\boldsymbol {v}_{SD}$, the Lagrangian velocity $\boldsymbol {v}_{L}$ was evaluated from the DNS velocity field for $\varepsilon =0.1$. The value of $\boldsymbol {v}_{L}(x,y)$ at each location $(x,y)$ was determined by computing the displacement $(\delta x,\delta y)$ of a tracer particle, located initially at $(x,y)$, over a cycle (i.e. from $t$ to $t + 2{\rm \pi}$), and the resulting velocity $\boldsymbol {v}_{L}(x,y)=(\delta x,\delta y)/(2 {\rm \pi})$, appropriately rescaled according to $\boldsymbol {v}_{L}/\varepsilon$, was then used to compute the streamlines and vorticity distributions shown in the last column of figure 4, to be compared with the asymptotic predictions shown in the adjacent column. As can be seen, the results are practically indistinguishable, especially for the cases $M=1$ and $M=2$, thereby giving additional confidence in the mathematical development. The somewhat larger departures found with $M=16$, characterized by relative differences in peak streamfunction in the inner and outer vortices of the order of 5 %, are to be expected, since, for these values of $\varepsilon =0.1$ and $M=16$, the relative ordering of the asymptotic development breaks down, in that the viscous term in (2.2) becomes smaller than the convective term. The quantification of these large-Womersley-number configurations can benefit from consideration of the double distinguished limit $\varepsilon \ll 1$ and $M \gg 1$ with ${\textit {Re}}_s=\varepsilon ^2 M^2 \sim 1$ proposed in the seminal analyses of Stuart (Reference Stuart1963, Reference Stuart1966) and Riley (Reference Riley1965, Reference Riley1967).

4. Fluid-particle drift for finite stroke lengths

The above velocity description, in which the Lagrangian mean motion is the result of the superposition of the steady-streaming and Stokes-drift velocity fields, is rigorously valid only in configurations with small stroke lengths $\varepsilon \ll 1$, with representative results presented earlier for $\varepsilon =0.1$ in figure 4. There is interest in testing the accuracy with which the asymptotic prediction $\boldsymbol {v}_{SS}+\boldsymbol {v}_{SD}$ describes the fluid-particle drift as the stroke length $\varepsilon$ increases to larger values. To that end, we computed numerically the trajectories of fluid particles undergoing multiple oscillatory cycles by integrating

(4.1)\begin{equation} \frac{{\rm d}\kern0.7pt \boldsymbol{x}_p}{{\rm d} t}=\varepsilon \boldsymbol{v}(\boldsymbol{x}_p,t), \end{equation}

subject to the initial condition $\boldsymbol {x}_p=\boldsymbol {x}_i$ at $t=t_i$, where $\boldsymbol {x}_p(t)$ represents the fluid-particle location scaled with $a$. The integrations employed the periodic Eulerian velocity $\boldsymbol {v}(\boldsymbol {x},t)$ obtained in DNS computations of pulsating flows with moderate stroke lengths $\varepsilon \sim 1$. Clearly, for a given initial location $\boldsymbol {x}_i$, the resulting trajectory $\boldsymbol {x}_p(t)$ depends on the specific selection of initial time $t=t_i$, so that some care must be taken when defining the mean Lagrangian drift when $\varepsilon$ is not small, as explained below. For a general discussion of Lagrangian mean flow pertaining to nonlinear waves, the reader is referred to the seminal paper of Andrews & McIntyre (Reference Andrews and McIntyre1978).

To illustrate the complications arising in defining the mean Lagrangian drift when $\varepsilon \sim 1$, we plot in figure 6 the results of a representative trajectory calculation, corresponding to oscillatory flow with $M=2$ and $\varepsilon =1$ about a cylinder array with $\ell =2$. Figure 6 shows the path followed by a fluid particle located at $\boldsymbol {x}_i=(-0.55,2.95)$ at $t=t_i={\rm \pi} /2$, corresponding to the instant of time when the outer velocity $u_\infty =\cos t$, decreasing, reaches a zero velocity $u_\infty =0$. For illustrative purposes, stars are used to mark the initial location $\boldsymbol {x}_i$ (red star) as well as the location $\boldsymbol {x}=(-2.62,2.87)$ (blue star) occupied by the fluid particle at time $t=3 {\rm \pi}/2$, when the outer velocity, now increasing from negative values, becomes zero again. The drift motion follows a recirculatory pattern, so that, after a large number of cycles, which would be proportional to $\varepsilon ^{-1}$ for $\varepsilon \ll 1$, the fluid particle returns to occupy a location close to (but not necessarily equal to) the initial location $\boldsymbol {x}_i$.

Figure 6. The grey curves represent the oscillatory trajectories determined numerically by integration of (4.1) with initial condition $\boldsymbol {x}_i=(-0.55,2.95)$ (marked with a red star) at $t_i={\rm \pi} /2$ for $M=2$, $\ell =2$ and $\varepsilon =1.0$. The blue star denotes the particle location at $t=3 {\rm \pi}/2$. The squares mark the fluid-particle locations $\boldsymbol {x}_p({\rm \pi} /2+2 {\rm \pi}n)$ (red squares) and $\boldsymbol {x}_p(3 {\rm \pi}/2+2 {\rm \pi}n)$ (blue squares) for $n=1,2,\ldots$, while the circles are the time averages evaluated with use of (4.2) for $t_i={\rm \pi} /2$ (red circles) and $t_i=3{\rm \pi} /2$ (blue circles).

Different options are available regarding the characterization of the Lagrangian drift. One could, for instance, consider the series of locations $\boldsymbol {x}_n=\boldsymbol {x}_p(t_i+2 {\rm \pi}n)$ occupied by the fluid particle at the end of subsequent cycles $n=1,2,\ldots$. This series, marked in figure 6 by red squares, serves to delineate the long-time drifting motion of the particle as it circles back towards its initial location following a large number of cycles. One can readily see a problem with this definition, in that, if we had considered instead the fluid particle located at $\boldsymbol {x}_i=(-2.62,2.87)$ (marked by the blue star) at $t_i=3 {\rm \pi}/2$, the path followed would be identical, but the Lagrangian drift described by the corresponding sequence of locations $\boldsymbol {x}_n=\boldsymbol {x}_p(t_i+2 {\rm \pi}n)$, indicated by blue squares, would be radically different, as seen in figure 6.

In trying to characterize the particle drift in a non-ambiguous way, it is therefore better to use instead the average location of the fluid particle during a given cycle $n$, computed according to

(4.2)\begin{equation} \boldsymbol{x}_n=\frac{1}{2{\rm \pi}} \int_{t_i+2 {\rm \pi}(n-1)}^{t_i+2 {\rm \pi}n} \boldsymbol{x}_p \,{\rm d}t. \end{equation}

As can be seen in figure 6, the values of $\boldsymbol {x}_n$ corresponding to $\boldsymbol {x}_i=(-0.55,2.95)$ and $t_i={\rm \pi} /2$, marked by red circles, and those obtained for $\boldsymbol {x}_i=(-2.62,2.87)$ and $t_i=3 {\rm \pi}/2$, marked by blue circles, describe the same path, thereby removing the above-mentioned arbitrariness.

As shown in the fourth column of figure 4, for $\varepsilon \ll 1$ the Lagrangian mean motion features recirculating vortices, whose centre $\boldsymbol {x}_c$ can be determined by computing the location where the Lagrangian streamfunction $\psi _{SS}+\psi _{SD}$ shows a local extremum. Similar recirculating patterns are found for $\varepsilon \sim 1$. In that case, the corresponding vortex centre can be obtained numerically by identifying the location $\boldsymbol {x}_c$ that satisfies $\boldsymbol {x}_n=\boldsymbol {x}_c$, so that the fluid particle describes exactly the same trajectory over subsequent cycles, with zero net drift.

The location of the vortex centre $\boldsymbol {x}_c$ of the Lagrangian mean flow is shown in figure 7 for oscillatory motion with infinitesimally small values of the stroke length $\varepsilon \ll 1$ and also with finite values $\varepsilon =(0.5,1.0,1.5)$. For the Womersley number $M=2$ considered in figure 7, there exists a single vortex, whose centre occupies a location that depends on the inter-cylinder spacing $\ell$. As can be seen, the results are in general agreement with those displayed in figure 2 for the steady-streaming motion, in that, as $\ell$ is reduced, the vortex centre migrates from a location near the ${\rm \pi} /4$ ray towards the vertical axis $x=0$. As expected, the DNS results for increasing stroke lengths $\varepsilon$ progressively depart from the $\varepsilon \ll 1$ predictions, with the vortex centre moving downwards while maintaining approximately the same horizontal location.

Figure 7. The variation with the inter-cylinder distance $\ell$ of the location of the Lagrangian vortex centre $\boldsymbol {x}_c$ for $M=2$ as determined in the limit $\varepsilon \ll 1$ and as determined from the DNS computations with $\varepsilon =(0.5,1.0,1.5)$. The symbols represent the results corresponding to $\ell =(1,1.25,1.5,1.75,2,2.5,3,4,6,8,10,15,\infty )$.

The increasing downward displacement of the vortex centre for increasing $\varepsilon$ shown in figure 7 is accompanied by a progressive distortion of the Lagrangian vortex. This is illustrated in figure 8 for $\ell =2$, with the vortex shape characterized by plotting the time-averaged path of fluid-particle trajectories initiated at points located at increasing vertical distances from the vortex centre, indicated in the figure caption. For each fluid particle, the plot shows a sequence of 80 cycles. Since the Lagrangian velocity is larger for larger $\varepsilon$ (i.e. $\boldsymbol {v}_{L} \propto \varepsilon$ for $\varepsilon \ll 1$, as demonstrated in figure 4), for the same number of cycles, the Lagrangian displacement increases with increasing $\varepsilon$, so that, for instance, the fluid particle closer to the vortex centre describes two laps for $\varepsilon =0.5$ and about 10 laps for $\varepsilon =1.5$.

Figure 8. Lagrangian mean motion for $\ell =2$ and $M=2$, including streamlines $\psi _{SS}+\psi _{SD}=\textrm {const.}$ with $\delta \psi =0.004$ for $\varepsilon \ll 1$ and time-averaged fluid-particle locations $\boldsymbol {x}_n$ for $\varepsilon =(0.5,1.0.1.5)$ computed using (4.2) for the trajectories determined by integrating (4.1) with initial condition $\boldsymbol {x}=\boldsymbol {x}_i$ at $t=0$. In computing the trajectories, the initial locations $\boldsymbol {x}_i$ were selected at fixed vertical distances $\delta y$ above the Lagrangian vortex centre $\boldsymbol {x}_c$, the latter indicated with an asterisk. Five different trajectories corresponding to $\delta y=(0.2,0.4,0.6,0.8,1.0)$ are plotted for $\varepsilon =0.5$ and $\varepsilon =1.0$, whereas, to avoid cluttering, only three trajectories corresponding to $\delta y=(0.2,0.6,1.0)$ are shown in the case $\varepsilon =1.5$.

The numerical results for $\varepsilon =(0.5,1.0,1.5)$ are to be compared with the Lagrangian streamlines computed in the limit $\varepsilon \ll 1$ with use of $\psi _{SS}+\psi _{SD}=\textrm {const}$. As can be seen, the Lagrangian vortex for $\varepsilon =0.5$ is almost indistinguishable from its $\varepsilon \ll 1$ counterpart and, even for the case $\varepsilon =1.0$, the asymptotic predictions provide a fairly good description of the circular drift motion. Departures are more pronounced for $\varepsilon =1.5$ as a result of the increasing nonlinearity. Contrary to the cases $\varepsilon =0.5$ and $\varepsilon =1.0$, for which all time-averaged locations corresponding to a given fluid particle closely lie along a well-defined closed path, for $\varepsilon =1.5$ the locations $\boldsymbol {x}_n$ are scattered within a band surrounding the vortex centre.

The comparisons presented in figures 7 and 8 indicate that the simple velocity description arising for $\varepsilon \ll 1$, in which the Lagrangian mean velocity is given by the sum of distinct steady-streaming and Stokes-drift components, can be used with unexpectedly good accuracy to quantify the fluid-particle drift in situations in which the stroke length is as large as the cylinder radius (i.e. order-unity values of $\varepsilon$) provided that the flow remains symmetric and periodic. In view of previous results pertaining to the single cylinder (Tatsuno & Bearman Reference Tatsuno and Bearman1990), increasing nonlinear effects can be expected to modify significantly the flow pattern depicted in figure 8 as the Reynolds number ${\textit {Re}}=\varepsilon M^2$ increases to sufficiently large values, with the associated Lagrangian motion eventually becoming chaotic, as the flow becomes turbulent; these additional nonlinear effects were not further investigated here.

5. Steady streaming in anharmonically oscillating flows

Most investigations of pulsating flows over cylinders consider outer streams with harmonically varying velocities, resulting in symmetric streaming flows with closed streamlines that are identical in all four quadrants. As shown by Davidson & Riley (Reference Davidson and Riley1972), the classical analysis can be extended to anharmonic flow by expressing the periodic outer velocity as a Fourier series $u_\infty = \sum _{n=1}^{\infty } \textrm {Re} (A_n \mathrm {e}^{\mathrm {i} n t} )$ involving the complex coefficients $A_n$. Correspondingly, the linear problem that arises at leading order in the limit $\varepsilon \ll 1$ can be solved by introducing Fourier-series expansions for the velocity $\boldsymbol {v}_0= \sum _{n=1}^{\infty } \textrm {Re} (A_n \mathrm {e}^{\mathrm {i} n t} \boldsymbol {V}_n )$. For the cylinder array, the complex function $\boldsymbol {V}_n$ corresponding to a given mode $n$ would be obtained by integrating (3.1a,b) subject to (3.2) for a Womersley number $M_n=(a^2 n \omega /\nu )^{1/2}$. In carrying the analysis to the following order, it is important to note that the forcing term $\langle \boldsymbol {v}_0 {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_0\rangle$ that determines the steady streaming through (3.3a,b) and the Stokes drift $\boldsymbol {v}_{SD}=\langle \int \boldsymbol {v}_0 \,\textrm {d}t \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_0 \rangle$ are obtained by time averaging the product of two Fourier series. Since the time average of the product of any two modes of different frequency is identically zero, the resulting functions become

(5.1a,b)\begin{align} \langle {\boldsymbol{v}}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{v}}_0\rangle=\frac{1}{2} \sum_{n=1}^{\infty} |A_n|^2 \,{\rm Re} ({\boldsymbol{V}_{\boldsymbol{n}}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{V}}_n^* ) \quad {\rm and} \quad {\boldsymbol{v}}_{{SD}}=\frac{1}{2} \sum_{n=1}^{\infty} \frac{|A_n|^2}{n} \,{\rm Im} ({\boldsymbol{V}}_n \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{V}}_n^* ), \end{align}

involving the sum of the separate contributions of each mode, with no inter-mode interactions. As a consequence, the steady streaming and Stokes drift generated by an anharmonic flow can be obtained simply as the sum of the corresponding steady-streaming and Stokes-drift velocities of each separate mode. Since each mode gives closed streamlines that are identical in all four quadrants, as those represented in figures 2 and 4, their linear superposition also gives symmetrical recirculatory patterns that are qualitatively similar to those obtained in the harmonic case, thereby maintaining the fore-and-aft symmetry of the flow. It can therefore be concluded that the description of the expected symmetry breaking arising in the presence of anharmonic flow requires consideration of the inter-mode interactions occurring at order $\varepsilon ^2$. These higher-order terms in the asymptotic expansion, which describe the flow asymmetries induced by anharmonic flow, have been computed for circular cylinders and spheres undergoing oscillations with $\varepsilon \ll 1$ and ${\textit {Re}}_s=\varepsilon ^2 M^2 \sim 1$ (Miyagi & Nakahasi Reference Miyagi and Nakahasi1975; Tatsuno Reference Tatsuno1981; Higa & Takahashi Reference Higa and Takahashi1987).

Many oscillatory flow phenomena of physiological interest display an anharmonic time dependence, that being, for example, the case of CSF flow along the spinal canal (Linninger et al. Reference Linninger, Tangen, Hsu and Frim2016). As revealed by magnetic resonance measurements of cardiac-driven motion (Coenen et al. Reference Coenen, Gutierrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019; Sincomb et al. Reference Sincomb, Coenen, Gutiérrez-Montes, Martínez-Bazán, Haughton and Sánchez2022), the flow rate exhibits a non-sinusoidal variation induced by the intracranial pressure, including a short period of fast caudal flow followed by a longer period of slow flow in the cranial direction. Since this pulsating stream interacts with nerve roots and ligaments that are aligned with the flow, a relevant question is whether such interactions can lead to the appearance of a longitudinal streaming motion, which could explain the enhanced transport rate previously observed (Stockman Reference Stockman2006, Reference Stockman2007).

To try to shed light on this matter, effects of anharmonicity were investigated in connection with pulsating flow over the streamwise cylinder array considered here. In view of the previous comments pertaining to flow over a cylinder, it can be expected that for $\varepsilon \ll 1$ the velocity corrections associated with the symmetry breaking are small, of order $\varepsilon ^2$ (Miyagi & Nakahasi Reference Miyagi and Nakahasi1975; Tatsuno Reference Tatsuno1981; Higa & Takahashi Reference Higa and Takahashi1987), so that the appearance of significant asymmetry, possibly leading to a non-zero streamwise flow rate, requires values of the stroke length comparable to the cylinder radius (i.e. $\varepsilon \sim 1$), a problem addressed here with use of DNS simulations. The integrations correspond to a cylinder array with $\ell =2$ for a simple two-term periodic ambient velocity $u_\infty =3 \cos (t)/4+\cos (2t)/4$, whose anharmonic temporal variation is shown in the inset in figure 9($a$).

Figure 9. Time-averaged DNS results corresponding to a cylinder array with $\ell =2$ and $M=(1,4,8,16)$ for the ambient periodic velocity $u_\infty =3 \cos (t)/4+\cos (2t)/4$ represented in the inset in panel ($a$). Panel ($a$) shows the variation with $\varepsilon$ of the rescaled streamwise flow rate $Q/\varepsilon ^2$, while panels ($b$)–($e$) represent streamlines (with spacing $\delta \langle \psi \rangle =0.001$ for $M=1$ and $\delta \langle \psi \rangle =0.006$ for $M=4,8$ and $16$) and vorticity contours for $\varepsilon =1.0$.

The time-averaged Eulerian velocity $\langle \boldsymbol {v} \rangle$ computed for $\varepsilon =1$ was used to determine the streamlines and vorticity shown for four different values of $M$ in figure 9($b$–$e$). The plots include the first two quadrants, as needed to illustrate the lack of fore-and-aft symmetry, which is less pronounced for $M=1$. For larger values of $M$, the time-averaged flow comprises two highly distorted vortices in the vicinity of the cylinder, surrounded by a region of nearly horizontal flow with velocities that decay slowly with distance. The comparison of the streaming results for $M=1$ and $M=16$ with those shown earlier in the second column of figures 4(a) and 4(c) for the harmonic case clearly indicate that the effects of anharmonicity are much more important for larger values of $M$, for which the outer vortex is replaced by a streamwise current, which is absent in the case $M=1$.

The streamline pattern shown in the plots for $M \ne 1$ is consistent with the existence of a non-zero streamwise flow rate $Q=\int _0^\infty \langle u \rangle (\ell,y) \,\textrm {d} y$ (or $Q=\int _1^\infty \langle u \rangle (0,y) \,\textrm {d} y$). The variation of $Q$ with $\varepsilon$, determined in the DNS integration from the value of the time-averaged streamfunction $\langle \psi \rangle$ in the far field, is shown in figure 9 for different values of $M$. The plot reveals that the proportionality $Q \propto \varepsilon ^2$, to be expected for $\varepsilon \ll 1$, continues to apply over the whole range of $\varepsilon$ considered in the DNS, for which the ratio $Q/\varepsilon ^2$ remains approximately constant. The negative value of $Q/\varepsilon ^2$, negligibly small for $M=1$, increases in magnitude for increasing $M$, reaching $Q/\varepsilon ^2 \simeq -0.58$ for $M = 16$.

6. Concluding remarks

The interaction of an oscillating stream with a streamwise linear array of cylinders gives rise to a stationary motion that has been quantified here for configurations with Womersley numbers $M$ of order unity and dimensionless stroke lengths $\varepsilon$ that are either $\varepsilon \ll 1$ or $\varepsilon \sim~1$, thereby yielding moderately small values of the Reynolds number ${\textit {Re}}=\varepsilon M^2=U_\infty a/\nu$, for which the flow remains two-dimensional, time-periodic and symmetric with respect to the centreline. For infinitesimally small values of $\varepsilon$, the Lagrangian mean motion is obtained as the sum of the steady-streaming and Stokes-drift components, which have been computed for different values of $M$ and of the inter-cylinder spacing $\ell$. The description has been validated by comparisons with results of DNS involving finite values of $\varepsilon$. The comparisons revealed, perhaps unexpectedly, that the simplified description for $\varepsilon \ll 1$ continues to give reasonably accurate predictions for the time-averaged Eulerian velocity and for the Lagrangian mean motion as the stroke length increases to values of order unity (see, in particular, the results shown for $\varepsilon =0.5$ in figures 3 and 8).

While most of the analysis focuses on oscillating streams with harmonic velocity, consideration is also given to the effects of anharmonicity, an analysis motivated by oscillatory flow phenomena of physiological interest. An important conclusion of our study is that the interaction of an anharmonic stream with a streamwise obstacle array can have a profound effect on the convective transport rate, especially in configurations with $\varepsilon \sim 1$ and large values of $M$, for which the presence of the array can be expected to induce a streamwise flow rate of order $U_\infty a$, corresponding to order-unity values of the dimensionless flow rate $Q$ shown in figure 9.

Further investigation is warranted to assess the importance of these effects in connection with the motion of CSF in the spinal canal, as needed to improve predictive capabilities of current flow and transport models (Sánchez et al. Reference Sánchez, Martínez-Bazán, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019; Sincomb et al. Reference Sincomb, Coenen, Gutiérrez-Montes, Martínez-Bazán, Haughton and Sánchez2022). To enable quantitative predictions, these future investigations should consider more realistic geometrical configurations, including annular models of the spinal canal with obstacles arranged longitudinally to represent the ventral and dorsal nerve roots (Stockman Reference Stockman2006, Reference Stockman2007). The results in § 5 suggest that the contribution of the induced Lagrangian motion to the streamwise transport rate is likely to be more prominent in the cervical region, where we find larger values of $\varepsilon$, while associated contributions in the lumbar region will be necessarily more limited.