Stationary flow driven by non-sinusoidal time-periodic pressure gradients in wavy-walled channels

The classical problem of secondary flow driven by a sinusoidally varying pressure gradient is extended here to address periodic pressure gradients of complex waveform, which are present in many oscillatory physiological flows. A slender two-dimensional wavy-walled channel is selected as a canonical model problem. Following standard steady-streaming analyses, valid for small values of the ratio ε of the stroke length of the pulsatile motion to the channel wavelength, the spatially periodic flow is described in terms of power-law expansions of ε, with the Womersley number assumed to be of order unity. The solution found at leading order involves a time-periodic velocity with a zero time-averaged value at any given point. As in the case of a sinusoidal pressure gradient, effects of inertia enter at the following order to induce a steady flow in the form of recirculating vortices with zero net flow rate. An improved two-term asymptotic description of this secondary flow is sought by carrying the analysis to the following order. It is found that, when the pressure gradient has a waveform with multiple harmonics, the resulting velocity corrections display a nonzero flow rate, not present in the single-frequency case, which enables stationary convective transport along the channel. Direct numerical simulations for values of ε of order unity are used to investigate effects of inertia and delineate the range of validity of the asymptotic limit ε≪1. The comparisons of the time-averaged velocity obtained numerically with the two-term asymptotic description reveals that the latter remains remarkably accurate for values of ε exceeding 0.5. As an illustrative example, the results of the model problem are used to investigate the cerebrospinal-fluid flow driven along the spinal canal by the cardiac and respiratory cycles, characterized by markedly non-sinusoidal waveforms.


Introduction
The interaction of an oscillating stream with no-slip boundaries gives rise to a time-averaged steady-streaming motion [1].The vast majority of studies dealing with such flows employ a sinusoidal signal to generate the primary oscillating flow.Thus, most investigations of flow over an oscillating solid body immersed in a stagnant fluid [2,3] or of flow over a fixed solid body placed in an oscillatory stream [4] assume sinusoidal oscillations, while for wall-bounded flows either the driving pressure gradient [5], the prescribed stroke volume [6], or the wall velocity [7] are assumed to vary sinusoidally in time.For geometrical configurations with streamwise symmetry, including circular cylinders [4] wavy-walled channels [8] and pipes [9], the resulting steady streaming develops in the form of recirculating vortices exhibiting symmetric closed streamlines.An important consequence is that a sinusoidal pressure gradient acting on a wall-bounded flow having streamwise symmetry, such as a wavy-walled channel [10,11], is unable to generate a steady streamwise flow rate.
Most oscillatory physiological flows driven by the cardiac and respiratory cycles display a periodic time dependence with multiple harmonics.For example, the intracranial pressure fluctuation driving the motion of cerebrospinal fluid (CSF) in the spinal canal has a waveform that includes three different peaks associated with the cardiac cycle, whose amplitudes decrease in a stepwise manner in a healthy individual [12].The flow induced in the canal has been shown to have a stead-ysteaming component [13] that contributes fundamentally to the Lagrangian transport rate [14].The presence of anatomical features obstructing the flow in the canal, including nerve roots and ligaments, has been shown to increase the magnitude of the resulting steady-streaming velocity [15].Previous analytical investigations addressing this aspect of the CSF motion have employed as a canonical flow configuration a wavy-walled channel, with the associated wavelength representing the intervertebral distance [10,11,16].These analyses, assuming a sinusoidally varying pressure gradient, showed that steady streaming occurs in the form of closed recirculating vortices confined to the channel cells, with no streamwise connectivity.Although it is unclear how a pressure gradient of complex waveform may change this symmetric recirculating flow pattern, it is worth recalling that the early analyses of secondary flows induced by cylinders and spheres undergoing non-sinusoidal oscillations [17][18][19] revealed the appearance of a small longitudinal flow rate aligned with the oscillation axis, not present in the sinusoidal case.Clearly, in wall-bounded flows the possible emergence of a net flow rate when the fluid is subject to a pressure gradient with complex waveform could have a large impact on the resulting longitudinal transport rate.This aspect of the problem is to be investigated below using the wavy-walled channel as illustrative canonical configuration, thereby extending the previous studies pertaining to a pure sinusoidal waveform [10,11,16].
The analysis considers the flow induced by a general time-periodic pressure difference with zero mean, for values of the oscillation period on the order of the viscous time across the canal.An asymptotic analysis is performed for values of the stroke length much smaller than the channel wavelength, a familiar limit [1] in which the problem is linear at leading order, resulting in an oscillatory velocity with zero mean that can be expressed in Fourier form [20].The first-order velocity corrections contain a stationary component in the form of closed symmetric recirculating vortices with no net flow rate.Symmetry breaking leading to a stationary streamwise flow rate arises at the following order in the expansion for small stroke lengths as a result of the interactions between different modes in the Fourier expansions.This finding is consistent with the previous analyses of the flow induced by oscillating circular cylinders and spheres [17][18][19].
The flow velocity predicted by the asymptotic analysis, including a small net flow rate varying quadratically with the stroke length, is to be compared below with results of direct numerical simulations.The comparisons reveal that the accuracy of the asymptotic results extends beyond the expected range of validity, with the predicted mean flow rate showing good agreement with the numerical results for values of the stroke length exceeding half of the channel wavelength.The influence of the waveform on the resulting net flow rate is illustrated by considering the pressure differences that characterize the cardiac and respiratory motion along the spinal canal, with the induced streamwise flow rate found to be more prominent in the former case.

Problem formulation
Consider a wavy channel of average gap size ℎ 0 filled with a Newtonian fluid of density ρ and kinematic viscosity ν.The channel is bounded by a flat surface and a wavy wall of wave length λ ≫ ℎ 0 , so that the resulting channel width ℎ is described by the periodic dimensionless function ℎ/ℎ 0 = H x′/λ , where x′ is the longitudinal distance measured along the channel.Although the analysis is carried out for wavy walls of general shape, the figures plotted below as well as the accompanying direct numerical simulations correspond to H = 1 + βcos 2πx′/λ , where β < 1 is the relative amplitude of the wall undulation (see Fig. 1).In previous steady-streaming analyses of wall-bounded flows [11], the wavy-walled channel has been reasoned to be a model for the subarachnoid space surrounding the spinal cord, for which ℎ 0 ≃ 1 − 4 mm, with the wave length λ ≃ 2 − 4 cm representing the distance between vertebrae, which cause canal obstructions that can be estimated to correspond to undulations in the range 0.2 ≲ β ≲ 0.4 [21,22].
We investigate the spatially periodic flow driven by a time-periodic pressure gradient with zero mean, resulting in a pressure difference per wavelength ΔpΠ ωt′ , where Δp is the characteristic amplitude, ω is the angular frequency and t′ represents the time.The 2π periodic dimensionless function Π ωt′ , describing the waveform of the pressure fluctuation, is expressible in the general Fourier-expansion form Π ωt′ = ∑ n = 1 ∞ ℜ A n exp inωt′ , where the order-unity coefficients A n are complex numbers and ℜ denotes the real part of a complex function.Since the flow is spatially periodic, it suffices to investigate the solution in the channel stretch 0 ≤ x′ ≤ λ, with the spatial pressure variation p′ satisfying p′ = 0 at x′ = 0 and p′ = ΔpΠ ωt′ at x′ = λ.The description employs cartesian coordinates x′, y′ , with y′ measured from the flat surface, and corresponding velocity components u′, v′ .
The resulting slender flow is characterized by longitudinal velocities of order u c = Δp/ ρωλ , as follows from a balance between the local acceleration ∂u′/ ∂t′ u c ω and the pressure gradient ρ −1 ∂p′/ ∂x′ Δp/ λρ , and much smaller transverse velocities of order v c = ℎ 0 /λ u c ≪ u c , as follows from the continuity balance ∂u′/ ∂x′ ∂v′/ ∂y′.Introducing the dimensionless variables reduces the conservation equations to (3) The velocity is periodic in x, i.e. u 0, y, t = u 1, y, t and v 0, y, t = v 1, y, t , and satisfies the no-slip boundary condition whereas the pressure p is identically zero at x = 0 and takes the value p = Π t at x = 1, with The sinusoidal case can be investigated by considering A 1 = 1 with A n = 0 for n ≥ 2. A pressure signal of complex waveform can be obtained by combining different modes with nonzero coefficients A n .For instance, the pressure function Π t corresponding to A 1 = 3/4 and A 2 = 1/4, used for some of the sample computations presented below, is shown in Fig.

2(a).
As can be seen, the solution depends on two flow parameters, namely, the Womersley number representing the ratio of the viscous time across the channel ℎ 0 2 /v to the characteristic oscillation time ω −1 , and the Strouhal number ε −1 , with representing the ratio of the stroke lengths u c /ω to the channel wavelength λ.The solution exhibits an additional dependence on the slenderness ratio ℎ 0 /λ ≪ 1.
The direct numerical simulations to be presented later employ the complete Eqs. ( 2)- (4).By way of contrast, the analytical development, using asymptotic expansions in powers of ε ≪ 1, will exploit the slenderness of the flow in simplifying the equations by assuming ℎ 0 /λ ≪ ε.Under those conditions, the slender-flow approximation, used in other analyses of steady streaming in channels [5][6][7][8]10,11,16], applies to the description of terms up to order ε 2 , with (4) reducing to ∂p/ ∂y = 0 and (3) yielding to be used together with (2) as starting point in the asymptotic analysis.
A key feature of the solution concerns the resulting mean flow rate Q = ∫ 0 H u dy, where represents the time-averaging operator.While Q was previously found to be identically zero for flow driven by a sinusoidal pressure gradient in infinite wavy-walled channels [10,11,16] and also in general channels with H 0 = H 1 [23], in the nonsinusoidal case treated here the flow rate will be shown to exhibit a small steady component of order ε 2 .

Asymptotic analysis for small stroke lengths
In the asymptotic limit ε ≪ 1 the slender-flow problem can be solved by substituting the asymptotic expansions u = u 0 x, y, t, τ + εu 1 x, y, t, τ + ⋯, v = v 0 x, y, t, τ + εv 1 x, y, t, τ + ⋯, and p = p 0 x, t, τ + εp 1 x, t, τ + ⋯ into ( 2) and ( 9) and solving sequentially the problems that arise at different orders in powers of ε.The solution, determined in [16] for the canonical case Π t = cos t up to terms of order ε, is extended here to order ε 2 for driving pressure differences of general waveform (6).The corresponding expansion for the mean flow rate with Q i = ∫ 0 H u i dy for i = 0,1, 2, …, will be shown to give Q 0 = 0 and Q 1 = 0, so that the analysis needs to be carried to order ε 2 to determine the first nonzero term Q 2 ≠ 0.

Leading-order solution
The linear problem that arises at zeroth order, involving the integration of ∂u 0 ∂x + ∂v 0 ∂y = 0 and ∂u with boundary conditions u 0 = v 0 = 0 at y = 0, H x , p 0 = 0 at x = 0 and p 0 = Π t at x = 1, can be solved exactly to give where with the auxiliary functions and involving the complex Womersley function evaluated with the frequency nω at the local height H x .A dummy integration variable y is introduced for convenience in ( 14).As is clear from ( 13), the leading-order solution has zero time-averaged velocities u 0 = v 0 = 0. Consequently, the leading-order volumetric flow rate also has a zero mean value Q 0 = 0. Sample flow rates Q 0 t computed from ( 19) for the pressure wave form of Fig.

First-order corrections
The solution at the following order ε involves the integration of with boundary conditions u 1 = v 1 = 0 at y = 0, H and p 1 = 0 at x = 0,1.The convective acceleration corresponding to the leading-order solution, driving the velocity corrections at this order, can be expressed in the form where for n ≥ 0 and F n = F −n * for n < 0, with the asterisk * denoting complex conjugates.Note that the second line of the above equation gives no contribution to the evaluation of F 0 and F ±1 .
With the exception of which is real, all F n x, y are in general complex functions.Also of interest is that in the sinusoidal case and all other F n x, y being identically zero.

Steady velocity corrections-The velocity corrections
u can be expressed as the sum of a steady part u 1 s , v 1 s , corresponding to a time-independent secondary flow, and an unsteady part u 1 u , v 1 u , with the latter satisfying The reduced equations for the steady components follow from taking the time average of ( 20) and ( 21) to give which can be integrated with boundary conditions u 1 s = v 1 s = 0 at y = 0, H and p 1 s = 0 at x = 0,1 to yield and where the pressure gradient can be shown to reduce to when account is taken of the condition H 0 = H 1 .The steady flow rate at this order is identically zero, i.e.Q 1 = ∫ 0 H u 1 s dy = 0, a result in agreement with previous findings regarding steady streaming in channels [10,11,16] and tubes [5].
Before we proceed further with the asymptotic description, it is instructive to analyze the characteristics of the streaming flow emerging at this order, described by the above Eqs.( 28)-( 30).Observation of (24) reveals that the steady-streaming velocities ( 28) and ( 29) result from nonlinear interactions of each Fourier mode with itself, a result encountered by Davidson and Riley [20] in their analysis of flow over oscillating cylinders.In the absence of inter-mode interactions, the associated stationary motion is just the linear superposition of the steady-streaming velocities associated with sinusoidal pressure gradients whose frequencies differ by an integer factor.Since each mode produces stationary vortices with no streamwise flow rate, their linear combination is also a periodic array of separate stationary vortices recirculating in each channel cell with no net streamwise flow rate.
To illustrate this feature of the solution, the streamlines and vorticity corresponding to the non-sinusoidal pressure of Fig. 2(a) are plotted in Fig. 3 for β = 0.4 and for two values of the Womersley number, α = 4 and α = 16.As can be seen, the flow structure, similar to that obtained in previous channel-flow analyses [11,16], exhibits recirculating vortices separated by the vertical streamlines x = 0,0.5,1.For α = 4 the flow in each half cell exhibits two counter-rotating vortices, whereas for α = 16 the flow develops an additional, much weaker vortex, located near the section with largest width.The streamlines are plotted using fixed increments δψ 1 of the streamfunction ψ 1 , defined in the usual way (i.e. by integration of ∂ψ 1 / ∂y = u 1 s and ∂ψ 1 / ∂x = − v 1 s with ψ 1 = 0 along the wall), so that the interline spacing provides a measure of the local velocity.To further quantify the motion, color contours are used to represent the associated vorticity Ω = ℎ 0 /λ 2 ∂v/ ∂x − ∂u/ ∂y, which reduces to Ω = − ∂u/ ∂y in the slender flow approximation employed here.
As previously mentioned, the streaming flow at this order has a zero flow rate s dy = 0, consistent with the recirculating patterns depicted in Fig. 3.A nonzero flow rate will be seen to emerge at the following order in the description, to be presented below.Similar higher-order steady-streaming analyses have been performed for flows over oscillatory circular cylinders and spheres [17][18][19], for which the complex waveform breaks the fore-and-aft symmetry of the flow.

Unsteady velocity corrections-
The unsteady component of the first-order velocity corrections can be expressed in the general form 31) where the complex functions U ˆn x, y and V ˆn x, y satisfy with boundary conditions U ˆn = V ˆn = 0 at y = 0, H and P ˆn = 0 at x = 0,1.The problem can be integrated to yield where the pressure gradient can be seen to reduce to when account is taken of the condition H 0 = H 1 .The functions G n and ∫ 0 y G n dy appearing in (33) are defined in ( 16) and ( 17), with n used to evaluate these expressions when n < 0.
The additional auxiliary functions T n and ∫ 0 y T n dy can be evaluated from which completes the description of the first-order corrections.It is worth noting that, in the sinusoidal case, the only nonzero terms in (31) are n = ± 2, so that the unsteady component of the velocity correction can be reduced with use of (26) to where ℑ denotes the imaginary part of a complex function.

Higher-order corrections
Just like the first-order corrections, the velocity corrections u 2 , v 2 emerging at order ε 2 include a time-independent component u 2 s , v 2 s = u 2 , v 2 , which can be obtained by time averaging the corresponding conservation equations to give to be integrated with boundary conditions u 2 s = v 2 s = 0 at y = 0, H and p 2 s = 0 at x = 0,1.The forcing term arises from time-averaging the convective acceleration at this order.In the sinusoidal case one finds ℱ = 0, because u 0 , v 0 ∝ e it and u 1 u , v 1 u ∝ e 2it , so that the solution to (38) reduces to u 2 s = v 2 s = 0. On the other hand, for a non-sinusoidal periodic waveform, the forcing term can be expressed in the form which can be evaluated by substituting the velocity components shown in ( 14) and (33).
Except for the forcing function in the momentum equation, the problem defined above for s is identical to that determining u 1 S , v 1 S , so that the solution can be trivially obtained by replacing F 0 with ℱ in ( 28) and ( 29) to give where the pressure gradient is given in this case by ( The value of u 2 s given in ( 41) can be used to compute s dy, which can be substituted into (11) to give describing the mean flow rate with small relative errors of order ε.As can be seen in ( 44), the value of Q is directly proportional to ε 2 and α 2 , with the Womersley number entering also through the forcing term ℱ.

Second-order steady streaming
The integral expressions ( 41)-( 43) can be used to evaluate the corrections to the streaming velocity v 2 s = u 2 s , v 2 s arising at order ε 2 .The solution was computed for different values of the channel undulation β and the Womersley number.Sample streamlines and associated vorticity contours are shown in Fig. 4 for β = 0.4 and two values of α.The solution, symmetric about x = 0.5, exhibits open streamlines aligned with the channel, whose walls delineate a stream tube carrying a nonzero flow rate.The flow structure is complicated by the presence of vortices with alternating circulation.For increasing α, these vortices approach the section of minimum width, with the velocity profile across the section of maximum width developing a nearly parabolic shape.
The variation with β and α of the net flow rate induced by the non-sinusoidal pressure difference of Fig. 2(a) was evaluated with use made of (44), giving the results shown in Fig. 5(a).The dependence on α 2 is seen to be non-monotonic, with Q increasing initially to reach a maximum at an intermediate value of α, beyond which Q decreases to eventually become negative.The dependence on β also is non-monotonic, with Q vanishing in the two limiting cases β ≪ 1, when the velocity approaches the Womersley solution with negligible convective acceleration, and 1 − β ≪ 1, when the pronounced wall undulation strangles the flow.According to the figure, values of in the range 0.7 < β < 0.8 result in the largest streaming motion for the specific periodic pressure considered in the computation.

Comparisons with direct numerical simulations
In order to validate the present theoretical model and assess its associated accuracy, the analytic predictions were compared with results of direct numerical simulations (DNS) employing the complete Navier-Stokes Eqs. ( 2)-( 4) to describe oscillatory flow with finite values of the stroke length ε in a slender channel with ℎ 0 /λ = 1/20 and total length 3λ.
The computations were carried out with the finite-volume solver Ansys Fluent (Release 20.2), assuring second-order accuracy in time and in space.A coupled algorithm was used for the pressure-velocity coupling.The computational domain was discretized using a structured uniform mesh.A grid sensitivity analysis was conducted to ensure the grid-size independence of the results.To that aim, a mesh initially containing 9×10 4 cells was systematically refined until no variation was obtained for the mean flow rate Q , in a channel with β = 0.3 and β = 0.4 for α = 6,8, 10 , resulting in a final configuration containing a total of 6×10 5 cells, which was used in all computations presented below.While the numerical evaluation of the integrals involved in the theoretical predictions takes seconds on a laptop computer, the DNS computations are significantly more expensive, with each simulation taking more than 60 h of CPU time to run 150 cycles, as needed to ensure convergence to a periodic solution with negligible inter-cycle variation, in a computational cluster using a total of 24 cores.
The converged numerical solution was used to compute the streaming flow by timeaveraging the periodic velocity.Predicted mean flow rates Q = ∫ 0 H u dy for ε = 0.1 and different values of α and β are represented as symbols in Fig. 5(a).As can be seen, the departures of the asymptotic predictions from the DNS results are very small for the conditions explored in the figure, thereby validating the theoretical development.Additional numerical results obtained for β = 0.4 and increasing values of ε are compared in Fig.
5 (b) with the asymptotic predictions Q /ε 2 = 0.016 α = 6 and Q /ε 2 = 0.031 α = 8 .The comparison reveals that the analytic results, which are strictly valid for asymptotically small values of ε ≪ 1, remain reasonably accurate for values of the stroke length exceeding ε = 0.5.This unexpected extended validity is consistent with previous findings pertaining to oscillating flow around circular cylinders [24].
Additional comparisons between the DNS and the theoretical predictions are shown in Fig. 6 for a channel with β = 0.4 and α = 6.Fig. 6(a) and (b) correspond to the steady-streaming velocity v 1 s = u 1 s , v 1 s and its correction v 2 s = u 2 s , v 2 s as obtained from the asymptotic analysis at order ε and ε 2 , respectively, using the driving pressure Π t = 3/4 cos t + 1/4 cos 2t shown in Fig. 2(a).The corresponding time-averaged velocity distribution obtained in DNS computations with ε = 0.1 is shown in Fig. 6(c).As can be seen in this last figure, the flow structure exhibits a slight aft-and-fore asymmetry, which is attributable to the multi-frequency character of the driving pressure.To further explore this issue, separate DNS computations were performed for the two sinusoidal signals Π t = 3/4 cos t and Π t = 1/4 cos 2t , with corresponding time-averaged velocities denoted by v I and v II , respectively.The sum of these two contributions, appropriately scaled with ε, is plotted in Fig. 6(d).The resulting symmetric flow structure is to be compared with the leading-order asymptotic prediction v 1 s shown in Fig. 6 (a).The results are nearly identical, with the peak stream function at the center of the vortices differing by less than 4%.
The observed differences between the DNS results in Fig. 6(c) and (d) are a consequence of the inter-mode interactions, which are responsible for the generation of the streamwise mean flow, described by the velocity function v 2 s = u 2 s , v 2 s in the asymptotic analysis at order ε 2 .These inter-mode interactions were quantified from the DNS results by plotting v − v I − v II in Fig. 6(e), with the scale ε 2 introduced for comparison with the velocity field v 2 s shown in Fig. 6(b).The small differences observed between the steady streaming resulting from intermode-interactions computed numerically (i.e. ( v − v I − v II /ε 2 ) and that predicted by the asymptotic analysis (i.e.v 2 s ), on the order of 10%, are compatible with the relative errors of order ε expected in the asymptotic description.Clearly, the satisfactory agreement displayed in Fig. 6 between the DNS results and the theoretical predictions provide additional confidence in the latter.

Cardiac and respiratory wave forms
The streaming flow correction v 2 s = u 2 s , v 2 s , identically zero in the sinusoidal case, depends fundamentally on the interactions of the different Fourier modes that define the shape of the periodic driving pressure Π t .In the asymptotic analysis, these interactions, which are apparent in (23), lead to the forcing term ℱ that ultimately determines the streaming flow ( 41)-( 44).Given the complexity of the associated algebra, it is not straightforward to ascertain how the interplay of the different modes dictates the final outcome.For example, for the two-mode pressure-wave form Π t = 3/4 cos t + 1/4 cos 2t , it was demonstrated in Fig. 5(a) that the resulting flow rate can be positive or negative depending on the value of α, but this finding might depend fundamentally on the pressure-wave form.Also, the linear dependence on α 2 revealed in (44) appears to indicate that pressure signals with a large high-frequency content (and therefore a higher effective Womersley number) might lead to more pronounced effects, but this has to be tested by investigating different waveforms, as done below by using as illustrative example the oscillatory motion in the spinal canal.
The oscillating flow of CSF is driven by the cardiac and respiratory cycles [25,26].
The cardiac-driven motion, induced by the pressure fluctuations generated in the cranial cavity by the inflow/outflow of arterial blood, is more pronounced in the cervical region [27,28], whereas the respiratory-driven motion, induced by abdominal pressure changes, is more pronounced in the lumbar region [29,30].The two periodic motions have very different periods, i.e. 2π/ω ≃ 1 s for the cardiac cycle and 2π/ω ≃ 5 s for the respiratory cycle.
Using the kinematic viscosity of CSF v = 0.7 × 10 −6 m 2 /s along with ℎ 0 ≃ 1 − 3 mm for the characteristic width of the spinal canal [31,32] yields corresponding Womersley numbers 3 ≲ α ≲ 9 for the cardiacdriven flow and 1.34 ≲ α ≲ 4.02 for the respiratory-driven flow.As revealed by magnetic resonance imaging (MRI), the resulting periodic motion is markedly non-sinusoidal, with typical normalized flow rates (instantaneous flow rate divided by its mean magnitude over a cycle) shown in Fig. 7.As indicated in the figure, the measurements of the cardiac-driven flow were taken at C3 level, whereas those of the respiratory-driven flow were taken at L1-L2 level, with results given for two different subjects.
In order to use the experimental measurements to investigate effects of complex waveforms, it is of interest to relate the shape of the oscillatory flow rate to the shape of the associated driving pressure difference (6).The computation of the coefficients A n involves the following steps.First, the flow-rate signal determined experimentally was normalized using its mean amplitude and the result was used to calculate the first 10 modes of the corresponding Fourier series Q 0 / Q 0 = ℜ ∑ n = 0 10 B n e int .The expression was then substituted into (19) to determine by inspection the corresponding Fourier coefficients of the pressure difference.The final values A n utilized in the computations were rescaled to give Π = 1, as is appropriate to enable comparisons between different waveforms.The resulting coefficients A n were computed for α = 6 (cardiac cycle) and α = 2.68 (respiratory cycle) with β = 0.3.Corresponding Π t functions are shown in separate panels of Fig. 7 along with the magnitude of the different coefficients A n .As can be seen, the pressure signal associated with the cardiac cycle exhibits pronounced peaks, also present in intracranial pressure measurements [12,33], whereas that of the respiratory cycle is more rounded.Correspondingly, the Fourier series of the cardiac-driven pressure signal shows a significant high-frequency content, which is not present in the respiratory-driven pressure, for which the higher-order amplitudes A 5 − A 10 are negligibly small.
To investigate the influence of the pressure waveform on streaming, the four different pressure signals Π t shown in Fig. 7 were used to compute the flow in the canonical channel-flow configuration considered here.The resulting secondorder steady-streaming velocity v 2 s = u 2 s , v 2 s is shown in the four panels on the right-hand side of Fig. 7, where a different scale is used for the stream function and vorticity of the cardiac and respiratory plots, as needed to accommodate their dissimilar magnitudes.As expected, since the Womersley number of the cardiac flow is larger, the associated streaming flow is much more vigorous, in agreement with the trends identified in Fig. 5(a).The general flow structure is similar to that depicted in Fig. 4, with the outer recirculating vortices of the cardiac cycle located closer to the section of minimum width, as can be expected from the previous results.There exist non-negligible inter-subject variability, with the first subject exhibiting much smaller cardiac-driven velocities but much larger respiratory-driven velocities.
For completeness, additional computations were carried out for different values of α using the pressure signals generated from the MRI flow-rate measurements, as delineated above.
In the range of Womersley numbers explored, corresponding to a canal width 1 ≤ ℎ 0 ≤ 3 mm with 2πω −1 = 1 s (cardiac) and 2πω −1 = 5 s (respiratory), the magnitude of the flow rate increases monotonically, as shown in Fig. 8.The large differences in Q revealed by the figure further underline the intersubject variability previously mentioned.As can be seen, the respiratory-driven flow rate of subject 1 is negative, but that of subject 2 is positive, indicating that relative changes in the frequency distribution of the non-sinusoidal signal modify the streaming flow in nontrivial ways that are difficult to predict.

Conclusion
The steady streaming associated with a waveform with multiple harmonics, previously analyzed in the context of the flow induced by oscillating cylinders and spheres [17][18][19], has been investigated here for wall-bounded flows by asymptotic and numerical methods, with the two-dimensional wavy-walled channel used as canonical model configuration.A perturbation analysis in the limit of small stroke lengths ε ≪ 1 provides closed-form analytic expression for the steady-streaming velocity v 1 s = u 1 s , v 1 s and its corrections v 2 s = u 2 s , v 2 s .
The latter, stemming from the interactions of the different Fourier modes that define the driving pressure gradient, carry a nonzero mean volumetric rate that breaks the fore-and-aft symmetry of the flow.The results of the asymptotic analysis are validated by comparisons with DNS computations.Somewhat unexpectedly, the asymptotic predictions for ε ≪ 1 are found to remain quantitatively accurate for values of ε > 0.5, well beyond the expected range of validity.Although a similar degree of agreement between asymptotic predictions and DNS computations has been found for ε 1 in other steady-streaming configurations, e.g. for flow about a linear array of cylinders [34], this result should be taken with caution, in that the relative departures can be dependent on the specific conditions.For instance, the agreement displayed in Fig. 5(b) might degrade for other values of α and β and also for pressure waveforms differing from that considered here, shown in Fig. 2(a).
The analysis indicates, in particular, that inter-frequency interactions leading to a net streamwise flow rate are more important for larger values of the Womersley number and intermediate values of the channel undulations.Sample computations using pressure gradients derived from in-vivo MRI measurements of the CSF flow rate in the spinal canal reveal large inter-subject variations.Also, the results suggest that the streaming flow resulting from inter-mode interactions is more pronounced for cardiac-driven flow, while it remains relatively small for respiratory-driven flow.Further steady-streaming studies addressing complex waveforms should account for the specific geometry of the spinal canal, thereby extending previous sinusoidal-pressure investigations [13,14,35].These future studies should investigate, in particular, whether the additional steady-streaming motion arising from inter-mode interactions can modify significantly the spatial structure of the mean Lagrangian flow, possibly providing connectivity between the different closed recirculating regions that have been predicted to appear along the spinal canal [31] In this context, effects of microanatomical features, such as nerve roots and ligaments, not accounted for in the previous work [31], also warrant further investigation.Schematic of the wavy channel of wave length λ, with an average channel width ℎ 0 and a relative amplitude of the wall undulation β.Streamlines and colour contours of vorticity corresponding to the second-order steadystreaming velocity v 2 s driven by a pressure difference Π t = 3/4 cos t + 1/4 cos 2t in a canal with β = 0.4 for α = 5 and α = 10.Fixed constant streamline spacings δψ 2 are used for the plots, with δψ 2 = 5 ⋅ 10 −4 for α = 5 and δψ 2 = 3 ⋅ 10 −3 for α = 10.Temporal of the normalized rate Q 0 / Q 0 induced by the cardiac cycle in the cervical region [31] and by the respiratory cycle in the lumbar region 2(a) are shown in Fig. 2(b).

Fig. 5 .
Fig. 5.The mean flow rate Q as obtained from the asymprotic analysis for ε&1 (curves) and from the DNS computations in a channel with ℎ a /λ = 1/20 (symbols).The plors show (a) the variation of Q / ε 2 α 2 for increasing α and different values of β (DNS results computed with = 0.1) and (b) a comparison for β = 0.4 of the asymprotic predictions Q /ε 2 = 0.0306 α = 8 and Q /ε 2 = 0.0156 α = 6 with the values obtained numerically for increasing values of ε.
[30]  as obtained for two different healthy subjects from magnetic resonance measurements, The pressure signals Π and associated Fourier coefficients A n shown in the accompanying panels are determined using (19) for β = 0.3 with α = 6 (cardiac) and α = 2.68 (respiratory), with the resulting signal normalized to give Π = 1.The panels on the right-hand side of the figure represent the streamlines and colour contours of vorticity corresponding to the second-order steady-streaming velocity v 2 s driven in the wavy-walled channel by the four pressure signals Π. Fixed constant streamfunction increments δψ 2 are used to draw the streamlines, with δψ 2 = 1 ⋅ 10 −2 for cardiac-driven flow and δψ 2 = 5 ⋅ 10 −5 for respiratory-driven flow.

Fig. 8 .
Fig. 8.The variation with α of rescaled mean flow Q /ε 2 determined for β = 0.3 using the cardiac and respiratory pressure signals derived from the normalized flow rates shown in Fig.7.