2.1. Model and equations

We consider a viscous Newtonian fluid with constant properties in a straight smooth rigid pipe of circular cross-section with a time-dependant bulk velocity (1.1). The flow is assumed to be incompressible and governed by the dimensionless Navier–Stokes equations (NSE)

(2.1)\begin{equation} \frac{\partial\pmb{u}}{\partial t} + \left(\pmb{u} \boldsymbol\cdot \boldsymbol{\nabla} \right) \pmb{u} ={-}\boldsymbol{\nabla} p + \frac{1}{{Re}}\nabla^{2}\pmb{u}+F_{d}\left(t\right) \pmb{e}_{z} \quad\text{and}\quad \boldsymbol{\nabla}\boldsymbol\cdot\pmb{u} = 0 . \end{equation}

Here, $\pmb {u}$ is the fluid velocity, $p$ the pressure and $F_{d}(t)$ the axial pressure gradient which drives the flow at the desired bulk velocity defined in (1.1). All variables in this study are rendered dimensionless using the pipe diameter ($D$), the time-averaged bulk velocity ($u_{s}$) and the fluid density ($\rho _{f}$). The equations are formulated in a cylindrical coordinate system $(r,\theta,z)$, where the velocity field has three components $\pmb {u}=(u_{r},u_{\theta },u_{z})$ in the radial, azimuthal and axial direction, respectively.

The laminar SW profile is unidirectional and only depends on the radial position and time $U\equiv u_z(r,t)$. In this work we refer to the laminar profile $U$ as the linear superposition of the SW profiles coming from the multiple harmonic components of (1.1).

The linearised Navier–Stokes equations (LNSE)

(2.2)\begin{equation} \frac{\partial \pmb{u}'}{\partial t} + \left(U \boldsymbol\cdot \boldsymbol{\nabla} \right) \pmb{u}' + \left(\pmb{u}' \boldsymbol\cdot \boldsymbol{\nabla} \right) U ={-}\boldsymbol{\nabla} p' + \frac{1}{{Re}}\nabla^{2}\pmb{u}' \quad\text{and}\quad \boldsymbol{\nabla}\boldsymbol\cdot\pmb{u}' = 0 \end{equation}

are obtained by decomposing the full velocity ($\pmb {u} = U\pmb {e}_{z} + \pmb {u}'$) and pressure ($p=P+p'$) fields into a laminar base flow (capital letters) and infinitesimal perturbations (prime superscript).

2.2. Bulk velocity waveforms

We design a generic waveform for the bulk velocity ($u_b$) to explore the effect its shape has on the stability of the corresponding laminar velocity profile. Brindise & Vlachos (Reference Brindise and Vlachos2018) recently showed experimentally that the slope and duration of the acceleration and deceleration phases are important for turbulence transition. Here we define waveforms whose characteristics we can control systematically. Specifically, we propose a certain way of defining the waveform by fixing only three parameters.

We first define six control points (black stars in figure 1), which represent the skeleton of our generic waveform. Their position is controlled by three parameters ($t_{ac}$, $t_{dc}$ and $t_{m}$). These parameters set the value of the coefficients $a_{n}$ and $b_{n}$ in (1.1) and do not affect the mean velocity ($u_{s}$) that defines Re. Then we define a spline, using a monotone piecewise cubic Hermite interpolating polynomia (Fritsch & Carlson Reference Fritsch and Carlson1980), that captures the position of the control points. We finally fit the spline using $N_{F}=30$ Fourier modes to obtain a smooth and periodic pulsation. In figure 1 we show the temporal evolution of $u_b$ and the corresponding laminar velocity profiles for a sine wave pulsation compared with four other cases. Namely, waveform 1 (WF1) and waveform 2 (WF2), the two base cases we use throughout our analysis, and waveform 3 (WF3) and waveform 4 (WF4), the two cases with the longest/shortest $t_{m}$ we consider in this work. We enforce an additional constraint to all waveforms investigated here, so that the bulk velocity never falls below $u_{b}=0$. This is inspired by cardiovascular flows, where the minimum bulk velocity in the larger vessels is close to zero (Bürk et al. Reference Bürk, Blanke, Stankovic, Barker, Russe, Geiger, Frydrychowicz, Langer and Markl2012). Although the bulk velocity is always positive, locally the velocity profile can have negative axial velocities (i.e. flow reversal) during some fraction of the pulsation, as exemplarily shown in figures 1(c) and 1(d).

Figure 1.

Definition of the generic waveform (WF) including five examples for $ {Wo}=11$. (a) Temporal evolution of the bulk velocity ($u_b$) over one pulsation period ($T$). The black stars denote the six control points, which define the waveform and the solid lines represent the 30 Fourier mode approximation of the spline that passes through those points. The five examples are WF1 ($t_{ac}=t_{dc}=0.05$, $t_{m}=0.45$), WF2 ($t_{ac}=t_{dc}=0.2$, $t_{m}=0.55$), WF3 ($t_{ac}=t_{dc}=0.1$, $t_{m}=0.4$), WF4 ($t_{ac}=t_{dc}=0.05$, $t_{m}=0.6$) and a single harmonic sine wave pulsation with $A=1$ as a reference. (b–d) The corresponding (laminar) velocity profiles ($U$) of the five waveforms defined in (a) at three different instants in time. Filled circles denote the existence and radial location of inflection points in the velocity profile.

All the waveforms we consider have an acceleration phase with a slope that is set by the parameter $t_{ac}$. Note that the total duration of the acceleration is $2 t_{ac}$ long. The bulk velocity remains in a high-velocity phase for the time span $t_{m}-t_{ac}-t_{dc}$. Then the pulsation enters a deceleration phase, whose slope is set by the parameter $t_{dc}$, so that the total duration of decelerates phase is $2 t_{dc}$ long. Finally the bulk velocity remains in a low-velocity phase for the rest of the period ($T-t_{m}-t_{ac}-t_{dc}$). The parameter $t_{m}\in [0,T]$ sets the maximum $ {Re}_{max}$ of the flow as

(2.3)\begin{equation} {Re}_{max} = {Re}\frac{T}{t_{m}}. \end{equation}

For $t_{m}={T}/{2}$ the waveform is symmetric and the high- and low-velocity phases have the same duration. For this specific choice, $ {Re}_{max}=2 {Re}$ and the minimum velocity is $u_{b}=0$. As $t_{m}\rightarrow 0$ the time the flow stays in a high-velocity phase decreases and $ {Re}_{max}$ increases. In the following $t_{ac}$, $t_{dc}$ and $t_{m}$ are normalised in terms of $T$ and they must satisfy

(2.4a–c)\begin{equation} t_{m}+ t_{ac} + t_{dc} < 1 ,\quad t_{ac}<0.5 \quad\text{and}\quad t_{dc}<0.5 . \end{equation}

2.3. Stability analysis

We study the linear stability of the laminar velocity profile ($U$) with two different approaches. First, we freeze $U$ at 200 equispaced instants in time and perform instantaneous LSA using the numerical method of Meseguer & Trefethen (Reference Meseguer and Trefethen2003). The method returns the complex eigenvalues $\pmb {\lambda }$ for each instantaneous velocity profile. If one of them has a positive real part, the profile is considered to be instantaneously linearly unstable. To measure the instantaneous level of instability we compute the maximum real part of the eigenvalues

(2.5)\begin{equation} \lambda_{max}\left(t\right)=\max \left( \mathbb{R} \left\{\pmb{\lambda} \left(t\right) \right\}\right) \end{equation}

for the instantaneous $U(r,t)$ at each discrete time step. Strictly, this approach is only valid in the quasi-steady limit, where $U(r,t)$ evolves much slower than the velocity perturbations ($\pmb {u}'$). It is, however, not a priori clear, how far these two time scales have to be apart for the quasi-steady assumption to remain approximately valid.

Second, we carry out a non-modal stability TGA. With this method we solve an optimisation problem to determine the perturbation with highest energy growth ($G$) out of all possible perturbations, in terms of perturbation shape ($\pmb {u}'_{0}=\pmb {u}'(t_{0})$) and initial time of perturbation (${t_{0}}/{T}$). In contrast to the LSA approach, this method allows both the laminar profile and the perturbations to evolve in time. The TGA yields the optimal transient energy growth

(2.6)\begin{equation} G_{TGA} = \max \left( \frac{E\left(t \right)}{E\left(t_{0} \right)} \right) = \max\left(\frac{\pmb{u}'\left(t\right)\boldsymbol\cdot \pmb{u}'\left(t\right) }{ \pmb{u}'\left(t_{0}\right) \boldsymbol\cdot \pmb{u}'\left(t_{0}\right) }\right) , \end{equation}

where $E= \pmb {u}'\boldsymbol {\cdot } \pmb {u}'$ is the perturbation energy. We refer the reader to Xu et al. (Reference Xu, Song and Avila2021) for further details.

For both approaches, we discretise (2.2) using a Fourier–Galerkin ansatz in $\theta$ and $z$, and a Chebyshev collocation method in $r$. For both we use $N_{r}=96$ radial points and consider only some azimuthal ($m$) and axial ($k$) wavenumbers. Our goal is to study the behaviour of helical perturbations. Thus, we only consider $m=1$ and $2\leq k \leq 4$ in steps of $\Delta k=0.5$, which correspond to the most amplified helical perturbations reported by Xu et al. (Reference Xu, Song and Avila2021). We set $\Delta t=0.0025 ({D}/{u_{s}})$ for both the integration of the laminar SW profile in the LSA and TGA, and the integration of perturbations in our TGA using (2.2). The used grid size and time step size are the same as in Xu et al. (Reference Xu, Song and Avila2021), who found them to be sufficient for a $ {Re}_{max}$ one order of magnitude higher than those we consider here.

2.4. DNS

We perform DNS of (2.1) using our open-source pseudo-spectral simulation code nsPipe (available at https://github.com/dfeldmann/nsCouette, López et al. Reference López, Feldmann, Rampp, Vela-Martín, Shi and Avila2020). In nsPipe, the governing equations (2.1) are discretised in cylindrical coordinates $(r,\theta,z)$ using a Fourier–Galerkin ansatz in $\theta$ ($N_{\theta }$ modes) and $z$ ($N_{z}$ modes) and high-order finite differences in $r$. Periodic boundary conditions are imposed in $\theta$ and $z$, and no-slip boundary conditions in the solid pipe wall. The discretised NSE are integrated forward in time using a second-order predictor–corrector method with variable time-step size ($\Delta t$). Further details about the numerical methods and functionalities of nsPipe are given in López et al. (Reference López, Feldmann, Rampp, Vela-Martín, Shi and Avila2020) and references therein.

Here, we perform DNS with different waveforms at $ {Re}=2000$ and $ {Wo}=11$ in a computational domain of size $({D}/{2}\times 2{\rm \pi} \times {50}D)$ using a resolution of $(N_{r}\times N_{\theta }\times N_{z})=(96\times 192\times 3000)$. For all waveforms we consider here, the largest instantaneous friction Reynolds number we have measured is $ {Re_{\tau }}=188$. In that case, the used resolution corresponds to a grid spacing in viscous units of $0.04 \le \Delta r^+\le 2.82$, $R^+\Delta \theta = 6.2$ and $\Delta z^+=6.3$, respectively. The time step size is automatically adapted to ensure numerical stability and accuracy and is always between $\Delta t=0.0009({D}/{u_s})$ and $0.0025({D}/{u_s})$. The code adjusts the value of $F_{d}(t)$ in (2.1) to enforce the desired bulk velocity $u_{b}(t)$, (1.1), at each time step.

3.1. Mechanism of the perturbation growth

We first performed an LSA for a single harmonic sine wave pulsation at $ {Re}=2000$, $ {Wo}=11$ and $A=1$ and observed that the velocity profile is instantaneously unstable for more than 50 % of the period. This can be seen in figure 4, where we show the maximum real part out of all the instantaneous eigenvalues ($\lambda _{max}$, see (2.5)) for this case. For most of the acceleration phase, $\lambda _{max}$ is constant and negative. This corresponds to the maximum eigenvalue of Hagen–Poiseuille flow at $ {Re}=2000$ (Meseguer & Trefethen Reference Meseguer and Trefethen2003). However, for the second half of the deceleration phase and the first half of the acceleration phase, $\lambda _{max}$ is positive. The crossover occurs at ${t}/{T}\approx {0.45}$, which is very close to the optimal time to trigger the helical perturbation (${t_0}/{T}\approx {0.5}$) found by Xu et al. (Reference Xu, Song and Avila2021) based on a TGA for the same values of Re, Wo and $A$. These results suggest that the helical perturbation actually takes advantage of the instantaneous linear instability of the laminar velocity profile.

Figure 4.

Laminar profile and instantaneous maximum eigenvalue $\lambda _{max}$ according to our LSA for a sine wave pulsation. In yellow the instantaneous laminar profiles $U(r,t)$ at $ {Re}=2000$, $ {Wo}=11$ and $A=1$. To not interfere with one another the profiles are scaled using a scalar with arbitrary units so the all time maximum is smaller than $t/T={0.15}$, because only the development of $U(r,t)$ in time is of interest. Circles denote the radial position ($r_{i}$) of inflection points in the velocity profile. Filled points correspond to inflection points that also satisfy the Fjørtoft criterion locally $({\partial ^2 U}/{\partial r^2})(U(r,t)-U(r_{i} ) )<0$. In red, the maximum real component out of all the instantaneous eigenvalues of the laminar profile is shown.

According to Miau et al. (Reference Miau, Wang, Jian and Hsu2017) and Nebauer (Reference Nebauer2019), the instantaneous linear instability of the SW velocity profile is related to the existence and characteristics (number or position ($r_{i}$)) of inflection points (${\partial ^2 U}/{\partial r^2}=0$). An inflection point is regarded as inviscidly unstable, when the Fjørtoft criterion

(3.1)\begin{equation} \frac{\partial^2 U}{\partial r^2}\left(U-U\!\left(r_{i}\right)\right)<0 \end{equation}

is satisfied locally (Schmid et al. Reference Schmid, Henningson and Jankowski2002). Nebauer (Reference Nebauer2019) already discussed that perturbations can sit on top of these inflection points and feed energy from them. In the following, we show how the helical perturbations take advantage of this mechanism to grow. To this end, we integrated the linearised equations (2.2) forward in time using the optimal helical perturbation according to our TGA as the initial condition. We computed the production ($P^\prime$) and dissipation ($D^\prime$) of the kinetic energy ($E$) contained in the perturbations as

(3.2)\begin{equation} P^{\prime} ={-}u_{r}^{\prime} u_{z}^{\prime}\frac{\partial U}{\partial r} \quad\text{and}\quad D^{\prime}={-}\frac{1}{{Re}}\boldsymbol{\nabla}\boldsymbol{u}^{\prime}:\boldsymbol{\nabla}\boldsymbol{u}^{\prime} \quad\text{where}\ \frac{\mathrm{d}E }{\mathrm{d} t} =P^{\prime} + D^{\prime} . \end{equation}

Note that $u_{r}^{\prime }$ and $u_{z}^{\prime }$ of the helical perturbation have an axial wavenumber of $k$, and the product between the two results in structures with axial wavenumber $2k$.

In figure 5 (and the Supplementary Movie available at https://doi.org/10.1017/jfm.2022.681) we show colour maps of $P^{\prime }$ and $D^{\prime }$ (only in Movie 1) structures in a $r$–$z$-plane, and on top, the instantaneous laminar velocity profile. Strong events of production and dissipation clearly follow the radial position of inflection points in the SW profile as time marches. At each time step, the strongest production events are always bigger than the strongest dissipation events (see the Supplementary Movie). This difference between the magnitude of production and dissipation structures explains the large growth of helical perturbations. Moreover, the fact that strong production events follow the inflection points, further suggests that the helical perturbations extract their growth from them.

Figure 5.

Link between inflection points in the SW profile ($U$) and production ($P^{\prime }$, see (3.2)) of kinetic energy contained in the helical perturbations in the $r$–$z$-plane at $\theta =0$. Results correspond to an integration of (2.2) using the optimal helical perturbation as initial condition for a sine wave pulsation at $ {Re}=2000$, $ {Wo}=11$ and $A=1$. Yellow lines represent $U\!(r,t)$ scaled in arbitrary units and circles represent existence and location ($r_{i}$) of inflection points. Yellow circles additionally satisfy the Fjørtoft criterion locally. We show two different instants of time: (a) mid deceleration phase at ${t_{0}}/{T}={0.5}$; (b) mid acceleration phase at ${t}/{T}={1}$. In both the production $P^{\prime }$ is normalised by the maximum value at ${t_{0}}/{T}={0.5}$, where the simulation was started.

3.2. Simple model for perturbation growth

Pulsatile pipe flow has at least two important time scales when it comes to the evolution of perturbations. One is the advective time scale (${D}/{u_s}$) and the other one is the pulsation period ($T={{\rm \pi} {Re}}/{2 {Wo}^2}$ in advective time units). As Cowley (Reference Cowley1987) mentioned, for sufficiently long periods (in terms of ${D}/{u_s}$), the perturbations would see a quasi-steady velocity profile. In that case, the perturbations would have enough time to grow on top of the instantaneous linear instability before the velocity profile changes a lot and becomes stable again.

In view of these findings, we propose that the energy growth ($G_{TGA}$) observed by Xu et al. (Reference Xu, Song and Avila2021), depends on how much and how long the SW velocity profile is linearly unstable. The instability of the velocity profile in turn depends on the existence of inflection points that satisfy the Fjørtoft criterion (Schmid et al. Reference Schmid, Henningson and Jankowski2002; Nebauer Reference Nebauer2019; Kern et al. Reference Kern, Beneitez, Hanifi and Henningson2021), as already discussed previously. With these ideas in mind we propose that from an LSA perspective the energy growth rate should scale like

(3.3)\begin{equation} \frac{E_{max}}{E_{0}}\propto G_{LSA}=e^{2\,\lambda_{i}\,T}, \end{equation}

where $\lambda _{i}$ is the time integral

(3.4)\begin{equation} \lambda_{i}=\frac{1}{T}\int_{t_{0}}^{t_{0}+\Delta t_{u}}\lambda_{max}\left(t\right)\,\mathrm{d}t, \end{equation}

for the time window $\Delta t_{u}$ where $\lambda _{max}>0$. The new parameter $\lambda _{i}$ is taken as a combined proxy for how much and how long the laminar profile $U$ is linearly unstable during one pulsation period.

In the following section, we compare the energy growth from our TGA with the hypothesis we propose in (3.3). We show that both calculations yield similar results, which indicates that the energy growth reported by Xu et al. (Reference Xu, Song and Avila2021) is not due to non-modal mechanisms, but to the instantaneous linear instability in the laminar profile.

4.1. Dependency with respect to the waveform

We first focus on $t_{m}$, which controls the asymmetry of the waveform. The smaller $t_{m}$ is, the shorter the high-velocity phase (in terms of $T$) and the larger $ {Re}_{max}$ become. From figure 6(a) it is clear that $\lambda _{i}$ increases monotonically as $t_{m}$ decreases; i.e. the waveform goes from a long high-velocity phase to a long low-velocity phase. This is because a longer low-velocity phase (smaller $t_m$) results in a longer fraction of the period $\Delta t_u$ where the profile is instantaneously unstable (figure 6b). Thus, the shorter $t_{m}$ is, the more unstable the laminar profile becomes. This conclusion is in good agreement with experimental findings of Brindise & Vlachos (Reference Brindise and Vlachos2018), who showed that flows with longer deceleration and longer low-velocity phases are more prone to transition by looking at turbulent kinetic energy and its production.

Figure 6.

(a) Relationship between eigenvalue proxy $\lambda _{i}$ (see (3.4)) and $t_{m}$. Cases correspond to $ {Re}=2000$, $ {Wo}=11$ and different lines indicate different $t_{ac}$ and $t_{dc}$. The lines correspond to the legend in panel (b). (b) Plot of $\Delta t_{u}$ or fraction of the period during which the laminar profile is instantaneously unstable for different $t_{m}$, $t_{ac}$ and $t_{dc}$. (c) Plot of $\Delta t_{u}$ for different waveforms with respect to Wo. The thin grey line corresponds to the lifetime of inflection points that satisfy the Fjørtoft criterion $\Delta t_{i}$ of WF2. (d) In colour $\Delta t_{i}$, and inflection point position span $\Delta r_{i}$ as an area, with respect to Wo for WF2.

In figure 6(a), we also show how $\lambda _{i}$ depends on the other two waveform parameters. Recall that $t_{ac}$ and $t_{dc}$ control the slope of the acceleration and the deceleration, but do not affect $ {Re}_{max}$. It is evident from figure 6(a), that $\lambda _{i}$ is inversely proportional to both parameters, implying that steeper acceleration and steeper deceleration both lead to more unstable flows. However, the sensitivity of $\lambda _{i}$ with respect to $t_{dc}$ is larger than the sensitivity with respect to $t_{ac}$. While increasing $t_{ac}$ by a factor of four decreases $\lambda _{i}$ by only 9 %, doing the same for $t_{dc}$ decreases $\lambda _{i}$ by 16 %.

4.2. Dependency with respect to Reynolds and Womersley numbers

Our hypothesis is that $G_{LSA}$ is proportional to the product of the period ($T$) and the proposed eigenvalue proxy ($\lambda _{i}$, (3.4)). Here we explore the dependency of $\lambda _{i}$ on Re and Wo and we compare $G_{LSA}$ with $G_{TGA}$ obtained by our TGA.

For both waveforms we have considered here, $\lambda _{i}$ grows with the Reynolds number (figure 7a). In the inviscid regime ($ {Re}\rightarrow \infty$), $\lambda _{i}$ approaches a given asymptotic value, which depends on the waveform.

Figure 7.

(a) Eigenvalue proxy $\lambda _{i}$ with respect to Re at $ {Wo}=11$ and two different waveforms (blue/green lines). (b) Eigenvalue proxy $\lambda _{i}$, (3.4), with respect to Wo at $ {Re}=2000$. (c) Energy growth $G_{LSA}$, see (3.3), with respect to Re at $ {Wo}=11$. (c) Energy growth $G_{LSA}$ with respect to Wo at $ {Re}=2000$. Orange lines correspond to the optimal transient growth $G_{TGA}$. Blue and vermillion lines correspond to waveform 1 with $t_{m}=0.45$ and $t_{ac}=t_{dc}=0.05$, whereas green lines correspond to waveform 2 with $t_{m}=0.55$ and $t_{ac}=t_{dc}=0.2$.

The dependence of $\lambda _{i}$ on the Womersley number is more complex. From $ {Wo}\gtrsim {2}$ onwards, $\lambda _{i}>0$ (figure 7b), but the exact value at which the flow becomes unstable, depends on both the Reynolds number and the waveform. Thereafter $\lambda _{i}$ increases with Wo, until it reaches a maximum around $ {Wo}\approx {11}$. The exact $ {Wo}$ and magnitude of this maximum depend on the waveform. If Wo is now further increased, $\lambda _{i}$ decreases again but remains positive for all parameters considered here.

The dependence of $\lambda _{i}$ on Wo is determined by the relationship between $\Delta t_{u}$ and Wo, as observed when comparing the curves in figures 7(b) and 6(c). It is the fraction of the period where the flow is unstable, $\Delta t_{u}$, what dictates the value of $\lambda _{i}$ with respect to Wo. In turn, as shown in figure 6(c), $\Delta t_{u}$ follows the trend of $\Delta t_{i}$. Here, $\Delta t_{i}$ is the fraction of the period, where the profile has only one inflection point, that additionally satisfies the Fjørtoft criterion. Thus, the presence of inflection points sets the fraction of the period where the laminar profile is unstable, which ultimately sets the level of instability $\lambda _{i}$.

For $ {Wo}\gtrsim {3}$ the velocity profiles exhibit inflection points for more than one-quarter ($\Delta t_{i} \gtrsim {T}/{4}$) of the pulsation period (figure 6c,d). With increasing Womersley number (${3}\lesssim {Wo}\lesssim 17$), the lifespan of the inflection points ($\Delta t_{i}$) increases. The inflection points appear close to the pipe wall at the early stages of deceleration (figure 4). During the rest of the deceleration and the subsequent low-velocity phase, the inflection points move towards the pipe centreline. However, before they are able to reach the centreline, they disappear during the acceleration phase. Their movement is restricted to a radial span $\Delta r_{i} = \max (r_{i})-\min (r_{i})$ that decreases with increasing Womersley number (figure 6d). For large Wo, the evolution of the velocity profile prevents the inflection points from approaching the centreline before they die. Already from $ {Wo}\approx 11$ onwards, they remain in the vicinity of the pipe wall ($\min (r_{i})< {D}/{4}$) and so do the perturbations that may grow on top of them. This in turn does not allow perturbations to access the more energetic flow in the central region of the pipe, resulting in a decreasing, but still bigger than zero, $\lambda _{i}$ for $ {Wo}>11$ (figure 7a). Note that at $ {Wo}>11$ perturbations can still extract energy from the inflection points but less efficiently than at $ {Wo} \approx 11$.

In order to characterise the dependency of $G_{LSA}$, (3.3), with respect to Re and Wo we combine our knowledge of $\lambda _{i}$ with the effect of the pulsation period $T={{\rm \pi} {Re}}/{2 {Wo}^{2}}$. We observe that for intermediate Womersley numbers, $G_{LSA}$ grows monotonically with Re (figure 7c). At sufficiently high Reynolds numbers, $\lambda _{i}$ is more or less constant (figure 7a) and $G_{LSA}$ ends up following the exponential relationship between $T$ and Re.

The combined effects of $\lambda _{i}$ and pulsation period set the point of maximum growth at $ {Wo}\approx 7.5$ (figure 7d). Depending on the waveform or the Reynolds number, the exact position of this maximum with respect to Wo can vary slightly. Interestingly the region close to $ {Wo}\approx 7.5$ matches the point of maximum transient growth for a flow driven with a sine wave pulsation (Xu et al. Reference Xu, Song and Avila2021), and it is close to the point of maximum growth of perturbations in pulsatile channel flow (Pier & Schmid Reference Pier and Schmid2017). It is at this particular Wo, where the competing effects of shorter pulsation periods (in terms of advective time units) and higher level of average instability of the laminar profile make the flow more susceptible for perturbations to grow.

In figures 7(c) and 7(d) we show that $G_{LSA}$, (3.3), approximates the optimal transient growth $G_{TGA}$, (2.6), reasonably well at several Re and Wo. Note that we do not use any constant to match the two lines, but directly show $G_{LSA}$ as computed according to (3.3). This shows that the energy growth of the helical perturbation is related to the instantaneous instability of the laminar profile and confirms our hypothesis. It is the instantaneous instability of the laminar profile what yields the outstanding perturbation growth observed in the TGA of Xu et al. (Reference Xu, Song and Avila2021), as long as $T$ is sufficiently long (in terms of advective time units). In our results, $G_{TGA}>G_{LSA}$, as expected. This means that, on top of the modal growth that comes from the instantaneous instability, the perturbation can further grow due to additional non-modal mechanisms.

4.3. Gradient descent

In this section, we exploit the knowledge gained on perturbation growth so far, to develop a simple parametric model. Specifically we model the dependency of the perturbation growth on the governing parameters by approximating (fitting) our two sets of LSA and TGA computations listed in table 1 with the expression

(4.1)$$\begin{gather} \log G_{g} = s \sigma {Re} + bi_{3}, \end{gather}$$

(4.2)$$\begin{gather}s = bi_{2} + w_{m} t_{m} + w_{ac} t_{ac} + w_{dc} t_{dc}, \end{gather}$$

(4.3)$$\begin{gather}\sigma = [w_{1} \left( {Wo}-b_{1}\right) + w_{2}\left( {Wo} -bi_{1}\right)^{2}] \exp{\left({-}w_{3}{Wo}\right)}. \end{gather}$$

The parameter $G_{g}$ is our approximation to $G_{LSA}$ or $G_{TGA}$. We use an exponential dependence on Re, an assumption motivated by (3.3), where we suggest that the perturbation growth scales exponentially with the product of the pulsation period $T$ and $\lambda _{i}$. As we show in figure 7(a), at a sufficiently high $ {Re}$, $\lambda _{i}$ reaches a constant value and $G_{LSA}$ (and, therefore, $G_{g}$) follows the exponential relationship with $ {Re}$ that comes from $T={{\rm \pi} {Re}}/{2 {Wo}^2}$ (figure 7b). Further, we assume that the slope of this relationship is given by the product of (4.2) and (4.3). The function $\sigma$ approximates the shape of $G_{LSA}$ with respect to Wo that is shown in figure 7(d). We now know that at $ {Wo}\approx 11$, $\lambda _{i}$ peaks, whereas $T$ always decreases as $ {Wo}$ increases. From these two ideas, we select a function $\sigma$ that goes to zero as $ {Wo}\rightarrow 0$ and $ {Wo} \rightarrow \infty$. The function $s$ accounts for the dependency of $G_{g}$ on the shape of the pulsation waveform (i.e. $t_{m}$, $t_{ac}$ and $t_{dc}$). From figures 6(a) and 6(b), $\lambda _{i}$ appears to depend linearly on each of the three parameters (at least for the values we consider here) and, henceforth, $G_{LSA}$ depends exponentially on each of them. Although the results in previous sections suggest cross-dependencies between the parameters, we ignore these in our model equation for the sake of simplicity.

We approximate $G_{g}$ by looking for the set of weights ($w_{i}$) and biases ($bi_{i}$) in (4.1)–(4.3) that minimise the error

(4.4)\begin{equation} \epsilon=\frac{1}{N} \sum^{N}_{n=1} \left(\log G-\log G_{g}\right)_{n}^2, \end{equation}

where $N$ is the total number of data items to fit (table 2). We produce two fits, one to the LSA results where $G=G_{LSA}$ and another to the TGA results where $G=G_{TGA}$. We initialise each fit with the vector $\pmb {x}^{0}$ that is filled with random guesses of weights and biases with values between zero and one. Then we use a gradient descent method to find the vector $\pmb {x}^{i}$ that minimises (4.4). Iterations are performed until a minimum is reached. At iteration $i$ the weights and biases are updated as

(4.5)\begin{equation} \pmb{x}^{i+1}=\pmb{x}^{i} - \eta \sum^{N}_{n=1} \frac{\mathrm{d}\epsilon^{n}}{\mathrm{d}\kern0.7pt \pmb{x}^{i}} , \end{equation}

where $\eta$ is a learning parameter that is dynamically adjusted so that the error $\epsilon (\pmb {x}^{i+1})<\epsilon (\pmb {x}^{i})$. We consider the case as converged when the error decreases to less than $10^{-10}$ for consecutive iterations.

Table 2.

Weights ($w$) and biases ($bi$) used to fit LSA and TGA results to (4.1)–(4.3).

The quality of the fit is visualised in figure 8. The horizontal axis represents the number of the case, where the list of cases is ordered in the sense of increasing first Re, then Wo and finally $t_{m}$, $t_{ac}$ and $t_{dc}$ as in table 1. This explains why the data appear in packets of functions that look similar to the shape of the function shown in figure 7(c).

Figure 8.

Results of our fit using the gradient descent method. (a) Fit of $G_{LSA}$ (see (3.3))–(4.1). Hollow circles are individual results of the LSA for all the parameters considered (first row in table 1). Black dots are the fit of the method. In each horizontal location we plot just one case with a given combination of Re, Wo, $t_{m}$, $t_{ac}$ and $t_{dc}$. The cases are ordered in increasing sense of first Re then Wo and finally $t_{m}$, $t_{ac}$ and $t_{dc}$. (b) Same fit but for $G_{TGA}$ (second row in table 1).

Note that the fit performs poorer as Re increases. We know that Re also has an effect on $\lambda _{i}$, as shown in figure 7(b). However, the formulation we have proposed ignores this, and other cross-dependencies, as it only considers the dependency of $G$ on Re coming from $T$. Despite these differences, the error is of the order of $10^{-4}$ for both fits.

The final weights of the two fits (table 2) show that $t_{m}$ has a greater effect on the energy growth compared with $t_{ac}$ and $t_{dc}$, because $|w_{m}|>|w_{dc}|>|w_{ac}|$. The smaller $t_{m}$, $t_{ac}$ and $t_{dc}$ are, the bigger the energy growth is, because $w_{m}< w_{ac}< w_{ac}<0$. Note that the two fits do not reach the same values for all weights and biases. This is due to the differences between the two sets of data, and the choices we use for the fitting formula. The weights are the key parameters, as they quantify the importance of some parameters with respect to others, and are found to be very similar for the fits to $G_{LSA}$ and $G_{TGA}$.

The fits allow us to generalise the dependency of $G_{LSA}$ and $G_{TGA}$ on all the control parameters for all the waveforms we have considered here. With this tool one can infer how much a helical perturbation can grow in a given pulsatile pipe flow by only knowing the waveform of the pulsation and the control parameters Re and Wo.

4.4. Physiological waveform

As a test to our observations, we study the behaviour of a laminar profile driven by a physiological waveform. To that end a particular signal presented in the physiological study of Bürk et al. (Reference Bürk, Blanke, Stankovic, Barker, Russe, Geiger, Frydrychowicz, Langer and Markl2012) was selected (figure 9a). In their study, they measured the mean velocity of blood flow at different sections of the aorta in several patients. In our study the signal for the descending aorta section of a young volunteer was fitted using $N_{F}=8$ Fourier modes. The resultant Fourier mode coefficients for the physiological waveform (1.1) are presented in table 3. We compute the waveform parameters by measuring the time span between the half point of the acceleration and deceleration ($t_{m}\approx 0.279$). In addition, we measure half the length of the acceleration ($t_{ac}\approx 0.108$) and the deceleration ($t_{dc}\approx 0.141$). By introducing these parameters to expression (4.1) and using the weights listed in table 2, we generated a guess of the energy growth on top of the laminar velocity profile in the physiological waveform (figure 9b,c). In addition, we performed LSA (solid blue line) and TGA (red dots) to probe the accuracy of the fits.

Figure 9.

Results and fit of the energy growth of helical perturbations on a laminar profile driven with a physiological waveform. (a) Physiological waveform defined by coefficients of table 3 in blue, compared with a sine wave pulsation in yellow. (b) Evolution of energy growth $G$ with Re at $ {Wo}=12$ for the physiological waveform. The blue solid line corresponds to $G_{LSA}$ and the blue dotted line is a guess using the expression (4.1) with the weights in the first row of table 2. Filled vermillion points correspond to TGA results and the dotted vermillion line corresponds to a guess using the expression (4.1) with the weights in the second row of table 2. (c) Evolution of energy growth $G$ with respect to Wo at $ {Re}=2000$ for the physiological waveform. Lines and symbols correspond to the same cases as in (b).

Table 3.

Fourier coefficients ($a_i$, $b_i$) used to approximate the physiological waveform of Bürk et al. (Reference Bürk, Blanke, Stankovic, Barker, Russe, Geiger, Frydrychowicz, Langer and Markl2012).

The fitting of expression (4.1) gives a good estimate of the growth of perturbations for the physiological waveform, even though all data used to obtain the fit coefficients are for $t_{m}\geq 0.3$ (table 1). Larger errors are found for $ {Wo}\leq 10$ and correspond to our definition of $\sigma$ in (4.3). Future analyses may improve this equation to better capture the behaviour at low Womersley numbers.

Figures 9(b) and 9(c) confirm some of the ideas proposed in this paper. First, the energy growth of the helical perturbations is related with the instantaneous instability of the laminar profile due to the presence of inflection points. This is confirmed by the similarity of $G_{TGA}$ with $G_{LSA}$. Second, the shape of the pulsation waveform has an important effect on the energy growth. This is shown in figure 6(a) and evident from the non-negligible values of $w_{m}$, $w_{ac}$ and $w_{dc}$ in our fits.