# Optimal energy growth in pulsatile channel and pipe flows

Benoît Pier, Peter J. Schmid

► To cite this version:

Benoît Pier, Peter J. Schmid. Optimal energy growth in pulsatile channel and pipe flows. Journal of Fluid Mechanics, 2021, 926, pp.A11. 10.1017/jfm.2021.702. hal-03084366v5

HAL Id: hal-03084366

https://hal.science/hal-03084366v5

Submitted on 7 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal energy growth in pulsatile channel and pipe flows

Benoît Pier¹† and Peter J. Schmid²

¹Laboratoire de mécanique des fluides et d'acoustique, CNRS, École centrale de Lyon, Université de Lyon 1, INSA Lyon, 36 avenue Guy-de-Collongue, 69134 Écully, France

²Department of Mathematics, Imperial College, South Kensington Campus, London SW7 2AZ, UK

(Received 25 March 2021; revised 30 July 2021; accepted 30 July 2021)

Pulsatile channel and pipe flows constitute a fundamental flow configuration with significant bearing on many applications in the engineering and medical sciences. Rotating machinery, hydraulic pumps or cardiovascular systems are dominated by time-periodic flows, and their stability characteristics play an important role in their efficient and proper operation. While previous work has mainly concentrated on the modal, harmonic response to an oscillatory or pulsatile base flow, this study employs a direct-adjoint optimisation technique to assess short-term instabilities, identify transient energy-amplification mechanisms and determine their prevalence within a wide parameter space. At low pulsation amplitudes, the transient dynamics is found to be similar to that resulting from the equivalent steady parabolic flow profile, and the oscillating flow component appears to have only a weak effect. After a critical pulsation amplitude is surpassed, linear transient growth is shown to increase exponentially with the pulsation amplitude and to occur mainly during the slow part of the pulsation cycle. In this latter regime, a detailed analysis of the energy transfer mechanisms demonstrates that the huge linear transient growth factors are the result of an optimal combination of Orr mechanism and intracyclic normal-mode growth during half a pulsation cycle. Two-dimensional sinuous perturbations are favoured in channel flow, while pipe flow is dominated by helical perturbations. An extensive parameter study is presented that tracks these flow features across variations in the pulsation amplitude, Reynolds and Womersley numbers, perturbation wavenumbers and imposed time horizon.

**Key words:** shear-flow instability

† Email address for correspondence: benoit.pier@ec-lyon.fr.
# 1. Introduction

Pulsatile flows are a common phenomenon in a variety of engineering flows, and they are ubiquitous in physiological configurations. The pulsatile flow through tubular geometries plays a key role in the haemodynamic system of many species as it is responsible for the transport of oxygenated blood to the organs and muscular tissue (Ku 1997; Pedley 2000). While in many of these configurations inertial effects are too weak to cause and sustain turbulent fluid motion, a variety of cardiovascular diseases can be linked to flow instabilities in the arteries (Chiu & Chien 2011). In addition, geometric modifications of the standard fluid-carrying vessels, such as stenoses, aneurysms or other pathologies, further amplify adverse flow effects and aggravate physiological consequences. For these reasons, a better understanding of pulsatile flows, and the perturbation dynamics they support, would be beneficial, if not mandatory, for improved diagnostics as well as the design of advanced medical devices.

Despite their importance in medical and engineering applications, pulsatile flows — and in particular their stability characteristics — have received far less attention than their steady analogues. Pulsatile flows comprise a steady as well as a time-periodic component. This is in contrast to oscillatory flows which consist of a harmonic part, but lack a steady background flow. The periodic time dependence precludes a standard modal approach, based on temporal Fourier normal modes, and instead calls for more complex methods, such as Floquet analysis. Furthermore, pulsatile flows are governed by a far larger suite of parameters than steady flows: besides the common Reynolds number $Re$ and the wavenumbers of the perturbations, pulsatile flows depend on the pulsation amplitudes and the non-dimensional frequency (the Womersley number $Wo$). For a non-modal analysis, the time horizon over which growth or decay is measured and the phase shift of the perturbation within a base-flow cycle have to be accounted for as well. Within this high-dimensional parameter space, a rich and varied perturbation dynamics can be observed, with important transitions between distinct flow behaviours.

The stability of pulsatile flow has been addressed by a few key studies that laid the foundation for our current understanding of its perturbation dynamics. An account of the pertinent body of literature has been presented in Pier & Schmid (2017) with emphasis on the modal treatment via Floquet analysis. A resume of earlier work on general time-periodic flows has been presented in Davis (1976). Further notable work by von Kerczek (1982) has built on this foundation and established a framework for the analysis of flows with a harmonic base flow. Generic configurations such as a Stokes layer (Blennerhassett & Bassom 2002) or channel and pipe flow with time-periodic pressure gradients (Thomas et al. 2011), have been investigated with modal techniques and have been mapped out as to their stability characteristics across a range of governing parameters. The influence of wall modifications, such as stenoses or aneurysms, on the overall stability behaviour has been addressed via numerical simulations (see, e.g., Blackburn et al. (2008); Gopalakrishnan et al. (2014)).

The role of pulsation in the transition from laminar to turbulent pipe flow has been recently investigated by Xu et al. (2017) and Xu & Avila (2018). These studies in particular concentrated on the emergence and life cycle of localised ‘puffs’, together with their role in triggering transition in the presence of a pulsating flow component, since the occurrence of turbulent bursts in each cycle has been found to be sensitive on flow parameters and configuration details. A strong influence of the Womersley number has been reported, and a distinct regime-switching across three proposed parameter regions has been observed (Xu et al. 2017). These experimental findings have been further corroborated by direct numerical simulations initiated by a localised perturbation (Xu & Avila 2018). The earlier numerical study by Tuzi & Blondeaux (2008) concluded that at moderate but subcritical
Reynolds numbers only parts of the harmonic cycle (around flow reversal) support turbulent flow via an instability and an associated break of the flow's symmetry.

While the early body of literature on time-periodic flows has concentrated on a modal (Floquet) approach, more recent studies have employed an initial-value perspective on the analysis of perturbation dynamics and energy growth. Biau (2016) has analysed the generic oscillatory Stokes layer as to its potential to support transiently growing perturbations over a forcing cycle. This study isolated the Orr mechanism as the dominant process by which energy amplification could be achieved efficiently for sufficiently high unsteady amplitudes. In particular the decelerating part of the forcing cycle has been identified as prone to strong non-modal growth. Complementary nonlinear simulations further verified that triggering by these mechanisms can yield subcritical transition to turbulence. A similar technique has been applied in a recent study by Xu et al. (2021) for oscillatory and pulsating pipe flow. Among others, they have reported that pulsating pipe flows are generally dominated by helical perturbations. In accordance with Biau (2016), a strong Orr-type mechanism has been found to dominate, once a threshold pulsation amplitude has been exceeded. Again, only half of the forcing cycle supported growth of the kinetic perturbation energy; disturbances have been observed to rapidly reach energy levels that facilitate a transition to turbulent fluid motion, often via localised disturbances.

These latter studies advocate the treatment of pulsatile flow as a generally time-dependent flow, distinct from a periodic Floquet ansatz. Over the past decades, the application of these non-modal techniques to hydrodynamic stability calculations has resulted in a more complete understanding of shear-driven instability phenomena. The generality of this approach (Schmid 2007) is well-suited for assessing pulsatile flow over a range of time scales, thus mapping out the optimal perturbation dynamics over partial and multiple pulsation cycles. This non-modal approach for time-dependent flows is based on a variational principle arising from a partial-differential-equation-constrained optimisation problem. It results in a direct-adjoint system of equations (Luchini & Bottaro 2014) that produce the maximum energy growth of perturbations over a prescribed time horizon. Time-dependent base flows are treated naturally within this formalism, and short-term energy amplification mechanisms, for example over a partial pulsation cycle, can be detected and extracted effectively. Over the past years, this computational framework has been successfully brought to bear on a variety of complex flow configurations (see, e.g., Magri (2019) and Qadri et al. (2021) for applications in reactive flows), and has furnished quantitative stability measures beyond the time-asymptotic limit and without the need for simplifying assumptions.

This article follows up on and extends earlier work (Pier & Schmid 2017) that demonstrated the influence of a pulsating flow component on the stability of channel flow via a linear (Floquet) and nonlinear analysis. In this present study, we focus on non-modal effects and the occurrence of transient energy amplification mechanisms under conditions that are asymptotically stable, both for rectangular channel and cylindrical pipe flows. The unsteady nature of the base flow lends itself to a formulation as a partial-differential-equation-constrained optimisation problem for the maximum energy gain which is subsequently solved by a variational approach based on direct-adjoint looping.

The main finding, and significance, of our investigation consists of the quantification of extremely large transient growth, brought on by the unsteady nature of the base flow. By considering both channel and pipe flows and carefully studying energy transfer mechanisms, we identify the fundamental mechanisms responsible for this huge growth, common to both geometries. This amplification potential translates directly into a strong sensitivity for the rise of coherent structures over one or many pulsatile cycles. While this feature of pulsatile flows has been observed and reported in previous studies, an encompassing treatment of this phenomenon, including its presence in parameter space and its manifestation in dominant
spatial structures, is still missing in the literature on unsteady flows. Our findings also have a direct connection to classifying transition scenarios in wall-bounded flows under the influence of cyclic base flow variations, thus extending the classical scenarios for steady flows and potential routes for the transition to turbulence occurring during part of the cycle.

Despite our attempt to analyse pulsating channel and pipe flows comprehensively, judicious choices had to be made to arrive at an emerging picture for the perturbation dynamics prevailing in these configurations. The ensuing parameter ranges have been selected to capture the most compelling and representative flow phenomena, while limiting our focus to flows encountered in physiological and medical situations. Haemodynamic applications, across a range of blood vessel geometries, are well covered by our choice of parameters. Nonetheless, configurations outside this parameter range are touched upon as well, to establish continuity or bifurcations in flow behaviour and to connect to other studies that investigate such parameter regimes in more detail, e.g., Xu et al. (2021).

The present paper represents the culmination of several years of work; a preliminary version of the main results has been presented at the 12th European Fluid Mechanics Conference in Vienna (Pier & Schmid 2018).

## 2. Flow configurations and governing equations

This investigation considers viscous incompressible flow through infinite channels and pipes of constant diameter. In this context, a flow is characterised by a velocity vector field $u(x, t)$ and a scalar pressure field $p(x, t)$ that depend on position $x$ and time $t$ and are governed by the Navier–Stokes equations

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u = \nu \Delta u - \nabla p, \qquad (2.1) $$

$$ 0 = \nabla \cdot u, \qquad (2.2) $$

where $\nu$ is the kinematic viscosity of the fluid, and the pressure has been redefined to eliminate the constant fluid density.

The channel-flow configuration calls for a formulation using Cartesian coordinates, while cylindrical coordinates are appropriate for pipe flows. In order to address both configurations with similar mathematical and numerical tools, we adopt a general formalism using three spatial coordinates $x_0, x_1, x_2$ and associated velocity components $u_0, u_1, u_2$. When analysing channel flow with respect to a Cartesian reference frame, the variables $x_0, x_1$ and $x_2$ denote wall-normal, streamwise and spanwise coordinates, respectively, while they stand for radial, streamwise and azimuthal coordinates when studying pipe flow in a cylindrical setting. Whatever the configuration, the flow domain corresponds to $|x_0| < D/2$ where $D$ is the channel or the pipe diameter, and no-slip boundary conditions prevail along the solid walls at $|x_0| = D/2$.

A formulation of the incompressible Navier–Stokes equations (2.1,2.2) in cylindrical coordinates comprises more terms than one in Cartesian coordinates. Nevertheless, the resulting equations have a very similar structure, and the above notations allow us to cast the governing equations into a single general system of partial differential equations, pertaining to both channel and pipe configurations, the details of which are given in Appendix A.

## 3. Base flows and non-dimensional control parameters

Pulsatile base flows driven by a spatially uniform and temporally periodic streamwise pressure gradient are obtained as exact solutions of the Navier–Stokes equations and consist of a velocity field in the streamwise $x_1$-direction with profiles that only depend on time $t$ and
on the wall-normal/radial coordinate $x_0$. Denoting by $\Omega$ the pulsation frequency, the base velocity profiles may be expanded as temporal Fourier series,

$$U_1(x_0, t) = \sum_n U_1^{(n)}(x_0) \exp(in\Omega t), \qquad (3.1)$$

and are associated with a periodic flow rate

$$Q(t) = \sum_n Q^{(n)} \exp(in\Omega t). \qquad (3.2)$$

In the above expressions, the conditions $Q^{(-n)} = [Q^{(n)}]^*$ and $U_1^{(-n)}(x_0) = [U_1^{(n)}(x_0)]^*$ ensure that all flow quantities are real (with $\star$ denoting a complex conjugate).

By invariance of these base flows in the streamwise $x_1$-direction, the different harmonics in the expansion (3.1) are not coupled through the nonlinear terms of the Navier-Stokes equations and the velocity components $U_1^{(n)}(x_0)$ are analytically obtained by solving simple differential equations derived for each harmonic component. The mean-flow component $U_1^{(0)}(x_0)$ displays a parabolic Poiseuille profile. For $n \neq 0$, following Womersley (1955), the profiles $U_1^{(n)}(x_0)$ are obtained in terms of Bessel functions in cylindrical coordinates corresponding to pipe flows, while they are obtained in terms of exponential functions in Cartesian coordinates corresponding to channel flows.

Pulsatile channel or pipe flows are characterised by the Womersley number

$$Wo \equiv \frac{D}{2} \sqrt{\frac{\Omega}{\nu}}, \qquad (3.3)$$

which is a non-dimensional measure of the pulsation frequency, and may be interpreted as the ratio of the pipe radius (or the channel half-diameter) to the thickness $\delta = \sqrt{\nu/\Omega}$ of the oscillating boundary layers developing near the walls. A pulsatile base flow is then completely specified by the Fourier components $Q^{(n)}$ of its flow rate (3.2), and the velocity profiles of the different harmonics (3.1) are obtained as

$$U_1^{(n)}(x_0) = \frac{Q^{(n)}}{A} W \left( \frac{x_0}{D/2}, \sqrt{n} Wo \right). \qquad (3.4)$$

In the above expression, $A$ denotes the relevant measure of the cross-section ($A = D$ for channels and $A = \pi D^2/4$ for pipes) and the function $W$ is the normalised velocity profile pertaining to each harmonic component. The analytic expressions of $W$ for channel and pipe flows are given in Appendix B.

In this investigation, we only consider pulsatile flows with a non-vanishing mean flow rate $Q^{(0)}$. Thus, the definition of the Reynolds number may be based on mean velocity $Q^{(0)}/A$, diameter $D$ and viscosity $\nu$, leading to

$$Re \equiv \frac{Q^{(0)}}{\nu} \text{ for channels and } Re \equiv \frac{Q^{(0)}}{\nu} \frac{4}{\pi D} \text{ for pipes.} \qquad (3.5)$$

Moreover, using $Q^{(0)}$ as reference, the flow rate waveform is completely determined by the non-dimensional ratios

$$\tilde{Q}^{(n)} \equiv \frac{Q^{(n)}}{Q^{(0)}}, \qquad (3.6)$$

corresponding to the amplitude (and phase) of the oscillating flow rate components ($n > 0$) relative to the mean flow.

In order to reduce the dimensionality of the control-parameter space for the rest of this
paper, we will only consider base flow rates with a single oscillating component

$$Q(t) = Q^{(0)}(1 + \tilde{Q} \cos \Omega t), \qquad (3.7)$$

where the pulsation amplitude $\tilde{Q} \equiv 2\tilde{Q}^{(1)}$ may be assumed real without loss of generality. Note that the theoretical and numerical methods developed for the present investigation are also suitable for studying the dynamics of pulsating base flows with higher harmonic content.

## 4. Mathematical formulation

This entire study considers the dynamics of small-amplitude perturbations developing in the basic pulsatile channel and pipe flows specified in the previous section. The incompressible Navier-Stokes equations are, therefore, linearised about these base flows. Considering that the base flows do not depend on the streamwise coordinate $x_1$ nor on the spanwise/azimuthal coordinate $x_2$, infinitesimally small velocity and pressure perturbations may thus be written as spatial normal modes of the form

$$\mathbf{u}^d(x_0, t) \exp i(\alpha_1 x_1 + \alpha_2 x_2), \qquad (4.1)$$

$$p^d(x_0, t) \exp i(\alpha_1 x_1 + \alpha_2 x_2), \qquad (4.2)$$

where $\alpha_1$ and $\alpha_2$ are streamwise and spanwise/azimuthal wavenumbers, respectively. Separation of total flow fields into basic and perturbation quantities and substitution of the expansions (4.1,4.2) into the governing equations (2.1,2.2) linearised about the relevant time-periodic base flow then yields a system of coupled linear partial differential governing equations of the form

$$A \partial_t \mathbf{q}(x_0, t) = \mathbf{L}(x_0, t) \mathbf{q}(x_0, t), \qquad (4.3)$$

where

$$\mathbf{q}(x_0, t) = \begin{pmatrix} u_0^d(x_0, t) \\ u_1^d(x_0, t) \\ u_2^d(x_0, t) \\ p^d(x_0, t) \end{pmatrix} \quad \text{and} \quad \mathbf{A} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (4.4)$$

Here, the superscript $d$ refers to components of the direct problem, to be distinguished from the adjoint variables below (4.6). The spatial differential operator $\mathbf{L}(x_0, t)$ in (4.3) is a 4-by-4 matrix and its coefficients involve $\partial_0$-differentiation, depend on the wavenumbers $\alpha_1$ and $\alpha_2$ as well as on the base velocity profiles $U_1(x_0, t)$; see Appendix C for explicit expressions of all these terms.

When studying transient growth effects and searching for optimal initial perturbations that are maximally amplified over a finite-time horizon, it is necessary to choose an appropriate measure of disturbance size (Schmid 2007). Using a classical energy-based inner product, the adjoint governing equations associated with the direct problem (4.3) are routinely obtained as

$$A \partial_t \mathbf{q}^{\dagger}(x_0, t) = \mathbf{L}^{\dagger}(x_0, t) \mathbf{q}^{\dagger}(x_0, t), \qquad (4.5)$$

and the adjoint differential operator $\mathbf{L}^{\dagger}(x_0, t)$ is explicitly given in Appendix D. In contrast with the direct equations (4.3), the adjoint equations (4.5) have negative diffusion coefficients and the adjoint fields

$$\mathbf{q}^{\dagger}(x_0, t) = \begin{pmatrix} u_0^a(x_0, t) \\ u_1^a(x_0, t) \\ u_2^a(x_0, t) \\ p^a(x_0, t) \end{pmatrix} \qquad (4.6)$$
are integrated backwards in time.

Denoting by $\{q(x_0, t_i), |x_0| < D/2\}$ an initial perturbation at time $t_i$, the evolution of this disturbance at subsequent times $t > t_i$ and the associated perturbation energy $E(t)$ are then obtained by solving the initial-value problem corresponding to (4.3) with $q(x_0, t_i)$ specified for $|x_0| < D/2$. The temporal evolution of the perturbation amplitude is then characterised by the ratio $E(t)/E(t_i)$, for $t > t_i$.

The maximum possible amplification of a disturbance over the interval $t_i < t < t_f$ is obtained as

$$G(t_i, t_f) = \max_{\{q(x_0, t_i)\}} \frac{E(t_f)}{E(t_i)}, \quad (4.7)$$

by optimising over all possible initial conditions at $t = t_i$. Note that, since the base flow is time-periodic, the amplification factor depends not only on the duration $t_f - t_i$ of the temporal evolution but also on the phase of its starting point $t_i$ within the pulsation cycle.

The particular initial condition at $t = t_i$ that achieves the largest amplification at $t = t_f$ is referred to as the optimal perturbation and the resulting flow fields at $t = t_f$ as the optimal response. In practice, the amplification factors $G(t_i, t_f)$ and associated optimal perturbations and responses are iteratively computed by successive direct–adjoint loops, consisting of temporal integration of the direct equations (4.3) from $t_i$ to $t_f$ and of the adjoint equations (4.5) from $t_f$ to $t_i$, using the numerical methods described in the next section.

In linearly stable configurations, all perturbations eventually decay and the maximal transient growth for given wavenumbers $\alpha_1$ and $\alpha_2$,

$$G^{\max}(\alpha_1, \alpha_2) = \max_{t_i, t_f} G(t_i, t_f; \alpha_1, \alpha_2), \quad (4.8)$$

is well defined and takes finite values. Obviously, $G^{\max}(\alpha_1, \alpha_2)$ also depends on the base flow configuration and its control parameters. For a given pulsating base flow, the largest possible transient amplification that may be achieved is obtained as

$$G_{\max}^{\max} = \max_{\alpha_1, \alpha_2} G^{\max}(\alpha_1, \alpha_2), \quad (4.9)$$

by considering all possible wavenumbers.

## 5. Numerical implementation

The direct and adjoint temporal evolution problems (4.3,4.5) are first order in time and involve spatial differential operators in the wall-normal $x_0$-coordinate.

For spatial discretisation we use a Chebyshev spectral method with collocation points spanning the whole diameter of the channel or the pipe. Whether considering channel or pipe flows, all computations are restricted to half of the domain, $0 \le x_0 \le D/2$, by taking into account the symmetry or antisymmetry of the different flow fields and using the associated discretised differential operators of corresponding symmetry. For channel flow calculations carried out in Cartesian coordinates, the parity of the different flow fields depends on the sinuous or varicose nature of the perturbation under consideration. Note that for all the channel flow configurations considered in this paper, the dynamics is dominated by sinuous perturbations. For pipe flow calculations carried out in cylindrical coordinates, it is the value (even or odd) of the azimuthal mode number that determines the parity of each of the different flow fields. Note that the singularities in the differential operators at the pipe axis ($x_0 = 0$) are only “apparent” (Boyd 2001): the exact solution is analytic at the axis even though the coefficients of the differential equations are not. Thus a consistent implementation
of the symmetry/antisymmetry conditions at the axis removes any apparent singularities and guarantees that the spectral method yields smooth solutions.

Time-marching of the direct and adjoint incompressible Navier–Stokes equations uses a second-order accurate predictor-corrector fractional-step method, derived from Raspo et al. (2002). In classical fashion, the maximal gain $G(t_i, t_f)$, together with optimal initial perturbation and final response, is then obtained by direct-adjoint loops, maximising the energy growth from $t = t_i$ to $t = t_f$. All subsequent quantities $G^{\max}$ and $G_{\max}^{\max}$ are derived from the gain $G$, by maximising over $t_i$ and $t_f$, and over $\alpha_1$ and $\alpha_2$.

Resorting to the general formulation of the governing equations detailed in Appendix A and taking advantage of the relevant symmetry properties of the different flow fields thus leads to a numerical implementation capable of handling all configurations of the present investigation.

This entire numerical solution procedure is a generalisation of an approach already used in our previous investigation (Pier & Schmid 2017), and its implementation in C++ is based on the ‘home-spun’ PackstaB library (Pier 2015, §A.6). The interested reader is referred to these references for further details of the general method.

## 6. Pulsating channel flow

The objective of the present section, which is the core part of the paper, is to investigate how the well known transient-growth properties of steady channel flow are modified by the presence of a pulsating base flow component. Starting with a steady Poiseuille flow, the approach consists of studying the influence of pulsation as the amplitude $\tilde{Q}$ is increased from 0 for different values of the Womersley number Wo.

First, we consider the growth rates $G$ of streamwise-invariant and spanwise-periodic streaks since they display the largest transient growth for Poiseuille flow. Then, the strikingly different behaviour observed for two-dimensional (spanwise-invariant) flows calls for a systematic computation of all possible three-dimensional perturbations. Having established the optimal amplification rates $G^{\max}$ that prevail over the whole wavenumber plane, we are then in a position to derive the maximal achievable energy amplification $G_{\max}^{\max}$ for a given pulsating base flow and to document its dependence on the pulsation amplitude $\tilde{Q}$, the Womersley number Wo and the Reynolds number Re. Finally, a detailed discussion of the energy transfer mechanisms allows us to highlight the various growth mechanisms that come into play during the different stages of the evolution and to explain the huge growth factors that are observed for pulsating flows, already for moderate pulsation amplitudes. We recall that sinuous perturbations prevail for all the situations investigated here; thus, all the results presented in this section correspond to flow fields of sinuous symmetry.

The vast parameter space of the problem requires a systematic exploration of the flow physics and a concentration on essential characteristics by a progressive compression of the governing parameters. To this end, we successively investigate the growth of streaks, two-dimensional and three-dimensional disturbances, before focusing on transient energy growth and the energy transfer mechanisms that accompany the observed amplifications. We conclude by isolating the shape and dynamics of the two- and three-dimensional structures that optimally exploit the unsteady background flow and thus exhibit maximal energy growth. Along this analysis, we present a sequential reduction of the parameter space, starting from the effect of cycle length and cycle phase, via spatial scales to the time horizon for optimal growth. Within each step, the essential features of the transient instability will be presented, before reducing the parameter dependency for the subsequent analysis. This section then culminates in the detailed examination of the most amplified disturbances, for the two- and three-dimensional case, under the influence of a pulsatile background flow.
<figure><img src="image_2.png" /><figcaption>Figure 1: Optimal energy growth for streaks with $\alpha_2 = 4$ and $\alpha_1 = 0$ in steady channel Poiseuille flow at $Re = 1000, 2000, \dots, 5000$. (a) Duration $t_f - t_i$ of growth phase measured in mean-flow advection time scale $\tau_Q$. (b) Rescaled growth factors $G/\text{Re}^2$ and $t_f - t_i$ measured in diffusion time scale $\tau_v = Re \tau_Q$.</figcaption></figure>

## 6.1. Growth of streaks

In steady channel Poiseuille flow, largest transient growth is known to occur for initial conditions which are spanwise periodic and consist of streamwise aligned vortices, thus corresponding to perturbations with $\alpha_1 = 0$ and $\alpha_2 \neq 0$. Figure 1(a) shows the optimal transient amplification at $Re = 1000, 2000, \dots, 5000$ computed for $\alpha_2 = 4$, which is near the most transiently amplified spanwise wavenumber. (Throughout this paper, length scales are non-dimensionalised with respect to the channel (or pipe) diameter $D$.) For a steady base flow, the energy growth factor $G(t_i, t_f)$ only depends on the duration $t_f - t_i$, here measured in mean-flow advection units $\tau_Q \equiv D^2/Q^{(0)}$. Replotting these data for $G/\text{Re}^2$ and measuring the duration $t_f - t_i$ in diffusion units $\tau_v \equiv D^2/\nu = \text{Re} \tau_Q$, the curves in figure 1(b) confirm the known scaling laws, leading to a maximum transient growth of $G^{\text{max}} \approx 1.1 \ 10^{-4} \text{Re}^2$ at $t^{\text{max}}/\tau_Q \approx 1.9 \ 10^{-2} \text{Re}$.

Adding to this steady base flow a pulsatile component of given amplitude and frequency, the transient growth properties are characterised by $G(t_i, t_f)$ which then depends both on the phase of the initial perturbation $t_i$ within the pulsation period $T \equiv 2\pi/\Omega$ and on the duration $t_f - t_i$ of the temporal evolution. For $Re = 2000$ and $Re = 5000$, figure 2 shows plots of the growth factors $G(t_i, t_f)$ for pulsation amplitudes $\tilde{Q} = 0.4$ and $1.0$ at $Wo = 10$. It is found that the amplitude $\tilde{Q}$ of the base flow modulation only weakly influences the streak growth. Even increasing $\tilde{Q}$ to values larger than unity (corresponding to negative flow rates during part of the pulsation cycle), does not significantly alter the distribution of $G(t_i, t_f)$: the maximum amplification remains at the same level and the growth hardly depends on the phase $t_i/T$. Thus streamwise-invariant ($\alpha_1 = 0$) perturbations appear to be almost unaffected by the time-dependent component of the base flow and to display a dynamics predominantly dictated by the time-averaged base flow. The discussion of energy transfer mechanisms in section §6.5 below will shed further light on this observation. Comparing figure 2(a) with 2(c), and 2(b) with 2(d), the similarity observed between plots at different $Re$ and same $\tilde{Q}$ also indicates that the scaling of $G$ with $Re^2$ remains valid for the transient growth of streaks in pulsating base flows.
<figure><img src="image_3.png" /><figcaption>Figure 2: Optimal transient amplification for streaks with $\alpha_2 = 4$ and $\alpha_1 = 0$ at (a,b) Re = 2000 and (c,d) Re = 5000. Pulsating channel flow at $\bar{W}_o = 10$ and (a,c) $\tilde{Q} = 0.4$ and (b,d) $\tilde{Q} = 1.0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Optimal energy growth for two-dimensional perturbations $\alpha_1 = 2$ and $\alpha_2 = 0$ in channel Poiseuille flow at Re = 1000, 2000, ..., 5000.</figcaption></figure>

## 6.2. Growth of two-dimensional perturbations

Two-dimensional spanwise invariant perturbations, corresponding to $\alpha_2 = 0$ and $\alpha_1 \neq 0$, exhibit much weaker transient amplification than streaks for the same steady Poiseuille flow. Figure 3 plots the transient growth properties prevailing for Poiseuille flow at Re = 1000, 2000, ..., 5000 for perturbations with $\alpha_1 = 2$ and $\alpha_2 = 0$, near the most unstable two-dimensional perturbation. Here, the maximal amplification $G^{\text{max}}$ scales linearly with the Reynolds number and reaches much lower values than those corresponding to streaks (see figure 1a); note that this maximal amplification is also reached for a much shorter time horizon.

The evolution of two-dimensional transient growth properties as the amplitude $\tilde{Q}$ of the pulsating component is increased is given in figure 4. After increasing $\tilde{Q}$, a second maximum emerges in the plot of $G$, located around $t_i/T = 0.2$ and $(t_f - t_i)/T = 0.5$. In contrast with
<figure><img src="image_5.png" /><figcaption>Figure 4: Optimal transient amplification for two-dimensional perturbations with $\alpha_2 = 0$ and $\alpha_1 = 2$ at (a-c) Re = 2000 and (d-f) Re = 5000. Pulsating base flow at $Wo = 10$ and (a,d) $\tilde{Q} = 0.04$, (b,e) $\tilde{Q} = 0.10$ and (c,f) $\tilde{Q} = 0.20$.</figcaption></figure>

the situation prevailing for streaks, this second maximum is seen to rapidly grow with $\tilde{Q}$ and to become the dominant feature, here already for $\tilde{Q} \approx 0.1$. While streaks display much larger transient growth for steady Poiseuille flow, these two-dimensional perturbations are found to become the most efficient optimal perturbations for pulsatile base flows, beyond some threshold value of the pulsating amplitude $\tilde{Q}$. This overwhelming growth of two-dimensional perturbations for pulsatile conditions will be explained in section §6.5, below, by detailed monitoring of the amplification process in comparison with the dynamics of temporal Floquet eigenmodes.

### 6.3. Growth of three-dimensional perturbations

The very different transient growth behaviour observed for streaky and two-dimensional perturbations calls for a systematic investigation in the entire $(\alpha_1, \alpha_2)$-wavenumber plane. For a given pulsating base flow, the optimal energy amplification $G^{\max}(\alpha_1, \alpha_2)$ (4.8) is obtained by maximising the transient growth $G(t_i, t_f; \alpha_1, \alpha_2)$ over all values of $t_i$ and $t_f$ at each prescribed wavenumber. We have systematically explored the control-parameter space spanning the ranges $1000 \le Re \le 5000$, $5 \le Wo \le 15$ and $0 \le \tilde{Q} \le 1$, and a few characteristic results are presented below.

The plot of $G^{\max}$ for steady Poiseuille flow ($\tilde{Q} = 0$) at Re = 4000 (figure 5a) confirms that strongest transient growth occurs for streaks (with $\alpha_1 = 0$) and that the largest value of $G^{\max} \approx 1763$ is reached at $\alpha_2 \approx 4.09$ (indicated by a black dot). Two-dimensional
<figure><img src="image_6.png" /><figcaption>Figure 5: Isolines of maximum energy growth $G^{\text{max}}$ in $(\alpha_1, \alpha_2)$-wavenumber plane for base flows with $\text{Re} = 4000$ (a–g) and $\text{Re} = 2000$ (h–j). (a,h) Steady Poiseuille flow ($\tilde{Q} = 0$), (b,d,f,i) pulsating base flows with $\tilde{Q} = 0.2$, (c,e,g,j) pulsating base flows with $\tilde{Q} = 0.5$. Womersley numbers: $\text{Wo} = 8$ in panels (b,c), $\text{Wo} = 10$ in panels (d,e,i,j), $\text{Wo} = 12$ in panels (f,g). The black dots indicate the wavenumbers where $G^{\text{max}}$ reaches its largest value.</figcaption></figure>

perturbations (with $\alpha_2 = 0$) experience growth factors that are two orders of magnitude smaller, with $G^{\text{max}} \approx 30$ for $\alpha_1 = 3.1$.

The distribution of maximal amplification factors $G^{\text{max}}$ in the $(\alpha_1, \alpha_2)$-plane evolves significantly as the amplitude $\tilde{Q}$ of the pulsating component is increased for a given pulsation frequency. Figures 5(b–g) reveal that, as $\tilde{Q}$ is increased, the maximum energy growth (indicated by a black dot) rapidly switches over from streaks to two-dimensional perturbations that experience growth factors sharply increasing with $\tilde{Q}$ while those experienced by streaks (along the $\alpha_2$-axis) do not much depend on $\tilde{Q}$ nor on $\text{Wo}$. Comparison of the results obtained with $\text{Wo} = 8$ (plots b,c), $\text{Wo} = 10$ (d,e) and $\text{Wo} = 12$ (f,g) demonstrates that the rate of
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of maximal transient energy amplification $G_{\max}^{\text{max}}$ with $\tilde{Q}$ for $5 \le W_o \le 15$ at (a) $R_e = 1000$, (b) $R_e = 2000$, (c) $R_e = 4000$ and (d) $R_e = 5000$.</figcaption></figure>

increase of $G_{\max}^{\text{max}}$ with $\tilde{Q}$ varies significantly with $W_o$ and is larger for lower values of the Womersley number.

Figures 5(h–j) illustrate the behaviour at $R_e = 2000$. For steady Poiseuille flow (figure 5h), the isolines of $G_{\max}^{\text{max}}$ display a similar structure as for $R_e = 4000$ (figure 5a) but with lower levels. After increasing the amplitude $\tilde{Q}$ of the pulsating flow component at $W_o = 10$, figures 5(i,j) show that two-dimensional perturbations again eventually dominate the response. However, at this lower Reynolds number, a larger value of $\tilde{Q}$ is required for the two-dimensional perturbations to emerge, and the increase of $G_{\max}^{\text{max}}$ with $\tilde{Q}$ also occurs at a lower rate. Thus $G_{\max}^{\text{max}}$ is found to reach values of the order of $10^5$ at $R_e = 2000$ for $\tilde{Q} = 0.5$ and $W_o = 10$ (figure 5j), while at $R_e = 4000$ values in excess of $10^{11}$ are observed (figure 5e).

### 6.4. Maximal transient growth

The maximal transient energy amplification achievable for a given base flow has been defined as $G_{\max}^{\text{max}}$ (4.9) and is derived by maximising $G^{\max}(\alpha_1, \alpha_2)$ over the entire wavenumber plane.

Figure 6 plots the evolution of $G_{\max}^{\text{max}}$ as the pulsation amplitude $\tilde{Q}$ is continuously increased for Womersley and Reynolds numbers in the range $5 \le W_o \le 15$ and $1000 \le R_e \le 5000$, respectively. At low values of $\tilde{Q}$, the pulsating flow component has a very weak influence and $G_{\max}^{\text{max}}$ remains near the value prevailing for steady Poiseuille flow at the same Reynolds number. For these low pulsation amplitudes, the optimal initial perturbation corresponds to streaks (with $\alpha_1 = 0$) and the associated growth duration $t_f - t_i$ remains very close to that prevailing for the equivalent mean Poiseuille flow (see also figure 8, below).

Beyond some critical value of $\tilde{Q}$, the amplification factor $G_{\max}^{\text{max}}$ starts to increase exponentially with $\tilde{Q}$, as illustrated by the nearly constant slopes in the logarithmic plots of figure 6 (note the different vertical scale used in plot (a) for $R_e = 1000$). This critical value $\tilde{Q}_c$ of the pulsation amplitude depends on $W_o$ and $R_e$ as shown in figure 7(a): increasing the Reynolds number is found to promote the two-dimensional perturbations which become the dominant feature already for $\tilde{Q} > 0.1$ around $R_e = 5000$. For $\tilde{Q} > \tilde{Q}_c$, the rate of the exponential
<figure><img src="image_8.png" /><figcaption>Figure 7: (a) Critical values $\tilde{Q}_c$ for transition between streaky and two-dimensional maximally amplified perturbations. (b) Exponential growth rate $\kappa$ of $G_{\max}^{\max}$ with $\tilde{Q}$ in the two-dimensional regime.</figcaption></figure>

growth of $G_{\max}^{\max}$ with $\tilde{Q}$ corresponds to the slopes seen in figure 6 and significantly increases as the Womersley number decreases. As a result, $G_{\max}^{\max}$ rapidly reaches ‘astronomical’ values, several orders of magnitude beyond the amplification rates prevailing for the corresponding steady Poiseuille flows. These exponential rates have been computed as

$$ \kappa \equiv \frac{1}{G_{\max}^{\max}} \frac{\partial G_{\max}^{\max}}{\partial \tilde{Q}}, \qquad (6.1) $$

and their variation with Re and Wo is given in figure 7(b). In this plot the values of $\kappa$ have been computed by taking the average over $\tilde{Q}_c < \tilde{Q} < \tilde{Q}_c + 0.1$, but note that the growth rate remains nearly constant over much larger intervals of $\tilde{Q}$ in the two-dimensional regime. Obviously, the growth rates are enhanced with the Reynolds number and they also significantly increase towards the lower Womersley numbers, corresponding to longer pulsation periods.

The regime change in the transient growth behaviour occurring for $\tilde{Q} = \tilde{Q}_c$ is further illustrated in figure 8 at Re = 4000. The evolution with pulsation amplitude $\tilde{Q}$ of the streamwise $\alpha_1$ and spanwise $\alpha_2$ wavenumbers associated with the maximally amplified optimal perturbations display a sharp transition from streaky ($\alpha_1 = 0, \alpha_2 \neq 0$) to two-dimensional ($\alpha_1 \neq 0, \alpha_2 = 0$) perturbations. For Wo = 6 and 8, the spanwise wavenumber here directly switches from $\alpha_2 \approx 4$ to 0. For higher values of Wo, however, a small range in $\tilde{Q}$ is observed where the maximally amplified perturbations consist of oblique waves with small but finite values of $\alpha_2$. This corresponds to configurations where the amplification of two-dimensional perturbations is still in competition with streaks, so that the maximum of $G^{\max}$ in the ($\alpha_1, \alpha_2$)-plane occurs slightly off the $\alpha_1$-axis, as illustrated by the black dot in figure 5(f) and corresponding dots in figures 8(a,b). It is then only for higher values of $\tilde{Q}$ that purely two-dimensional ($\alpha_2 = 0$) perturbations prevail.

The transition from streaky to two-dimensional maximally amplified perturbations is also accompanied by a significant change in the duration of the growth phase $t_f - t_i$, shown in figure 8(c) in mean-flow advection units $\tau_Q$ and in figure 8(d) in units of the pulsation period $T$. These two time scales are associated with different dynamical features and related as $Wo^2 T = \pi^2 / 2 Re \tau_Q$. At weak pulsation amplitudes $\tilde{Q}$, the duration $t_f - t_i$ remains very close to the value prevailing for streaks developing in the equivalent steady Poiseuille flow, here approximately $75 \tau_Q$ at Re = 4000 (compare with figure 1). At higher pulsation amplitudes, when two-dimensional perturbations dominate, maximal amplification occurs over intervals $t_f - t_i$ that approximately correspond to half a pulsation period, $T/2$. Thus, the transition
<figure><img src="image_9.png" /><figcaption>Figure 8: Characterisation of the maximally amplified optimal perturbations as the pulsation amplitude $\tilde{Q}$ is increased for $W_o = 6, 8, 10, 12, 14$ at $Re = 4000$. (a) Streamwise wavenumber $\alpha_1$, (b) spanwise wavenumber $\alpha_2$, (c) duration of transient growth $t_f - t_i$ measured in mean-flow advection units $\tau_Q$ and (d) in pulsation periods $T$.</figcaption></figure>

from streaky to two-dimensional perturbations also coincides with a change in the dynamical time scale: from streak growth essentially dictated by the mean flow to two-dimensional perturbations strongly amplified over half a pulsation cycle.

### 6.5. Discussion of energy transfer mechanisms

In this final subsection on channel flows, we investigate the energy production and dissipation mechanisms in order to explain the different transient-growth scenarios that have been identified.

Following the notations introduced in section §4, we consider a perturbation of the form

$$ \mathbf{u}(x_0, t) \exp i(\alpha_1 x_1 + \alpha_2 x_2) + \text{c.c.} \quad \text{with} \quad \mathbf{u}(x_0, t) = \begin{pmatrix} u_0(x_0, t) \\ u_1(x_0, t) \\ u_2(x_0, t) \end{pmatrix}, \qquad (6.2) $$

using a complex-valued three-dimensional velocity vector $\mathbf{u}(x_0, t)$. Such a perturbation is associated with an instantaneous kinetic energy per unit volume of

$$ E(t) = \frac{1}{D} \int_{-D/2}^{+D/2} e(x_0, t)dx_0, \qquad (6.3) $$

where

$$ e(x_0, t) = \|\mathbf{u}(x_0, t)\|^2 \equiv \mathbf{u}(x_0, t) \cdot [\mathbf{u}(x_0, t)]^* \qquad (6.4) $$
represents the local energy density. Thus, the temporal energy variation,

$$ \frac{dE(t)}{dt} = \frac{1}{D} \int_{-D/2}^{+D/2} (\partial_t \mathbf{u}(x_0, t) \cdot [\mathbf{u}(x_0, t)]^* + \mathbf{u}(x_0, t) \cdot [\partial_t \mathbf{u}(x_0, t)]^*) dx_0, \quad (6.5) $$

follows from the dynamics of $\mathbf{u}(x_0, t)$, governed by the Navier-Stokes equations linearised about the pulsating base flow (4.3). Separating terms due to interaction with the base flow from those involving viscous dissipation leads to

$$ \frac{dE(t)}{dt} = \Pi(t) - \Theta(t), \quad (6.6) $$

where

$$ \Pi(t) = \frac{1}{D} \int_{-D/2}^{+D/2} \pi(x_0, t) dx_0 \quad \text{and} \quad \Theta(t) = \frac{1}{D} \int_{-D/2}^{+D/2} \theta(x_0, t) dx_0, \quad (6.7) $$

with

$$ \pi(x_0, t) = - \frac{\partial U_1(x_0, t)}{\partial x_0} \left( u_0(x_0, t)[u_1(x_0, t)]^* + [u_0(x_0, t)]^* u_1(x_0, t) \right) \quad (6.8) $$

and

$$ \theta(x_0, t) = 2\nu \left( \|\partial_0 \mathbf{u}(x_0, t)\|^2 + (\alpha_1^2 + \alpha_2^2) \|\mathbf{u}(x_0, t)\|^2 \right). \quad (6.9) $$

The term $\pi(x_0, t)$ accounts for energy transfer between the pulsating base flow and the perturbation: it essentially represents energy production due to base-flow shear, but negative values may occur and its profile across the channel crucially depends on the relative phases of $u_0(x_0, t)$ and $u_1(x_0, t)$.

Another quantity of interest is the instantaneous growth rate

$$ \sigma(t) = \frac{1}{2E(t)} \frac{dE(t)}{dt} = \frac{\Pi(t) - \Theta(t)}{2E(t)}, \quad (6.10) $$

particularly relevant during phases of near-exponential amplification.

Close monitoring of the spatiotemporal development of the base-flow interaction $\pi(x_0, t)$ and the dissipation $\theta(x_0, t)$ terms will clarify the amplification mechanisms that govern the different stages of the dynamics.

We focus on two characteristic configurations that have already been discussed: pulsating base flows at Re = 4000 and Wo = 10 with two different pulsation amplitudes, $\tilde{Q} = 0.1$ and $\tilde{Q} = 0.2$, associated with streaky and two-dimensional maximally amplified perturbations, respectively.

### 6.5.1. Streaky maximally amplified optimal perturbation

For the lower pulsation amplitude of $\tilde{Q} = 0.1$, a maximal amplification of $G_{\max}^{\max} = 1.77 \ 10^3$ is achieved from $t_i = 0.166T$ to $t_f = 1.389T$ for streamwise invariant and spanwise periodic perturbations with $\alpha_1 = 0$ and $\alpha_2 = 4.073$. The associated temporal evolution of the perturbation energy $E(t)$ is shown in figure 9(a), with the corresponding instantaneous growth rate $\sigma(t)$ in figure 9(b). Here, the transient growth is seen to follow the classical pattern prevailing for steady Poiseuille flow: a strong and very short initial boost for $t_i < t < t_\star = 0.175T$ (blue parts of the curves), followed by a phase of gradually weakening growth for $t_\star < t < t_f$ (in red) towards the maximum response. And indeed, these curves in figure 9(a,b) are almost identical to the accompanying insets that correspond to the maximally amplified perturbations for steady Poiseuille flow at the same Reynolds number, characterised by $\alpha_1 = 0$, $\alpha_2 = 4.088$, $G_{\max}^{\max} = 1.76 \ 10^3$. This evolution is the result of energy production $\Pi(t)$ and dissipation $\Theta(t)$, shown in figure 9(c). As can be seen by plotting these
<figure><img src="image_10.png" /><figcaption>Figure 9: Temporal development of maximally amplified perturbation at $Re = 4000$, $Wo = 10$ and $\tilde{Q} = 0.1$. Optimal perturbation is streamwise invariant with $\alpha_1 = 0$ and $\alpha_2 = 4.073$. (a) Evolution of energy $E(t)$ from $t_i = 0.166T$ to $t_f = 1.389T$, leading to $G_{\max}^{\max} = 1.7710^3$. (b) Corresponding instantaneous growth rate $\sigma(t)$. (c) Energy production $\Pi(t)$ (solid line) and dissipation $\Theta(t)$ (dashed line). (d) Production and dissipation terms relative to instantaneous energy. Insets in panels (a) and (b) correspond to maximally amplified streaks for steady Poiseuille flow at same Reynolds number.</figcaption></figure>

quantities relative to the instantaneous energy in figure 9(d), viscous dissipation plays here a minor part in the transient growth throughout the entire process from $t_i$ to $t_f$.

The temporal evolution of the spatial structure of the maximally amplified streaky perturbation is illustrated in figure 10. Selected snapshots correspond to the thick black dots in figure 9: $t_i = 0.166T$ optimal initial perturbation (thick blue curves); $t = 0.169T$ (thin blue curves); $t_\star = 0.175T$ at maximal instantaneous growth (thick green curves); $t = 0.500T$ (thin red curves); $t_f = 1.389T$ optimal response (thick red curves). In order to enable comparison of these profiles throughout the temporal evolution, they have here all been normalised to unit total energy. As expected, the initial perturbation consists in streamwise aligned vortices, that fill the entire channel cross-section, with a vanishing streamwise velocity component: see thick blue curves in figures 10(a-c) and corresponding vector plot in figure 10(e). Transient amplification promotes streamwise velocity while reducing wall-normal and spanwise velocity components, leading to a final response that solely consists of streamwise velocity: see thick red curves in figures 10(a-c) and $u_1$-isolines in figure 10(f). The energy production profiles $\pi$ shown in figure 10(d) result from the interaction of base flow shear with $u_0$ and $u_1$, and are therefore significant only around $t_\star = 0.175T$ (green curve), while displaying vanishing levels near $t_i$ and $t_f$. Dissipation profiles $\theta$ (not shown) remain at small values throughout the entire evolution.

Clearly, in this regime, the oscillating component of the base flow has a very weak influence, the amplification process operates as for the equivalent steady Poiseuille flow by
<figure><img src="image_11.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f)</figcaption></figure>

Figure 10: Evolution of spatial structure of maximally amplified streaks for $Re = 4000$, $Wo = 10$ and $\tilde{Q} = 0.1$. Spatial profiles of flow fields (normalised to unit total energy) over half-channel $0 \le x_0 \le D/2$ at different snapshots: $t_i = 0.166T$ optimal initial perturbation (thick blue lines); $t = 0.169T$ (thin blue lines); $t_\star = 0.175T$ at maximum growth rate (thick green lines); $t = 0.500T$ (thin red lines); $t_f = 1.389T$ optimal response (thick red lines). Envelope of (a) wall-normal $|u_0(x_0, t)|$, (b) streamwise $|u_1(x_0, t)|$ and (c) spanwise $|u_2(x_0, t)|$ velocity perturbations. (d) Energy production $\pi(x_0, t)$ terms. Snapshots of velocity fields in half-channel over two spanwise wavelengths ($\lambda_2 = 2\pi/\alpha_2$): (e) vector plot of $(u_0, u_2)$ for initial perturbation at $t_i$ and (f) equispaced isolines of streamwise component $u_1$ of response at $t = t_f$.

redistributing streamwise momentum by streamwise vortices, and the dynamics is essentially dictated by the lift-up effect. This insensitivity to the pulsating base-flow component explains why the maximal growth factors $G_{\text{max}}^{\text{max}}$ prevailing for streaky optimal perturbations remain at nearly constant level as $\tilde{Q}$ is increased, as observed in figure 6.
<figure><img src="image_17.png" /><figcaption>Figure 11: Temporal development of maximally amplified perturbation at Re = 4000, Wo = 10 and $\tilde{Q} = 0.2$. Optimal perturbation is two-dimensional with $\alpha_1 = 2.619$ and $\alpha_2 = 0$. (a) Evolution of energy $E(t)$ from $t_i = 0.168T$ to $t_f = 0.698T$, leading to</figcaption></figure>

$G_{\max}^{\max} = 5.48 \ 10^4$. (b) Corresponding instantaneous growth rate $\sigma(t)$. (c) Energy production $\Pi(t)$ (solid line) and dissipation $\Theta(t)$ (dashed line). (d) Production and dissipation terms relative to instantaneous energy. Growth occurs in two stages: phase I (in blue) for $t_i < t < t_\star$ and phase II (in red) for $t_\star < t < t_f$, separated by stall at $t_\star = 0.251T$ (green vertical line).

### 6.5.2. Two-dimensional maximally amplified optimal perturbation

A markedly different scenario prevails at higher pulsation amplitudes when the largest transient amplification is achieved for two-dimensional (streamwise periodic and spanwise invariant) perturbations.

As already shown in figure 5(d), for a pulsation amplitude of $\tilde{Q} = 0.2$, the maximally amplified optimal initial perturbation at Re = 4000 and Wo = 10 occurs for $\alpha_1 = 2.619$ and $\alpha_2 = 0$ and leads to an amplification of $G_{\max}^{\max} = 5.48 \ 10^4$ from $t_i = 0.168T$ to $t_f = 0.698T$. The temporal evolution of the perturbation energy is given in figure 11(a), with the associated instantaneous growth rate in figure 11(b). The transient growth that occurs over the interval $t_i < t < t_f$ is here seen to develop in two stages: first a relatively short period (phase I, blue curves) of rapid growth followed by a longer interval (phase II, red curves) of weaker but almost constant growth. Between these two stages, the amplification stalls and the instantaneous growth displays a minimum value, which is here slightly negative around $t_\star = 0.251T$. This two-stage evolution results from production and dissipation contributions, as illustrated in figure 11(c,d): a first peak in $\Pi(t)$ during phase I is responsible for the rapid growth of the perturbation, followed by a sustained nearly exponential increase during phase II. The contribution of the relative dissipation $\Theta(t)/2E(t)$ is significant only at the very beginning, before rapidly dropping to low values.

The mechanisms responsible for the growth of the perturbation differ in both phases,
<figure><img src="image_18.png" /><figcaption>Figure 12: Evolution of spatial structure of maximally amplified (two-dimensional) optimal perturbation for Re = 4000, Wo = 10 and $\bar{Q} = 0.2$. Spatial profiles of flow fields (normalised to unit total energy) over half-channel $0 \le x_0 \le D/2$ at different snapshots: $t_i = 0.168T$ optimal initial perturbation (thick blue lines); $t = 0.186T$ and $t = 0.210T$ in phase I (thin blue lines); $t_\star = 0.251T$ at stall (thick green lines); $t = 0.300T$ and $t = 0.500T$ in phase II (thin red lines); $t_f = 0.698T$ optimal response (thick red lines).</figcaption></figure>

(a) Envelope of wall-normal $|u_0(x_0, t)|$ and (b) streamwise $|u_1(x_0, t)|$ velocity perturbations. (c) Energy production $\pi(x_0, t)$ and (d) dissipation $\theta(x_0, t)$ terms. Snapshots of $(u_0, u_1)$ velocity fields in half-channel over two streamwise wavelengths ($\lambda_1 = 2\pi/\alpha_1$):

(e) initial perturbation at $t = t_i$ and (f) response at $t = t_f$.

as illustrated by the profiles in figure 12. These plots show the evolution of the spatial distribution of various fields by selected snapshots, corresponding to the thick black dots in figure 11. Note that the perturbation profiles have again been normalised to unit total energy. The envelopes of wall-normal $|u_0|$ and streamwise $|u_1|$ velocity components in figure 12(a,b) show that the initial perturbation at $t_i$ (thick blue curves) is localised toward
<figure><img src="image_19.png" /><figcaption>Figure 13: Temporal energy evolution of maximally amplified perturbation (blue and red curves) compared with least stable normal mode (thick black curve), at $Re = 4000$, $Wo = 10$ and $\dot{Q} = 0.2$, for $\alpha_1 = 2.619$. Normal-mode energy exponentially decays in the long term, according to a Floquet multiplier of $\mu_F = 0.0766$, but displays intracyclic amplification by a factor of $G_{nm} = 4.30 \ 10^3$, approximately during half a pulsation period. Optimal perturbation is amplified by $G_{max}^{max} = 5.48 \ 10^4$. Phase II (in red) closely follows normal mode while phase I (in blue) is very similar to optimal growth prevailing for steady Poiseuille flow at same parameters (inset). Grey sinusoidal curve represents base flow rate $Q(t)$ (not to scale).</figcaption></figure>

the wall with maximum amplitude around $x_0 = 0.4D$, before spreading over the entire channel cross-section in the subsequent evolution. The spatial distribution of the base-flow interaction terms $\pi(x_0, t)$ (shown in figure 12c) reveals that the driving mechanism is strong and concentrated around $x_0 = 0.4D$ in phase I (blue curves) while weaker and evenly spread out in phase II (red curves). In contrast, plots of $\theta(x_0, t)$ (figure 12d) show that dissipation is only significant in the initial stages for $t \approx t_i$ and reduced to a very thin boundary layer near the wall throughout the rest of the evolution. The vector plot of the initial perturbation in a streamwise channel cross-section (figure 12e) highlights the flow structures concentrated near the wall and characteristically tilted upstream. In contrast, the final response (figure 12f) fills the entire channel.

Finally, we compare the dynamics of the present maximally amplified optimal perturbation with the development of the temporal normal mode prevailing for the same pulsating base flow at the same spatial wavenumbers. Such normal modes have been extensively computed and characterised in our previous investigation (Pier & Schmid 2017), using both Floquet eigenproblems and linearised direct numerical simulations.

All pulsatile base flows under consideration here are linearly stable so that temporal eigenmodes decay in the long term. The thick black curve in figure 13 shows the temporal evolution of perturbation energy for the least stable normal mode at $Re = 4000$, $Wo = 10$ and $\dot{Q} = 0.2$, with $\alpha_1 = 2.619$ and $\alpha_2 = 0$. The negative mean slope in this logarithmic plot confirms the decay, governed by a Floquet multiplier of $\mu_F = 0.0766$. Thus, the perturbation energy of this normal mode is reduced by a factor of $\mu_F^2$ after each pulsation period. However, within each pulsation cycle, significant modulation occurs. This intracyclic growth and decay has been shown to approximately coincide with base-flow deceleration and acceleration phases (Pier & Schmid 2017), as indicated by the grey sinusoidal line representing $Q(t)$. Here, the normal mode displays an intracyclic amplification of $G_{nm} = 4.30 \ 10^3$.

Comparison of optimal-perturbation and normal-mode energy curves reveals that, during phase II ($t_\star < t < t_f$, red part of curve in figure 13), the optimal perturbation closely follows
the normal-mode dynamics. And, indeed, optimal perturbation and normal mode also display very similar flow fields during that interval.

In the initial phase I ($t_i < t < t_\star$, blue part of curve in figure 13), however, the maximally amplified perturbation takes advantage of the optimal initial condition responsible for the initial boost in the response through the Orr mechanism. The amplification during phase I is almost identical to the maximal growth experienced by a two-dimensional optimal initial condition for steady Poiseuille flow at the same Reynolds number and same streamwise wavenumber, shown in the inset in figure 13.

These considerations reveal that the maximally amplified two-dimensional perturbation is an optimal combination of Orr mechanism (phase I) and intracyclic normal-mode growth over half a pulsation cycle (phase II). Growth during phase I is essentially determined by the equivalent steady Poiseuille base flow: the resulting amplification factor therefore scales approximately linearly on Re while being largely independent of Wo and $\tilde{Q}$. In contrast, growth during phase II closely follows the intracyclic amplification of the associated temporal eigenmodes, and the magnitude of this intracyclic growth has been shown to strongly depend on Wo and $\tilde{Q}$: whatever the Womersley number, it increases almost exponentially with $\tilde{Q}$, and the increase is fastest at the lower values of Wo (Pier & Schmid 2017). This exponential growth with $\tilde{Q}$ explains why two-dimensional optimal perturbations always eventually prevail over streaky perturbations, as observed in figure 6.

The maximal growth factor $G_{\max}$ is obviously always larger than either contribution of phase I or of phase II to the total growth. But while the contribution of phase I remains at moderate levels (one or two orders of magnitude, as for steady Poiseuille flow), it is phase II that is responsible for the huge amplification factors prevailing as the modulation amplitude $\tilde{Q}$ increases. As a result, except for weak pulsation amplitudes, the Orr mechanism only contributes a small factor to the maximal growth $G_{\max}$, while most of the amplification process is due to modal growth during base-flow deceleration.

# 7. Pulsating pipe flow

After having presented detailed results for channel flows, we now turn to the transient growth properties of pulsating flows through circular pipes. The organisation of this section is similar to the previous one. However, since most features are equivalent, many details may be omitted here.

By adopting the general formulation appropriate for both Cartesian and cylindrical coordinates, the analysis of pulsating pipe flows is carried out with the same numerical codes as previously used for pulsating channel flows. Due to periodicity in the azimuthal coordinate, the wavenumber $\alpha_2$ only takes integer values, but otherwise the numerical implementation proceeds as for a Cartesian formulation. Recall that the apparent singularity at the pipe axis ($x_0 = 0$) resolves itself by taking advantage of the symmetry properties relevant for each flow component, since all flow fields are either symmetric or antisymmetric in the radial coordinate $x_0$.

## 7.1. *Transient growth of streaks and helical perturbations*

Since for steady Hagen-Poiseuille flow, streamwise invariant streaks with $\alpha_2 = 1$ and $\alpha_1 = 0$ undergo the largest non-modal growth, we first consider the transient amplification features prevailing for the same type of perturbations developing in pulsating pipe flows. Figure 14 shows the amplification factors $G(t_i, t_f)$ obtained at $Wo = 10$, for $Re = 2000$ and $5000$, $\tilde{Q} = 0.4$ and $1.0$. The control parameters are the same as those used in figure 2 for pulsating channel flow, and it is observed that the transient growth properties are very similar.
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Figure 14: Optimal transient amplification for streaks with $\alpha_2 = 1$ and $\alpha_1 = 0$ at (a,b) Re = 2000 and (c,d) Re = 5000. Pulsating pipe flow at $Wo = 10$ and (a,c) $\tilde{Q} = 0.4$ and (b,d) $\tilde{Q} = 1.0$.</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>Figure 15: Optimal transient amplification for helical perturbations with $\alpha_2 = 1$ and $\alpha_1 = 2$ at (a–c) Re = 2000 and (d–f) Re = 5000. Pulsating pipe flow at $Wo = 10$ and (a,d) $\tilde{Q} = 0.1$, (b,e) $\tilde{Q} = 0.4$ and (c,f) $\tilde{Q} = 0.6$.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 16: Maximum energy growth $G^{\text{max}}$ for pulsating pipe flows, over $0 < \alpha_1 < 6$ at $\tilde{Q} = 0.0, 0.1, 0.2, \dots, 1.0$. Here (a–f) Re = 4000 and (g–i) Re = 2000; (a–c,h) Wo = 10, (d,g) Wo = 8, (e) Wo = 12 and (f,i) Wo = 14; (a) $\alpha_2 = 0$, (b,d–i) $\alpha_2 = 1$ and (c) $\alpha_2 = 2$.</figcaption></figure>

For streamwise periodic ($\alpha_1 \neq 0$) perturbations, the least stable temporal modes correspond to helical perturbations with $\alpha_2 = 1$. Investigation of transient growth characteristics for $\alpha_1 \neq 0$ also confirms that perturbations with $\alpha_2 = 1$ dominate over axisymmetric perturbations ($\alpha_2 = 0$) as well as over those of higher azimuthal order ($\alpha_2 \ge 2$).

Figure 15 illustrates the transient growth properties for $\alpha_2 = 1$ and $\alpha_1 = 2$ at $Wo = 10$ and $Re = 2000$ and $Re = 5000$ as the amplitude $\tilde{Q}$ of the pulsating base flow component is increased. As for pulsating channel flow, a second maximum emerges that rapidly dominates the dynamics beyond some value of the base flow modulation amplitude. This maximum is again located near $t_i/T = 0.2$ and $(t_f - t_i)/T = 0.5$ and corresponds thus to amplification over half a pulsation cycle.
## 7.2. Optimal growth at given wavenumbers

The maximal transient growth $G^{\max}$, computed by optimisation of $G(t_i, t_f)$ over all values of $t_i$ and $t_f$ for fixed wavenumbers $\alpha_1$ and $\alpha_2$, is shown in figure 16. In each plot the evolution of $G^{\max}$ curves for $0 < \alpha_1 < 6$ is illustrated as $\tilde{Q}$ is increased from $\tilde{Q} = 0$ to $\tilde{Q} = 1$ in steps of 0.1. Panels (a–c) compare the growth of axisymmetric perturbations, $\alpha_2 = 0$ in panel (a), with that of helical perturbations, $\alpha_2 = 1$ in panel (b) and $\alpha_2 = 2$ in panel (c). Clearly, under pulsating flow conditions, axisymmetric initial conditions undergo transient amplification that is not much larger than for the equivalent steady Poiseuille flow, as demonstrated by the nearly overlapping curves in figure 16(a). Non-axisymmetric perturbations, however, experience transient amplification that rapidly grows with $\tilde{Q}$, and strongest growth occurs for $\alpha_2 = 1$ (figure 16b). Computation of $G^{\max}$ for all $\alpha_2 \le 6$ (results not shown) reveals that the same scenario prevails at higher azimuthal order, but the rate of increase of $G^{\max}$ with $\tilde{Q}$ is significantly weaker for higher $\alpha_2$.

Evolution of the growth characteristics for Re = 4000 with different Womersley numbers, Wo = 8 in panel (d), Wo = 12 in panel (e) and Wo = 14 in panel (f), confirms again that largest amplification factors occur for lower pulsation frequencies, i.e., longer pulsation cycles.

Finally, values obtained at lower Re = 2000 for Wo = 8 in panel (g), 10 in panel (h) and 14 in panel (i) show that the general trend is similar but with lower values of $G^{\max}$, as expected for lower Re.

## 7.3. Maximal amplification

Finally, the maximal amplification $G_{\max}^{\text{max}}$ achievable for a given pulsating pipe flow is obtained by optimising $G^{\max}(\alpha_1, \alpha_2)$ over all streamwise wavenumbers $\alpha_1$ and azimuthal mode numbers $\alpha_2$. Figure 17 shows the variation of $G_{\max}^{\text{max}}$ as the pulsation amplitude $\tilde{Q}$ is increased for Womersley numbers in the range $5 \le \text{Wo} \le 15$ and Re = 2000, 3000, 4000 and 5000. The behaviour is again found to be similar to that prevailing for pulsating channel flows: at low pulsation amplitudes, $G_{\max}^{\text{max}}$ hardly departs from the value corresponding to the equivalent steady Poiseuille flow; beyond a critical value $\tilde{Q}_c$ of the pulsation amplitude $\tilde{Q}$, transition to approximately exponential growth of $G_{\max}^{\text{max}}$ with $\tilde{Q}$ takes over. The results shown in figure 17(a) perfectly match those of figure 4(a) of Xu et al. (2021), for the subset of control parameter values that is common to both studies. This agreement further validates our methods.

The variation with Wo and Re of the critical value $\tilde{Q}_c$ for crossover between the two regimes is shown in figure 18(a). In the exponential regime prevailing for $\tilde{Q} \ge \tilde{Q}_c$, the growth rates $\kappa$ corresponding to the slopes in figure 17 are given in figure 18(b), computed according to equation (6.1). For a given Reynolds number, the curves of $\tilde{Q}_c$ in figure 18(a) are seen to display a minimum for moderate values of the Womersley number, while they increase both for large and small values of Wo. The increase of $\tilde{Q}_c$ with Wo for $\text{Wo} \ge 10$ is strongest at lower values of the Reynolds number. In contrast, for $\text{Wo} \le 10$ the values of $\tilde{Q}_c$ depend much less on Re.

Comparison of the values of $\tilde{Q}_c$ and $\kappa$ for pipe flows (figure 18) with those prevailing for channel flows shown in figure 7, reveals that pipe flows require larger pulsation amplitudes to switch to the regime with exponentially growing amplification factors $G_{\max}^{\text{max}}$. This is especially true for lower Womersley numbers (see also figure 20 below with additional data for Wo = 3). Also, while the growth rates $\kappa$ display very similar trends for both channel (figure 7b) and pipe configurations (figure 18b), the values for pipe flows are approximately half those of channel flows.

The streamwise wavenumber $\alpha_1$ associated with the most amplified perturbation as the
<figure><img src="image_31.png" /><figcaption>Figure 17: Evolution of maximal transient energy amplification $G_{\text{max}}^{\text{max}}$ for pulsating pipe flows with $\tilde{Q}$ for $W_o = 5, 6, 7, \dots, 15$ at (a) $Re = 2000$, (b) $Re = 3000$, (c) $Re = 4000$ and (d) $Re = 5000$.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 18: (a) Critical values $\tilde{Q}_c$ for transition between streaky and helical maximally amplified perturbations in pulsating pipe flows. (b) Exponential growth rate $\kappa$ of $G_{\text{max}}^{\text{max}}$ with $\tilde{Q}$ in the helical regime.</figcaption></figure>

pulsation amplitude $\tilde{Q}$ is varied for a range of Womersley numbers is monitored in figure 19(a) and (b) for $Re = 2000$ and $Re = 4000$, respectively. These plots demonstrate that the regime change occurring at $\tilde{Q}_c$ is indeed associated with a jump in streamwise wavenumber from $\alpha_1 = 0$ for $\tilde{Q} < \tilde{Q}_c$ to finite $\alpha_1$-values for $\tilde{Q} > \tilde{Q}_c$. In contrast with channel flows, however, for all pulsating pipe flow configurations investigated here, the optimal perturbations always occur with azimuthal mode number $\alpha_2 = 1$. Thus the critical value $\tilde{Q}_c$ always corresponds to a transition from streaky ($\alpha_1 = 0, \alpha_2 = 1$) to helical ($\alpha_1 \neq 0, \alpha_2 = 1$) optimal perturbations, at the same azimuthal mode number.

This regime change is also associated with a discontinuity in the duration of the growth phase $t_f - t_i$ for the optimal amplification process, as illustrated in figure 19(c,d) for $Re =$
<figure><img src="image_33.png" /><figcaption>Figure 19: Characterisation of the maximally amplified optimal perturbations as the pulsation amplitude $\tilde{Q}$ is increased for $Wo = 6, 8, 10, 12, 14$. Streamwise wavenumber $\alpha_1$ at (a) Re = 2000 and (b) Re = 4000. Duration of transient growth $t_f - t_i$ at Re = 4000 (c) measured in mean-flow advection units $\tau_Q$ and (d) in pulsation periods $T$.</figcaption></figure>

4000. The optimised duration $t_f - t_i$ is given in mean-flow advection units $\tau_Q$ in figure 19(c) and in units of the pulsation period $T$ in figure 19(d). These plots illustrate that pulsating pipe flows display similar transient dynamics as channel flows: for $\tilde{Q} < \tilde{Q}_c$, the optimal duration $t_f - t_i$ remains close to the value prevailing for the average parabolic flow profile; for $\tilde{Q} > \tilde{Q}_c$, when helical perturbations dominate, maximal amplification occurs over intervals corresponding approximately to half a pulsation period. Thus our results confirm the findings of Xu *et al.* (2021) that helical perturbations dominate the transient growth at large pulsation amplitudes. By our detailed comparison of channel and pipe configurations at moderate pulsation frequencies, we have been able to highlight the fundamental growth mechanisms, which are common to both geometries.

## 8. Conclusion

Considering pulsating flows through both channels and pipes, we have investigated the non-modal transient energy amplification resulting from optimal initial conditions. Our study has systematically covered the pulsating base flows for $1000 \le Re \le 5000, 5 \le Wo \le 15$ and $0 \le \tilde{Q} \le 1$.

While channel and pipe flows display quite different linear modal stability characteristics, their non-modal transient growth features are found to be very similar in situations that are linearly stable. Optimal energy growth occurs according to two distinct scenarios. At weak pulsation amplitudes $\tilde{Q}$, the behaviour is similar to that resulting from the equivalent steady
<figure><img src="image_34.png" /><figcaption>Figure 20: Evolution of maximal transient energy amplification $G_{\text{max}}^{\text{max}}$ with $\tilde{Q}$ for $Wo = 3$, 5, 10, 15 and 20 at $Re = 5000$, for (a) channel and (b) pipe configurations.</figcaption></figure>

Poiseuille flow, and the oscillating flow component appears to have only a small effect. Beyond a critical value of $\tilde{Q}$, however, transient growth increases exponentially with $\tilde{Q}$ and reaches astronomical values, already for moderate pulsation amplitudes. In this latter regime, optimal growth mainly occurs over half a pulsation period, during the slow part of the pulsation cycle, and closely follows the intracyclic amplification of the associated Floquet eigenmodes. We have previously shown (Pier & Schmid 2017) that the intracyclic modulation amplitudes derived from temporal normal modes may be huge, even for linearly decaying eigenmodes. The maximal transient amplification factors $G_{\text{max}}^{\text{max}}$ computed in the present investigation are even larger since they take advantage of both this normal-mode intracyclic growth and non-modal Orr-type amplification, which contributes in the early stage of the growth process.

These findings have been firmly established by a comprehensive investigation of pulsating flows over the range $5 \le Wo \le 15$, deemed to be the most relevant for applications in the haemodynamic context. In order to explore the expected behaviour beyond that frequency range, figure 20 shows the maximal achievable transient growth $G_{\text{max}}^{\text{max}}$ for channel and pipe flows at $Re = 5000$, extending the results of figures 6(d) and 17(d) by including data at lower and higher pulsation frequencies, $Wo = 3$ and $Wo = 20$, respectively. At high pulsation frequencies, it is observed that the pulsating component is rather inefficient in producing $G_{\text{max}}^{\text{max}}$ factors beyond those prevailing for steady base flows, a result closely related to the fact that high-frequency pulsation also has a strong stabilising effect on modal temporal growth rates. In the low frequency regime, the curves for $Wo = 3$ indicate that strong growth is still possible but requires larger pulsation amplitudes $\tilde{Q}$. When lowering $Wo$, the critical value $\tilde{Q}_c$ for onset of the exponential regime increases moderately for channel flows and significantly for pipe flows ($\tilde{Q}_c \approx 0.98$ for $Wo = 3$). These plots are in agreement with observations already made by Xu *et al.* (2021) for pulsating pipe flows. Concerning the spatial structure of the optimal perturbations, the results of figure 20 follow the same scenario as previously discussed: while streaky perturbations prevail for $\tilde{Q} < \tilde{Q}_c$, at larger pulsation amplitudes two-dimensional and helical perturbations dominate, respectively, in channel and pipe flows.

It should be noted that a major difference between channel and pipe flows concerns their linear modal instability features. Indeed, for channel flows there exists a critical Reynolds number beyond which linear instability occurs. This is well known for steady Poiseuille channel flow, and the dependence of this critical Reynolds number with the pulsating flow parameters has been extensively discussed in our previous work (Pier & Schmid 2017). By contrast, steady pipe Poiseuille flow remains linearly stable, whatever the Reynolds number. For time-periodic base flows, linear instability has been found for purely oscillating pipe flows (Thomas *et al.* 2011). However, the presence of a non-vanishing mean flow rate appears
to have a stabilising effect and all pulsating pipe flows considered in the present study are far
from temporal instability.

Another difference is that two-dimensional (spanwise invariant and streamwise periodic) perturbations are the most unstable or the least stable for channel flows, whereas the leading linear instability for pipe flows occurs for helical modes with $\alpha_2 = 1$ and $\alpha_1 \neq 0$, which dominate over perturbations of higher azimuthal order ($\alpha_2 \ge 2$) as well as over axisymmetric ($\alpha_2 = 0$) ones. It is found that this remains true for pulsating pipe flows.

While channel flows are rapidly dominated by two-dimensional sinuous perturbations, pipe flows are dominated by helical perturbations in similar pulsating flow regimes. For pipe flows, axisymmetric perturbations never prevail. But note that the Cartesian equivalent of axisymmetric perturbations are two-dimensional perturbations of varicose symmetry, which never prevail either. The closest equivalent to a two-dimensional sinuous perturbation in a cylindrical geometry is a helical perturbation (with $\alpha_2 = 1$).

This study gives a detailed and comprehensive perspective on the perturbation dynamics in pulsatile channel and pipe flow, treating these configurations within a time-dependent, initial-value problem formalism and thus avoiding restrictive assumptions of a modal, time-asymptotic approach. This analysis identified a rich perturbation behaviour driven by parametric and transient excitation over one or multiple forcing cycles and the dominance of an Orr-type amplification mechanism at early times that acts efficiently and selectively across a significant parameter range, once a critical pulsation amplitude has been surpassed.

Our study lays the foundation for a future analysis of pulsating base flows with higher harmonic content, such as blood flow rates resulting from the cardiac pulse. The present approach could also be generalised to take into account compliant walls or to address nonlinear fluid effects.

## Acknowledgements

Institut national de physique nucléaire et de physique des particules (IN2P3-CNRS) is warmly thanked for very generously providing computing resources. Without access to its Centre de calcul allowing to run thousands of independent single-processor jobs it would have been very difficult to carry out the present research.

The Royal Society is gratefully acknowledged for its financial support through an International Exchanges Cost Share award (IEC/R2/170133).

Early stages of this work were carried out while B.P. was visiting the BP Institute, enjoying a French Government fellowship at Churchill College, Cambridge University.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. General formulation of the Navier–Stokes equations

In order to handle both Cartesian and cylindrical formulations of the governing Navier-Stokes equations (2.1,2.2), a general set of spatial coordinates $x_0, x_1, x_2$ and associated velocity components $u_0, u_1, u_2$ is used. These correspond to either wall-normal, streamwise and spanwise directions for channel flows, or radial, streamwise and azimuthal directions for pipe flows, respectively. Using these coordinates and velocity fields, the components of the
incompressible Navier-Stokes equations become

$$ \partial_t u_0 + (\mathbf{u} \cdot \nabla)u_0 - \left[ \frac{1}{x_0} u_2^2 \right] = \nu \left( \Delta u_0 + \left[ \frac{1}{x_0^2} (-u_0 - 2\partial_2 u_2) \right] \right) - \partial_0 p, \quad (A1) $$

$$ \partial_t u_1 + (\mathbf{u} \cdot \nabla)u_1 = \nu \left( \Delta u_1 \right) - \partial_1 p, \quad (A2) $$

$$ \partial_t u_2 + (\mathbf{u} \cdot \nabla)u_2 + \left[ \frac{1}{x_0} u_0 u_2 \right] = \nu \left( \Delta u_2 + \left[ \frac{1}{x_0^2} (-u_2 + 2\partial_2 u_0) \right] \right) - \left[ \frac{1}{x_0} \right] \partial_2 p, \quad (A3) $$

$$ 0 = \partial_0 u_0 + \left[ \frac{1}{x_0} u_0 \right] + \partial_1 u_1 + \left[ \frac{1}{x_0} \right] \partial_2 u_2, \quad (A4) $$

with the notations

$$ \partial_0 \equiv \partial_{x_0}, \quad \partial_1 \equiv \partial_{x_1}, \quad \partial_2 \equiv \partial_{x_2} \quad (A5) $$

and

$$ \mathbf{u} \cdot \nabla \equiv u_0 \partial_0 + u_1 \partial_1 + \left[ \frac{1}{x_0} \right] u_2 \partial_2, \quad (A6) $$

$$ \Delta \equiv \partial_{00} + \left[ \left( \frac{1}{x_0} \right) \partial_0 \right] + \partial_{11} + \left[ \left( \frac{1}{x_0^2} \right) \partial_{22}. \quad (A7) $$

In these expressions, the terms enclosed in boxes are only present in the formulation using cylindrical coordinates and pertaining to the pipe flow configuration. Resorting to such a general formalism is particularly useful when developing numerical codes to solve both channel and pipe flows: the boxed terms may be switched on or off depending on the flow configuration.

## Appendix B. Analytic expressions of the pulsating base flow profiles

For pulsating base flows prevailing in infinite channels or pipes, the harmonic components $U_1^{(n)}(x_0)$ of the streamwise velocity fields (3.1) display profiles following the shape function $W(\xi, \omega)$ with $\xi = 2x_0/D$ and $\omega = \sqrt{n} \text{Wo}$.

When considering channel flows in Cartesian coordinates, the oscillating velocity profiles are analytically obtained in terms of hyperbolic functions

$$ W(\xi, \omega) \equiv \left( \frac{\cosh(\sqrt{i}\xi\omega)}{\cosh(\sqrt{i}\omega)} - 1 \right) \left/ \left( \frac{\tanh(\sqrt{i}\omega)}{\sqrt{i}\omega} - 1 \right) \right., \quad (B1) $$

for $|\xi| \le 1$ and $\omega \ne 0$, while the steady component is parabolic,

$$ W(\xi, \omega = 0) \equiv \frac{3}{2} (1 - \xi^2). \quad (B2) $$

When considering pipe flows in cylindrical coordinates, the velocity profiles involve Bessel functions

$$ W(\xi, \omega) \equiv \left( \frac{J_0(\sqrt{-i}\xi\omega)}{J_0(\sqrt{-i}\omega)} - 1 \right) \left/ \left( \frac{2}{\sqrt{-i}\omega} \frac{J_1(\sqrt{-i}\omega)}{J_0(\sqrt{-i}\omega)} - 1 \right) \right., \quad (B3) $$
with $J_0$ and $J_1$ denoting the classic Bessel functions of the first kind. The steady component is again parabolic,

$$W(\xi, \omega = 0) \equiv 2(1 - \xi^2). \qquad (B 4)$$

All the profiles above are normalised so that their cross-sectional average equals unity. Thus, the pulsating base flow velocity components (3.4) are simply obtained by multiplying these profiles with the flow rate coefficients $Q^{(n)}$.

## Appendix C. Linear governing equations of direct problem

In the direct formulation of the incompressible Navier–Stokes equations (4.3), the spatial differential operator $L(x_0, t)$ is a 4-by-4 matrix of the form

$$L(x_0, t) \equiv \begin{pmatrix} L_{00}(x_0, t) & 0 & L_{02}(x_0, t) & -\partial_0 \\ L_{10}(x_0, t) & L_{11}(x_0, t) & 0 & -i\alpha_1 \\ L_{20}(x_0, t) & 0 & L_{22}(x_0, t) & -\frac{1}{x_0}i\alpha_2 \\ -\partial_0 - \frac{1}{x_0} & -i\alpha_1 & -\frac{1}{x_0}i\alpha_2 & 0 \end{pmatrix}. \qquad (C 1)$$

Its coefficients involve $\partial_0$-differentiation, depend on the spatial wavenumbers as well as on the base flow velocity profiles. Their explicit expressions are the following:

$$L_{00}(x_0, t) = -i\alpha_1 U_1(x_0, t) + v \left( \Delta - \frac{1}{x_0^2} \right), \qquad (C 2)$$

$$L_{11}(x_0, t) = -i\alpha_1 U_1(x_0, t) + v\Delta, \qquad (C 3)$$

$$L_{22}(x_0, t) = -i\alpha_1 U_1(x_0, t) + v \left( \Delta - \frac{1}{x_0^2} \right), \qquad (C 4)$$

$$L_{10}(x_0, t) = -\partial_0 U_1(x_0, t), \qquad (C 5)$$

$$L_{20}(x_0, t) = -L_{02}(x_0, t) = \boxed{2v \frac{i\alpha_2}{x_0^2}}, \qquad (C 6)$$

with

$$\Delta \equiv \partial_{00} + \left[ \frac{1}{x_0} \frac{\partial_0}{\partial t} - \alpha_1^2 - \frac{1}{x_0^2} \alpha_2^2 \right]. \qquad (C 7)$$

## Appendix D. Adjoint problem

In the adjoint formulation of the incompressible Navier–Stokes equations (4.5), the spatial differential operator $L^\dagger(x_0, t)$ is obtained as

$$L^\dagger(x_0, t) \equiv \begin{pmatrix} L_{00}^\dagger(x_0, t) & L_{01}^\dagger(x_0, t) & L_{02}^\dagger(x_0, t) & -\partial_0 \\ 0 & L_{11}^\dagger(x_0, t) & 0 & -i\alpha_1 \\ L_{20}^\dagger(x_0, t) & 0 & L_{22}^\dagger(x_0, t) & -\frac{1}{x_0}i\alpha_2 \\ -\partial_0 - \frac{1}{x_0} & -i\alpha_1 & -\frac{1}{x_0}i\alpha_2 & 0 \end{pmatrix}, \qquad (D 1)$$
with

$$L_{00}^{\dagger}(x_0, t) = -i\alpha_1 U_1(x_0, t) - v \left( \Delta - \boxed{\frac{1}{x_0^2}} \right), \quad (D 2)$$

$$L_{11}^{\dagger}(x_0, t) = -i\alpha_1 U_1(x_0, t) - v\Delta, \quad (D 3)$$

$$L_{22}^{\dagger}(x_0, t) = -i\alpha_1 U_1(x_0, t) - v \left( \Delta - \boxed{\frac{1}{x_0^2}} \right), \quad (D 4)$$

$$L_{01}^{\dagger}(x_0, t) = \partial_0 U_1(x_0, t), \quad (D 5)$$

$$L_{20}^{\dagger}(x_0, t) = -L_{02}^{\dagger}(x_0, t) = \boxed{2v \frac{i\alpha_2}{x_0^2}}. \quad (D 6)$$

REFERENCES

BIAU, D. 2016 Transient growth of perturbations in Stokes oscillatory flows. *J. Fluid Mech.* **794**, R4.

BLACKBURN, H. M., SHERWIN, S. J. & BARKLEY, D. 2008 Convective instability and transient growth in steady and pulsatile stenotic flows. *J. Fluid Mech.* **607**, 267–277.

BLENNERHASSETT, P. J. & BASSOM, A. P. 2002 The linear stability of flat Stokes layers. *J. Fluid Mech.* **464**, 393–410.

BOYD, J. P. 2001 *Chebyshev and Fourier spectral methods*. New York: Dover.

CHIU, J.-J. & CHIEN, S. 2011 Effects of disturbed flow on vascular endothelium: pathophysiological basis and clinical perspectives. *Physiol. Rev.* **91**, 327–387.

DAVIS, S. H. 1976 The stability of time-periodic flows. *Annu. Rev. Fluid Mech.* **8**, 57–74.

GOPALAKRISHNAN, S. S., PIER, B. & BIESHEUVEL, A. 2014 Dynamics of pulsatile flow through model abdominal aortic aneurysms. *J. Fluid Mech.* **758**, 150–179.

VON KERCZEK, C. H. 1982 The instability of oscillatory plane Poiseuille flow. *J. Fluid Mech.* **116**, 91–114.

KU, D. N. 1997 Blood flow in arteries. *Annu. Rev. Fluid Mech.* **29**, 399–434.

LUCHINI, P. & BOTTARO, A. 2014 Adjoint equations in stability analysis. *Annu. Rev. Fluid Mech.* **46**, 493–517.

MAGRI, L. 2019 Adjoint methods as design tools in thermoacoustics. *Appl. Mech. Rev.* **71**, 020801.

PEDLEY, T. J. 2000 Blood flow in arteries and veins. In *Perspectives in Fluid Dynamics* (ed. G. K. Batchelor, H. K. Moffatt & M. G. Worster), pp. 105–158. Cambridge: Cambridge University Press.

PIER, B. 2015 Dynamique des écoulements ouverts : instabilités et transition, courbure et rotation. Habilitation à diriger des recherches, École centrale de Lyon & Université de Lyon, https://hal.archives-ouvertes.fr/tel-01107517.

PIER, B. & SCHMID, P. J. 2017 Linear and nonlinear dynamics of pulsatile channel flow. *J. Fluid Mech.* **815**, 435–480.

PIER, B. & SCHMID, P. J. 2018 Optimal energy growth in pulsatile channel and pipe flows. 12th European Fluid Mechanics Conference (EFMC12).

QADRI, U., MAGRI, L., IHME, M. & SCHMID, P. J. 2021 Using adjoint-based optimization to enhance ignition in non-premixed jets. *Proc. Roy. Soc. A* **477**, 20200472.

RASPO, I., HUGUES, S., SERRE, E., RANDRIAMAMPIANINA, A. & BONTOUX, P. 2002 A spectral projection method for the simulation of complex three-dimensional rotating flows. *Computers and Fluids* **31**, 745–767.

SCHMID, P. J. 2007 Nonmodal stability theory. *Annu. Rev. Fluid Mech.* **39**, 129–162.

THOMAS, C., BASSOM, A. P., BLENNERHASSETT, P. J. & DAVIES, C. 2011 The linear stability of oscillatory Poiseuille flow in channels and pipes. *Proc. R. Soc. Lond. A* **467**, 2643–2662.

TUZI, R. & BLONDEAUX, P. 2008 Intermittent turbulence in a pulsating pipe flow. *J. Fluid Mech.* **599**, 51–79.

WOMERSLEY, J. R. 1955 Method for the calculation of velocity, rate of flow and viscous drag in arteries when pressure gradient is known. *J. Physiol.* **127**, 553–563.

XU, D. & AVILA, M. 2018 The effect of pulsation frequency on transition in pulsatile pipe flow. *J. Fluid Mech.* **857**, 937–951.

XU, D., SONG, B. & AVILA, M. 2021 Non-modal transient growth of disturbances in pulsatile and oscillatory pipe flows. *J. Fluid Mech.* **907**, R5.
XU, D., WARNECKE, S., SONG, B., MA, X. & HOF, B. 2017 Transition to turbulence in pulsating pipe flow. *J. Fluid Mech.* **831**, 418–432.