Lubrication analysis of peristaltic motion in non-axisymmetric annular tubes

Abstract This paper addresses peristaltic flow induced in a non-axisymmetric annular tube by a periodic small-amplitude wave of arbitrary shape propagating axially along its inner surface, assumed to be a circular cylinder. The study is motivated by recent in vivo experimental observations pertaining to the flow of cerebrospinal fluid along the perivascular spaces of cerebral arteries. The analysis employs the lubrication approximation, describing low-Reynolds-number peristaltic flow in the long-wavelength approximation. Closed-form analytic expressions are derived for the average pumping rate in infinitely long tubes and also in tubes of finite length. Consideration is also given to the transverse motion arising in non-axisymmetric tubes. For small-amplitude waves, the solution is reduced to the integration of a parameter-free Stokes-flow problem, which is solved for relevant cross-sectional shapes, with closed-form analytical results derived for thin canals.


Introduction
The transport of fluid associated with the progression of contraction/expansion waves along a tube plays a fundamental role in many physiological phenomena and technological applications (Jaffrin & Shapiro 1971). In many situations, the Reynolds number is small and the wavelength is large compared with the tube transverse dimension, so that the motion can be described in the lubrication approximation. The seminal analysis of Shapiro, Jaffrin & Weinberg (1969), corresponding to peristaltic motion driven by a train of sinusoidal waves travelling along an infinitely long tube of circular section, has been † Email address for correspondence: als@ucsd.edu extended to account for effects of different cross-sectional shapes (Rath 1982) and finite tube lengths (Li & Brasseur 1993), for example.
Peristaltic pumping driven by arterial pulsations due to the heartbeat has been reasoned to play a role in the transport of cerebrospinal fluid (CSF) along the perivascular spaces that surround cerebral arteries (Rennels et al. 1985;Stoodley et al. 1997;Hadaczek et al. 2006). The recent in vivo particle image velocimetry (PIV) measurements in mice of Mestre et al. (2018) have demonstrated that the resulting pulsatile motion of the CSF exhibits a net flow rate in the direction of the arterial blood flow -see also Raghunandan et al. (2021) for more recent experiments. The experiments also revealed that the periarterial space is a mostly unobstructed non-axisymmetric annular canal whose width is comparable to the arterial radius (Min Rivas et al. 2020). As reasoned by Thomas (2019), because of the strong asymmetry of the cross-section, the associated flow characteristics and resulting pumping rate can be expected to differ significantly from those predicted for an axisymmetric annulus (Bilston et al. 2003;Wang & Olbricht 2011). The effects of asymmetry have been recently addressed using an inner circle (representing the artery) and an outer eccentric ellipse as an approximate model for the cross-sectional shape of the periarterial space. Attention has been given to the associated hydraulic resistance, which was evaluated by Tithof et al. (2019), and to the periodic flow induced by a small-amplitude sinusoidal wave, which was described by Carr et al. (2021) by integrating numerically the full Navier-Stokes equations over five cycles in a tube of length equal to two wavelengths. The latter analysis numerically explored the dependence of the pumping rate on the cross-sectional shape and pointed out the inherent three-dimensionality of the flow, arising in the absence of axial symmetry. These features of the solution are further investigated here by means of analytical methods that exploit simplifications arising for low-Reynolds-number slender flows. Despite the numerous simplifying approximations introduced in the description, the results have potential applicability in connection with the flow in the cerebral perivascular system, as discussed at the end of the paper.
As in the previous numerical simulations of Carr et al. (2021), the present work considers peristaltic flow in a non-axisymmetric annular canal. Cylindrical coordinates (x , r , θ) and corresponding velocity components (u , v , w ) are to be used in the analysis. As indicated in figure 1, the canal is bounded externally by a fixed cylinder of radius r e (θ ), while the inner surface, representing the arterial wall supporting a travelling peristaltic wave, has a circular cross-section of radius r a , comparable in magnitude to the outer radius r e , given by where t is the time, and κ and ω are the wavenumber and angular frequency, with 2π/κ representing the wavelength, assumed to be much longer than the characteristic canal transverse dimension (i.e. κr o 1). The wave is described by the general 2π-periodic function T (ξ ) of the characteristic coordinate ξ = κx − ωt . This wavefunction satisfies 2π 0 T dξ = 0, so that r o is the time-averaged radius. In periarterial flow, the specific shape of T (ξ ), schematically represented in figure 1, is related to the cardiac pulsation, as revealed by in vivo temporal measurements of arterial diameters in the mouse brain (Mestre et al. 2018), which also showed the relative wave amplitude ε to be small, on the order of ε ∼ 0.02.

The infinitely long tube
We begin by considering the case of an infinitely long tube undergoing peristaltic motion, for which the 2π-periodic variations in the rescaled coordinates t = ωt and x = κx can be described in terms of the single characteristic coordinate ξ = κx − ωt . As discussed by Shapiro et al. (1969), most of the results, including the associated volume flow rate Q (ξ ), are in general applicable to tubes whose length L is a multiple of the wavelength (i.e. L = 2πj/κ, with j = 1, 2, 3, . . .).
The fluid in the annular canal has density ρ and viscosity μ. To focus more directly on the conditions of interest for periarterial CSF flow, it is assumed that the characteristic viscous time r 2 o /(μ/ρ) satisfies r 2 o /(μ/ρ) ω −1 , so that the motion is quasi-steady and with negligibly small inertia. The resulting Stokes equations further simplify as a result of the slenderness condition r o κ −1 , in that, with errors of order (κr o ) 2 , longitudinal derivatives can be neglected when computing the viscous stresses. At the same level of approximation, the transverse variations of the pressure, which are a factor (κr o ) 2 smaller than the corresponding longitudinal variations, can be discarded when computing the axial flow.
The dimensionless formulation employs the coordinates ξ = κx − ωt and r = r /r o , along with the order-unity velocity components u = u /(εω/κ), v = v /(εωr o ) and w = w /(εωr o ), and corresponding flow rate Q = Q /(εr 2 o ω/κ). The accompanying spatial pressure variations p are expressed in the form where the functionp(ξ, r, θ) is introduced to measure the small pressure variations, of order (κr o ) 2 1, found across the section of the canal, which enter in the description of the transverse motion. In the lubrication approximation adopted here, the problem reduces to the integration of

The axial velocity
Following standard practice (White 2006), the axial velocity can be expressed in the compact form u = −(dp/dξ)U, where the function U(r, θ; a) satisfies the Poisson problem 1 r The accompanying volume flow rate becomes in terms of the dimensionless hydraulic resistance Simple solutions to (3.1) are available when the tube is a concentric circular annulus of outer radius r e = b (White 2006) and also for thin annular tubes of general shape satisfying r e − a 1, when the solution can be approximated by the expressions Numerical integration is in general needed to solve (3.1) for non-axisymmetric tubes with r e − a ∼ 1, with results relevant to the periarterial cross-section computed by Tithof et al. (2019). As in Carr et al. (2021), illustrative results will be given below for two different outer boundaries r e (θ ), depicted at the top of figure 2, namely, a coaxial elliptical cylinder with semi-axis lengths br o and cr o , a configuration investigated by Shivakumar (1973), and a circular cylinder of radius br o whose centre is displaced a distance cr o from that of the inner circular cylinder, a configuration first studied by Piercy, Hooper & Winny (1933). In the latter case, the approximate expression yields very accurate results provided that b is not too large (White 2006).  The function U and associated hydraulic resistance R depend on the characteristic coordinate ξ through the inner radius a = 1 + εT (ξ ). The dependence is weak for cases with ε 1, to be investigated below, for which expanding R in a Taylor series leads to where both evaluated at a = 1, are, respectively, the hydraulic resistance of the tube with constant inner radius r a = r o and the relative variation of the hydraulic resistance in response to changes of the inner radius. Selected values of R o and Δ, determined numerically with use of (3.1) and (3.3), are shown in figure 2 for the eccentric and elliptical annuli. For the former, figure 2 also shows predictions obtained using the simple approximate formulae stemming from (3.5).

The pumping rate
The determination of the axial pressure gradient dp/dξ , which enters as a factor in u = −(dp/dξ)U, begins by integrating (2.5) across the tube section with the boundary conditions (2.6) to give upon substitution of (3.2) and a = 1 + εT (ξ ). Integrating the above equation once gives involving the unknown integration constant C λ . Since 2π 0 T dξ = 0, taking the average of the above equation over a cycle provides for the time-averaged pumping rate at any location. To facilitate the notation, averages of 2π-periodic functions of a single independent variable are denoted throughout the paper by · = 2π 0 (·) dy/(2π), where y is a generic dummy integration variable, i.e. y = ξ in (4.3).
The value of the constant C λ in (4.3) can be obtained by multiplying (4.2) by R and integrating the resulting equation between ξ = 0 and ξ = 2π to give For generality, following the work of Shapiro et al. (1969), in writing (4.4) we have assumed that the periodic streamwise pressure gradient dp/dξ may have a steady component, associated with a pressure rise δp λ per wavelength. The above results, which are valid for ε ∼ 1, can be simplified in the small-amplitude limit ε 1. Using the expression (3.6) in (4.4) and neglecting terms of order ε 2 and smaller provides C λ = 2πε(Δ + 1/2) T 2 − δp λ /(2πR o ), which can be substituted into (4.3) to finally give for the mean pumping rate. The simple theoretical formula (4.5) is valid for infinitely long tubes and also for tubes with length L = 2πj/κ, with j = 1, 2, 3, . . .. There is interest in comparing the resulting prediction with the recent direct numerical simulations (DNS) of Carr et al. (2021), corresponding to a tube of length L = 4π/κ with δp λ = 0. The DNS considered a sinusoidal wave T = sin ξ with relative amplitude ε = 0.02 for two different geometries, namely, an eccentric annulus with b = √ 2.4 and eccentricity 0 ≤ c ≤ 0.349, and an elliptical annulus with bc = 2.4 and 1 ≤ b/c ≤ 1.667. The mean flow rate was evaluated for 10 different configurations, with results represented in Carr et al. (2021) in terms of the dimensionless variable Q * = Q /(1.4πr 2 o ω/κ), related to our steady pumping rate by Q = 1.4πQ * /ε. The 10 different points are shown in figure 2 along with the theoretical prediction Q = πεΔ, obtained from (4.5) for T 2 = sin 2 ξ = 1/2 and δp λ = 0 using the values of Δ corresponding to b = √ 2.4 (eccentric circles) and bc = 2.4 (elliptical tube), respectively. As can be seen, the analytical predictions agree reasonably well with the numerical results, with relative departures remaining below 15 %. The observed differences can be attributable to the fact that the asymptotic analysis employs the long-wave limit κr o 1, while the wavelength in the computations was only moderately long, yielding values of κr o 0.5. The relative differences between the computations and the theoretical predictions remain of the order of (κr o ) 2 , as is consistent with the accuracy of the asymptotic development leading to (4.5).
The above problem can be solved analytically for thin annular tubes of general shape, defined by the width distribution h(θ ) = r e − 1 1. The simplified solution for the axial motion, given earlier in (3.4a,b), can be written in the form U o = y(h − y)/2 and R −1 o = ( 2π 0 h 3 dθ)/12, involving the radial distance to the inner cylinder y = r − 1. These simplified expressions can be used in (5.4) to evaluate the function q, so that (5.3) takes the reduced form ∂V ∂y Integrating with the boundary condition V = −1 at y = 0 yields when evaluated at y = h, where V = 0. The analysis continues by noting that, for h 1, the radial velocities and the small radial variation of the pressure can be neglected in (5.2) when computing the azimuthal velocity, Substituting the associated volume flux h 0 W dy = −(dP/dθ)h 3 /12 into (5.8) yields a second-order equation for the pressureP. A first integral gives − dP dθ which can be used to evaluate the azimuthal pressure gradient dP/dθ, thereby completing the determination of the azimuthal velocity (5.9) and of the radial velocity, the latter being obtained from (5.7) with the use of (5.9). The constant of integration in (5.10) is written in terms of the angular location θ o at which the azimuthal volume flux h 0 W dy vanishes. Its value can be obtained by integrating (5.10) a second time after dividing by h 3 and imposing the condition thatP must be periodic (i.e.P(0) =P(2π)). Note that for cross-sections that exhibit symmetry, there is more than one possible value of θ o , i.e. θ o = (0, π) for the eccentric annulus and θ o = (0, π/2, π, 3π/2) for the elliptical annulus.
In general, numerical integration is needed to solve (5.1)-(5.5) for cross-sections with r e ∼ 1, with illustrative results given in figure 3 for the two geometries considered here. In each case, the parameters b and c defining the shape are selected to be those employed in the previous simulations (Carr et al. 2021), thereby facilitating comparison. Figure 3 shows the projection of the streamlines onto the tube section, i.e. the curves that are tangent to the transverse velocity (v, w) =Ṫ (V, W). Since both transverse velocity components are proportional to the functionṪ (κx − ωt ), the resulting streamline pattern is independent of time and identical at all tube sections. The projected streamlines are perpendicular to the inner cylinder and tangent to the outer surface (except at the critical points θ = θ o ), a flow structure also revealed in the previous computations of Carr et al. (2021).

Tubes of finite length
Peristaltic motion in tubes of arbitrary length was first investigated by Li & Brasseur (1993). Unlike the infinitely long tube, for which the single characteristic variable ξ = κx − ωt suffices to describe the temporal and streamwise dependence of the resulting periodic flow, for a tube of arbitrary finite length one needs to consider in general two different independent variables, because the solution is 2π-periodic in the rescaled time t = ωt but non-periodic in the spatial variable x = κx . It is seen that the use of ξ = κx − ωt and t = ωt as independent variables facilitates the description of the flow rate Q(ξ, t) and pressure p(ξ, t).
The initial steps of the analysis leading to (4.1) parallel those presented above, with the pressure determined from where = κL is the dimensionless tube length and δp L (t) = δp L /[μεω/(κr o ) 2 ] is the dimensionless pressure difference between the inlet x = 0 (ξ = −t) and the outlet x = L (ξ = − t). Integrating (6.1) once gives involving the periodic function C L (t). Multiplying by R and integrating a second time yields after imposing the boundary conditions at both ends. As expected, since R, T anḋ T = dT /dξ are 2π-periodic functions of ξ , the expression (6.3) naturally reduces to (4.4) with δp λ = δp L /j when L is a multiple of the wavelength, so that = κL = 2πj, with j = 1, 2, 3, . . . . Substituting (6.3) into (6.2) provides the flow rate Q(x, t), whose dependence on the streamwise coordinate x enters through the wavefunction T (ξ ) = T (x − t). Taking the time average of (6.2) at any given section x reveals that 1 2π 2π 0 Q dt = C L − επ T 2 (6.4) in terms of the mean value C L . The solution simplifies for small peristaltic amplitudes ε 1. Using (3.6) in (6.3), substituting the result into (6.2), and taking the time average gives, with small errors of order ε 2 , 1 2π for the time-averaged flow rate. As expected, the above equation reduces to (4.5) for = 2πj, when δp L = jδp λ and −t −t T dξ = 0. For a sinusoidal wave T = sin ξ in a tube with equal pressure at both ends, (6.5) reduces to 1 2π to be compared with the result 2π 0 Q dt/(2π) = πεΔ obtained earlier for an infinitely long tube. The comparison indicates that the effect of non-integral length / = 2πj is to reduce the level of pumping, in agreement with previous findings (Li & Brasseur 1993).

Conclusions
A complete analytical description has been given for the lubrication flow induced in non-axisymmetric tubes by a periodic wave propagating axially along its inner surface, assumed to be a circular cylinder. The explicit expressions derived for the pumping rate, involving the dimensionless hydraulic resistance R o and its relative variation in response to changes of the inner radius Δ, can be used to assess the role of peristaltic pumping in applications involving transport along asymmetric annuli. In particular, the results given in § 6 can in principle be used in connection with the problem of CSF flow in the perivascular system, with the overall network of perivascular spaces surrounding arteries assumed to be composed of multiple peristaltic elements of length L -possibly different for different elements -connected at bifurcating nodes. Since the wavelength 2π/κ of the arterial pulsations is very large compared with the length L of each peristaltic element, the expression given above in (6.5) for the relation between δp L and Q could be simplified by using the limit = κL 1. Combining the δp L -Q relation corresponding to the different peristaltic elements with the continuity condition at each bifurcating node provides a set of equations that determine the distribution of flow and pressure throughout the periarterial system. Simplified computational approaches like this can be instrumental in understanding the overall features of CSF flow, complementing on-going efforts .
Funding. This work was supported by the National Science Foundation through grant number 1853954. W.C. acknowledges partial support from the 'Convenio Plurianual Comunidad de Madrid -Universidad Carlos III de Madrid' through grant CSFLOW-CM-UC3M.
Declaration of interests. The authors report no conflict of interest.