Transient linear stability of pulsating Poiseuille flow using optimally time-dependent modes

Abstract

Time-dependent flows are notoriously challenging for classical linear stability analysis. Most progress in understanding the linear stability of these flows has been made for time-periodic flows via Floquet theory focusing on time-asymptotic stability. However, little attention has been given to the transient intracyclic linear stability of periodic flows since no general tools exist for its analysis. In this work, we explore the potential of using the recent framework of the optimally time-dependent (OTD) modes (Babaee & Sapsis, Proc. R. Soc. Lond. A, vol. 472, 2016, 20150779) to extract information about both the transient and the time-asymptotic linear stability of pulsating Poiseuille flow. The analysis of the instantaneous OTD modes in the limit cycle leads to the identification of the dominant instability mechanism of pulsating Poiseuille flow by comparing them with the spectrum and the eigenmodes of the Orr–Sommerfeld operator. In accordance with evidence from recent direct numerical simulations, it is found that structures akin to Tollmien–Schlichting waves are the dominant feature over a large range of pulsation amplitudes and frequencies but that for low pulsation frequencies these modes disappear during the damping phase of the pulsation cycle as the pulsation amplitude is increased beyond a threshold value. The maximum achievable non-normal growth rate during the limit cycle was found to be nearly identical to that in plane Poiseuille flow. The existence of subharmonic perturbation cycles compared with the base flow pulsation is documented for the first time in pulsating Poiseuille flow.

1. Introduction

Since the seminal work of Reynolds, to whom we owe the ubiquitous Reynolds number, on transition to turbulence in pipes (Reynolds 1883) and the discovery of the famed inflection point criterion by Rayleigh (Rayleigh 1879) that are the cornerstones for hydrodynamic stability theory, the field has seen a great deal of development. A central concept in stability analysis is to consider the linearisation of the operator or Jacobian about a base flow and compute its eigenvalues that govern whether infinitesimal perturbations will ultimately grow or decay. This approach has had tremendous success in the analysis of the asymptotic fate of perturbations in steady parallel flows from the early solution of the Rayleigh–Bénard convection (Rayleigh 1916) to the analysis of the viscous linear stability of parallel flow leading to the Orr–Sommerfeld (OS) operator (Orr 1907; Sommerfeld 1908) that has been the object of numerous analytical and numerical studies (Orszag 1971; Drazin & Reid 1981). Despite the progress, it was only relatively recently that the persistent mismatch between the predicted onset of linear instability in many shear flows and the experimental and numerical evidence of transition at much lower Reynolds numbers (Klebanoff, Tidstrom & Sargent 1962; Avila et al. 2011) was explained. In fact, for non-normal operators such as the linearised Navier–Stokes operator, the focus on its spectrum alone blinds the investigation for growth mechanisms exploiting the non-orthogonality of its eigenbasis; instead, the pseudospectra need to be analysed (Reddy & Henningson 1993; Trefethen et al. 1993). This new perspective on the stability of fluid flow puts the focus on transient amplification of initial disturbances to levels at which nonlinear effects take over that bypass the modal path to turbulence (Henningson, Lundbladh & Johansson 1993). This realisation has led to the development of non-modal stability theory that centres the analysis on a linear initial value problem avoiding the restrictions imposed by the normal mode ansatz. For a detailed review on the topic see Schmid (2007). Besides linear stability methods, a number of approaches based on nonlinear theory have been developed. For a recent review on these methods see Kerswell (2018).

In recent years, with the increased availability of computer power, global instability theory has gained interest (Huerre & Monkewitz (1990), see Chomaz (2005) or Theofilis (2011) for reviews of recent developments). This method does not rely on spatial homogeneity of the flow, allowing for much more geometric flexibility at the cost of larger computations that require specialised tools (Theofilis 2003). The combination of the linear stability analysis framework with matrix-free timestepping methods routinely used in flow solvers (Barkley, Blackburn & Sherwin 2008; Bagheri et al. 2009a) has paved the way for global stability and transient growth analyses of increasing complexity and fidelity (Blackburn, Barkley & Shervin 2008; Monokrousos et al. 2010; He et al. 2017; Tsigklifis & Lucey 2017; Chauvat et al. 2020).

A common feature of most of the studied flow cases is that, despite the spatial complexity, they rely on a steady base flow in order to perform the stability analysis. This limits the direct application of the methods to sub-critical Reynolds numbers where flows are naturally stable (Mack & Schmid 2011; He et al. 2017). In unstable scenarios, it requires artificial stabilisation methods to recover a steady base flow such as fixed-point iterations (reviewed in Knoll & Keyes 2004) or relaxation term methods such as selective frequency damping (Åkervik et al. 2006; Loiseau et al. 2014; Chauvat et al. 2020) or, for oscillating flows, the use of time-averaged data (Hammond & Redekopp 1997). In fact, unsteady flows are notoriously challenging to analyse and progress in the understanding of their linear stability characteristics has been slow (Schmid & Henningson 2001). Significant steps towards this goal have been made only for time-periodic flows using energy theory (Davis & Von Kerczek 1973), non-modal stability analysis (Xu, Song & Avila 2021) and especially Floquet theory (Davis 1976). The latter transforms the time-periodic problem into a linear eigenvalue problem by assuming periodic solutions of the same frequency as the base flow that can be solved with standard stability tools for steady states. While this technique has been applied to a wide range of flows (Von Kerczek 1982; Barkley & Henderson 1996; Marques & Lopez 1997; Pier & Schmid 2017), it has the fundamental shortcoming that, although it distinguishes between stable and unstable flows in the time-asymptotic limit, it ignores intracyclic growth and decay that can be substantial enough to trigger turbulence via nonlinear mechanisms.

Recently, the optimally time-dependent (OTD) modes have been proposed as a new framework for the analysis of linear perturbations in time-dependent systems (Babaee & Sapsis 2016). The method generates a time-evolving orthonormal basis of a finite-dimensional subspace around an arbitrary time-dependent trajectory that spans the instantaneously most unstable directions of the tangent space. Following their proposal, a series of works have established a number of characteristics of the OTD modes. In particular, it was shown that they converge exponentially fast to the most unstable eigendirections of the Cauchy–Green tensor, making them well suited for the reduced-order computation of finite-time Lyapunov exponents (Babaee et al. 2017), and that they are identical to continuously orthonormalised Gram–Schmidt vectors for a specific choice of parameters (Blanchard & Sapsis 2019a), which has proved to be a relevant connection for the application of OTD modes. The strength of the OTD basis vectors is their capacity to span the instantaneously most unstable subspace of user defined size in the (possibly infinite-dimensional) phase space and to capture both modal and non-modal growth within that subspace. Furthermore, the OTD basis vectors are flow invariant (Babaee & Sapsis 2016) and are dynamically consistent with the full dynamics (Farazmand & Sapsis 2016; Blanchard, Mowlavi & Sapsis 2018). These characteristics have led to several applications ranging from control of linear instabilities (Blanchard et al. 2018; Blanchard & Sapsis 2019c), prediction of dynamical events in a statistical framework (Farazmand & Sapsis 2016), the computation of sensitivities (Donello, Carpenter & Babaee 2020) to edge tracking (Beneitez et al. 2020), leveraging the ability of the OTD modes to follow the linear dynamics even along a chaotic trajectory.

In these works, the OTD basis was used primarily as a stable and efficient numerical tool to obtain a basis of the most unstable subspace, but the structure of the OTD modes themselves has not been in the focus. In fact, a linearisation around a complex, nonlinear state makes the physical interpretation of the OTD modes difficult. In order to develop a better understanding of the physical relevance of the OTD modes, it is useful to apply the framework to a simpler case in which the base flow is better understood. A good candidate for this purpose is pulsating Poiseuille flow. In this canonical time-periodic parallel flow configuration, the constant pressure gradient of plane Poiseuille flow is modulated harmonically around the mean. Due to the geometric simplicity, the availability of a general analytical solution of the incompressible Navier–Stokes equations for the base flow and the wealth of data in the literature, especially on its Floquet stability, this flow case is ideal as a benchmark for transient linear stability analysis using the OTD framework. In particular, the configuration is not only unsteady but non-autonomous via the forced pressure oscillations making its instantaneous stability properties entirely inaccessible to linear stability analysis via the global mode approach.

The aim of this paper is to give an overview of the theoretical understanding of the OTD framework and to illustrate its potential as a general tool for the linear stability analysis of time-dependent flows, in particular for time-periodic configurations. The general purpose implementation in the spectral element code Nek5000 (Fischer, Lottes & Kerkemeier 2008) is applied to pulsating Poiseuille flow to feature the capabilities and assess the limitations of the method in order to provide a guide for the application of the framework to other flow cases. The instantaneous growth rates and spatial structure of the OTD modes are analysed, the interaction between pulsation and non-normal growth potential is explored and the existence of subharmonic perturbation cycles is documented.

The remainder of this paper is organised as follows. In § 2 we review the OTD framework theory before introducing the formulation for the incompressible Navier–Stokes equations used in this work (§ 3). We then introduce the flow case used to illustrate and validate the implementation in Nek5000 (§ 4) and present the results of the transient linear stability analysis based on the OTD modes (§ 6). A discussion of the main findings and concluding remarks are gathered in § 7.

2. OTD modes

2.1. Problem specification

We consider a general $n$-dimensional dynamical system

(2.1)\begin{equation} \frac{\partial z}{\partial t} = f(z,t),\quad z \in \mathbb{R}^{n},\enspace t \in [t_0, t_0 + T], \end{equation}

where $f : \mathbb {R}^{n} \times [t_0, t_0 + T] \to \mathbb {R}^{n}$ is a nonlinear but sufficiently smooth vector field. We denote by $z(t):=z(t,z_0,t_0)$ a solution of (2.1) with the initial conditions $z(t_0) = z_0$ as the state of a trajectory of the system at time $t$. Infinitesimal perturbations $q(z(t),t)$ about the nonlinear trajectory $z(t)$ are governed by the linearised dynamics

(2.2)\begin{equation} \frac{\partial q}{\partial t} = L(z(t),t) q,\quad q \in \mathbb{R}^{n},\enspace \mbox{with}\ L(z(t),t) := \boldsymbol{\nabla}_z f(z(t),t), \end{equation}

where $L(z(t),t)$ is the Jacobian of the vector field $f$.

The general solution of (2.2) can be expressed using the propagator $\boldsymbol {\nabla } F_{t_0}^{t}$ that maps the perturbations at the initial time $t_0$ to $t$ along the linearised dynamics given by

(2.3)\begin{equation} q(t) = \boldsymbol{\nabla} F_{t_0}^{t} q_0, \end{equation}

where $q_0 = q(t_0)$.

2.2. Preliminaries

When analysing systems of the form (2.2) for their stability properties, the dominant direction of the tangent space is easily recovered e.g. via the impulse response technique. If the goal is to also extract information about subdominant directions, the naive approach to simply integrate a set of vectors using the equations in time is in general doomed to fail. In practice, almost all initial conditions will asymptotically align with the most unstable direction while their magnitude grows or decays exponentially, leading to numerical instability and large errors (Wolf et al. 1985). In order to avoid both the collapse of all perturbations onto the least stable direction as well as numerical over/underflow, the OTD modes were introduced in Babaee & Sapsis (2016) by including an orthonormality constraint into the linear evolution equations. The resulting evolution equations for the OTD basis vectors $q_i$ for $i=1,\ldots ,r$ are given by

(2.4)\begin{equation} \frac{\partial q_i}{\partial t} = Lq_i - \sum_{j=1}^{r} \left( \langle Lq_i, q_j \rangle - \varPhi_{ij} \right) q_j ,\end{equation}

where $\langle \cdot , \cdot \rangle$ is an appropriate inner product, $\varPhi _{ij} \in \mathbb {R}^{r \times r}$ is a skew–symmetric but otherwise arbitrary matrix. We have omitted the explicit time and trajectory dependence of the linearised operator for clarity. The $r$-dimensional subspace spanned by the orthonormal basis vectors $q_i$ is called the OTD subspace and spans the instantaneously most unstable directions in phase space (Babaee & Sapsis 2016). We can subsequently define a reduced operator $L_r \in \mathbb {R}^{r \times r}$ by orthogonal projection of the full operator $L$ onto the OTD subspace

(2.5)\begin{equation} (L_r)_{ij} := \langle Lq_i, q_j \rangle . \end{equation}

The reduced operator thus gives us a computationally tractable approximation of the action of the full Jacobian which is unfeasible in high-dimensional systems (Babaee et al. 2017).

Combining the set of OTD basis vectors in a matrix $\mathbb {R}^{n \times r} \ni Q = [q_1, \ldots , q_r]$, (2.4) can be recast in the following form:

(2.6)\begin{equation} \frac{\partial Q}{\partial t} = LQ - Q(L_r - \varPhi), \end{equation}

where the reduced operator has the form $L_r = Q^{\textrm {T}} L Q$ and $(\cdot )^{\textrm {T}}$ denotes matrix transposition.

From a geometrical point of view, it is instructive to consider (2.4) as the solution to the minimisation of the norm

(2.7)\begin{equation} G \left( \frac{\partial Q}{\partial t} \right) := \sum_{i=1}^{r} \left\| \frac{\partial q_i}{\partial t} - Lq_i \right\|^{2} \end{equation}

under the constraint of orthonormality, i.e. $\langle q_i, q_j \rangle = \delta _{ij}$, where $\delta _{ij}$ is the Kronecker delta. The OTD basis vectors are therefore obtained not by optimising the individual vectors $q_i$ but rather by optimising their rate of change to minimise the Euclidean distance between $\dot{q}_i$ and the image of $q_i$ through the action of the linear operator $L$. It is under this perspective that the name ‘optimally time-dependent’ modes is defined, in sharp contrast to many methods aiming at finding modal decompositions that are themselves optimal in some sense such as proper orthogonal decomposition (POD) or dynamic mode decomposition (DMD) modes. The OTD basis vectors in themselves lack a clear physical interpretation. Instead, they serve as a numerically stable way of obtaining an orthonormal basis of the subspace spanning the most unstable directions in phase space.

The evolution equation for the OTD basis vectors also has a geometrical interpretation, as pointed out in Babaee & Sapsis (2016). The first term that is subtracted from the general linearised dynamics is $QL_r = QQ^{\textrm {T}}\cdot LQ$ which corresponds to an orthogonal projection of the action of the full linear operator onto the OTD subspace. The resulting rate of change is therefore constrained to the orthogonal complement of the OTD subspace. This condition is in fact already sufficient to evolve the basis vectors and corresponds to the trivial choice of $\varPhi$. This is related to the fact that in an $n$-dimensional phase space the orthogonal complement of any basis vector $q_i$ is an $(n-1)$-dimensional hyperplane. The rotation matrix $\varPhi$ therefore describes the permissible movements of the basis vectors relative to each other without violating the orthonormality constraint. More precisely, the off-diagonal terms encode the $n-1$ possible rotations within their respective orthogonal complements. This gives an intuitive motivation for the fact that the space spanned by the OTD basis vectors is unaffected by the choice of $\varPhi$ as long as $\varPhi = -\varPhi ^{\textrm {T}}$ (Babaee & Sapsis 2016). Previous works also sometimes include the rotation matrix in the definition of the reduced operator (Farazmand & Sapsis 2016; Blanchard & Sapsis 2019c).

A further step in the understanding of the OTD framework was taken in Blanchard & Sapsis (2019a), who consider the link to Gram–Schmidt vectors. They show that on choosing the rotation matrix as

(2.8) \begin{equation} \varPhi_{ij} = \begin{cases} -\langle Lq_j,q_i\rangle & \mbox{if } j < i,\\ 0 & \mbox{if } j=i, \\ \langle Lq_j,q_i\rangle & \mbox{if } j > i,\end{cases} \end{equation}

the OTD basis vectors become identical to the continuously orthonormalised Gram–Schmidt vectors and the evolution equation for the $i$th OTD basis vector reduces to

(2.9) \begin{equation} \frac{\partial q_i}{\partial t} = Lq_i - \langle Lq_i,q_i \rangle q_i - \sum_{j=1}^{i-1} (\langle Lq_i,q_j \rangle + \langle Lq_j,q_i \rangle)q_j . \end{equation}

Apart from the theoretical appeal, this choice of $\varPhi$ also has advantageous numerical properties since the evolution equations (2.9) are lower triangular (notice the difference in the limit of the sum compared with (2.4)) and can therefore be efficiently solved through forward substitution (Blanchard & Sapsis 2019a). The original OTD formulation is not hierarchical in the sense that adding a basis vector requires the recomputation of all vectors. The Blanchard & Sapsis (2019a) formulation overcomes this issue.

Another feature of this formulation is that the first OTD basis vector is effectively only constrained to change orthogonal to itself but is otherwise free to follow the linearised dynamics. Based on our earlier discussion, this means that it will quickly converge to the dominant direction of the tangent space and subsequently follow it. The first OTD basis vector is therefore covariant with the dynamics, i.e. is a solution of (2.2), while all other basis vectors are not, since they are constrained to the orthogonal complement of the previous basis vectors. The direct connection between the OTD basis vectors and covariant vectors is still lacking in the literature.

The OTD basis vectors, with exception of the first one, are therefore not amenable to physical interpretation and it is more useful to consider the most unstable directions within the OTD subspace spanned by the columns of $Q$, which contain information about the physical instabilities of the flow. In order to recover the $r$ most unstable directions $U$ in phase space, we can compute an eigendecomposition of the reduced operator $L_r=V\varLambda V^{-1}$ where the columns of $V\in \mathbb {C}^{r \times r}$ contain the eigenvectors of $L_r$ and $\varLambda$ is a diagonal matrix with the corresponding eigenvalues $\lambda _i$ ordered by decreasing real part, i.e. $\mathrm {Re} (\lambda _1) \ge \mathrm {Re} (\lambda _2) \ge \cdots \ge \mathrm {Re} (\lambda _n)$. As pointed out in Babaee & Sapsis (2016), the corresponding directions $U$ are then obtained by projecting the basis vectors $Q=[q_1, \ldots , q_r]$ as

(2.10)\begin{equation} U = Q V. \end{equation}

Furthermore, we can compute the instantaneous growth rates $\sigma _i$ in the subspace that is obtained by computing the eigenvalues of the symmetrised operator $L_\sigma = (L_r + L_r^{\textrm {T}})/2$ that contain information about the instantaneous non-normal growth potential in the subspace that is generated if the eigendirections of $L_r$ are not orthogonal. In particular, the dominant eigenvalue $\sigma _{max}$ is the numerical abscissa and corresponds to the maximum growth rate in the subspace.

We will in the following use the term ‘OTD modes’ only to refer to the projected OTD basis vectors $u_i$, $i = 1,\ldots , r$, that are the columns of $U$.

2.3. OTD basis and finite-time Lyapunov exponents

Lyapunov exponents are a central tool in the analysis of dynamical systems and describe the asymptotic rate of separation in time of two points initially infinitesimally close in phase space (Wolf et al. 1985). Dealing with finite time horizons and transient phenomena, the finite-time Lyapunov exponents can be used instead. They are defined as the eigenvalues of the Cauchy–Green strain tensor $C_{t_0}^{t}$ defined as

(2.11)\begin{equation} C_{t_0}^{t} = (\boldsymbol{\nabla} F_{t_0}^{t})^{\textrm{T}} \boldsymbol{\nabla} F_{t_0}^{t}. \end{equation}

The Cauchy–Green strain tensor is symmetric positive definite and thus has an orthogonal eigenbasis $\xi _i$ with positive eigenvalues $\mu _i$.

While it follows from the equivalence between (continuously orthonormalised) Gram–Schmidt vectors and the OTD basis vectors that the OTD modes ultimately converge to the dominant eigendirections of the asymptotic Cauchy–Green tensor, it was independently shown by Babaee et al. (2017) that this convergence is exponentially fast (under mildly restrictive conditions) and that it is possible to compute the finite-time Lyapunov exponents (FTLEs) as a byproduct of the computation of the OTD basis at negligible extra cost, independently of the dimensionality of the considered system.

In this work we use the method for the computation of FTLEs given in Blanchard & Sapsis (2019a) that uses a classic result for the relation between FTLEs and Gram–Schmidt vectors. The leading FTLEs can then be computed as

(2.12)\begin{equation} \mu_i = \frac{1}{t-t_0} \int_{t_0}^{t} \langle Lq_i,q_i \rangle \,{\textrm{d}}\tau = \frac{1}{t-t_0} \int_{t_0}^{t} (L_{r})_{ii} \,{\textrm{d}}\tau. \end{equation}

If the OTD modes are used to compute the dominant FTLEs, a general guideline is to use the Blanchard & Sapsis formulation and recompute the FTLEs considering subsets of $L_r$. This is possible due to the hierarchical nature of the formulation. The FTLEs can be computed online, i.e. alongside the base flow trajectory, or offline in a separate computation, provided the reduced operator is sampled with sufficiently high frequency.

The FTLEs computed over a period of a limit cycle are equal to the real part of the Floquet exponents, i.e. the linear temporal growth rates over one period (Huhn & Magri 2020). In the rest of the paper, the FTLEs $\mu _i$ are computed exclusively over a period of the limit cycle and are therefore synonymous with the real part of the corresponding Floquet exponent. Note that the FTLEs and the Floquet exponents can be computed using the OTD framework because they depend only on the subspace and do not require the computation of covariant vectors. Since the OTD basis vectors (except the first) are not covariant with the dynamics, they do not automatically generate the Floquet vectors which are the special case of covariant vectors for time-periodic flows (Kuptsov & Parlitz 2012).

2.4. Computational cost and choosing the OTD subspace dimension $r$

The computational cost of computing the OTD modes was discussed in Babaee et al. (2017); we summarise the results here to motivate the choices for the general purpose implementation of the OTD modes in Nek5000 that includes a spatially localised version.

In general, the computation of an $r$-dimensional OTD basis requires the computation of the base flow trajectory (2.1) coupled with $r$ solutions of the linearised problem with an added forcing term (2.4), regardless of the chosen formulation. All of these computations involve the full system and are typically of similar cost. While the base flow and the OTD basis can be computed sequentially for time-invariant systems using the Blanchard & Sapsis formulation (2.9), it should be noted that, in the case of a time-dependent dynamical system, (2.4) involves a linearisation around the time-dependent base flow implying that all $r+1$ equations need to be solved simultaneously or the data saved for reuse. If the full system is high-dimensional, especially the simultaneous alternative can be very memory intensive.

One possibility to reduce the cost of the method, both in terms of computation time and storage, is to restrict the linear solves to part of the base flow domain. This option is discussed in more detail in § A.4 together with a description of an efficient implementation in Nek5000 and its validation. For the flow case considered in this work, however, the base flow is analytical (see § 4) and the classical version is employed.

While the appropriate size of the OTD subspace is dependent on the particular application and a general rule is elusive, a few considerations on the nature of the system at hand can guide the choice in practice.

Blanchard et al. (2018), who use the OTD subspace to control linear instabilities, suggest choosing the number of OTD basis vectors $r$ such that

(2.13)\begin{equation} r \geq \max (\mbox{dim } \mathcal{E}_q, \mbox{dim } \mathcal{E}_q^{s}), \end{equation}

where $\mathcal {E}_q$ and $\mathcal {E}_q^{s}$ refer to the unstable parts of the eigenspaces of the full operator $L$ and its symmetric part $(L+L^{\textrm {T}})/2$, respectively. This choice is motivated by the requirement to span the directions responsible for both modal ($\mathcal {E}_q$) and, in the case of non-normal operators, non-modal ($\mathcal {E}_q^{s}$) instabilities.

It is worth mentioning that the most challenging situation for the OTD modes is the case of eigenvalue crossings in which the most unstable direction is temporarily ambiguous. It is therefore advisable to compute at least one extra mode to mitigate the effects of such crossings on the results of interest.

2.5. Boundary and initial conditions for the OTD modes

The OTD equations require the same boundary conditions as the linearised problem. The simplest approach is to initialise all perturbation fields with random noise and orthonormalise the vectors by a Gram–Schmidt process to form an orthonormal basis. Note that physical boundary conditions of the problem do not need to be satisfied from the start, but that the OTD basis will quickly adapt, as pointed out in Babaee & Sapsis (2016). The advantage of random noise initialisation is that it ensures coverage of the broadest possible range of initial perturbation frequencies for the flow to pick up on, thus circumventing the need for prior knowledge of the structure of the instabilities of interest. This is the approach followed in most of this work.

In general, the choice of initial conditions is heavily case dependent and in concrete applications it is often beneficial to include existing knowledge of the system into the initial conditions to accelerate convergence of the OTD basis.

In particular, if several of the dominating eigenvalues of the system have a similar magnitude or the system has several eigenvalues close to the real axis (neutral stability), the convergence to the OTD subspace can be slow. In view of these difficulties, an initialisation with random noise is generally not a good choice in order to achieve fast convergence; Instead, an educated initial guess can decisively reduce the time needed for the OTD subspace to align with the dominant directions in the tangent space.

Another possibility is to initialise the OTD subspace with the leading eigenvectors of the symmetric part of the full operator $(L+L^{\textrm {T}})/2$, thus aligning it with the instantaneously fastest growing directions. The trade-off with this strategy is that it requires the computation of (part of) the spectrum of the full operator which can be computed using variants of the Arnoldi algorithm even for large systems.

3. Governing equations

The focus of the present work is to illustrate the OTD framework by applying it to incompressible fluid flow that is governed by the non-dimensional incompressible Navier–Stokes equations

(3.1)\begin{equation} \left.\begin{gathered} \frac{\partial \boldsymbol{U}_b}{\partial t} ={-} (\boldsymbol{U}_b \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U}_b - \boldsymbol{\nabla} p_b + \frac{1}{\mbox{Re}} \nabla^{2} \boldsymbol{U}_b + f_b , \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U}_b = 0, \end{gathered}\right\} \end{equation}

where $\boldsymbol {U}_b = (U_x,U_y,U_z)^{\textrm {T}}$ is the base flow velocity vector, $p_b$ is the pressure, $f_b$ is an external forcing term and $Re = U_\infty L/\nu$ is the Reynolds number based on the reference velocity $U_\infty$, the reference length scale $L$ and the kinematic viscosity $\nu$, supplemented by the appropriate boundary conditions.

The linearised Navier–Stokes equations governing the evolution of an infinitesimal perturbation $\boldsymbol {u} = (u,v,w)$ with the associated pressure field $p$ on top of a base flow $\boldsymbol {U}_b$ are given by

(3.2)\begin{equation} \left.\begin{gathered} \frac{\partial \boldsymbol{u}}{\partial t} ={-} (\boldsymbol{U}_b \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} - (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{U}_b - \boldsymbol{\nabla} p +\frac{1}{\mbox{Re}} \nabla^{2} \boldsymbol{u} + f, \\ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gathered}\right\} \end{equation}

where the same non-dimensionalisation as for the nonlinear equations is applied. For our analysis we allow time-dependent solutions of the Navier–Stokes equations as base flow.

When (3.1)–(3.2) are solved numerically, they need to be discretised in space and time. In order to apply the OTD formalism to the resulting discretised Navier–Stokes equations, we transform them into dynamical system form. In Nek5000, this is achieved by projecting the solutions onto a divergence-free space, thus removing the explicit dependence on pressure. Combining the three components of the $i$th velocity perturbation into a single vector $q_i$, (3.2) can be expressed for each perturbation as a forced dynamical system

(3.3)\begin{equation} \frac{\partial q_i}{\partial t} = L_{{LNS}}q_i + f_i,\quad i = 1, \ldots, r, \end{equation}

where $L_{{LNS}}$ is the discrete linearised operator. Note that $L_{{LNS}}$ is intrinsically time dependent via the base flow.

Comparing (3.3) with the Blanchard & Sapsis formulation (2.9) shows that the additional constraint for the OTD equations can be introduced directly via the external forcing term

(3.4)\begin{equation} f_i ={-} \langle L_{{LNS}} q_i, q_i \rangle - \sum_{j=1}^{i-1} (\langle L_{{LNS}} q_i, q_j \rangle + \langle L_{{LNS}} q_j, q_i \rangle) q_j, \end{equation}

where the inner product $\langle \cdot , \cdot \rangle$ is the standard energy norm over the computational domain $\varOmega$ where the energy of a perturbation $q_i = (u,v,w)$ is computed as

(3.5)\begin{equation} \langle q_i,q_i \rangle = \int_\varOmega [u^{2} + v^{2} + w^{2}]\,{\textrm{d}}\varOmega. \end{equation}

The standard formulation of the OTD equations (2.4) can be obtained in a similar fashion.

Note that the choice of inner product defines the resulting OTD basis and the derived quantities. Therefore, especially in situations like compressible flows where no obvious physical choice exists (Colonius et al. 2002), care must be taken in the definition of the inner product and the subsequent interpretation of the OTD modes.

4. Pulsating Poiseuille flow

4.1. Flow case

In order to illustrate the implementation of the OTD methodology in an unsteady setting, we consider pulsating plane Poiseuille flow. This flow case exhibits a parallel, temporally periodic, streamwise and spanwise independent base flow in the axial direction of the form $\boldsymbol {U}_p = (U_x(\kern0.07pt y,t),0,0)$ that satisfies the incompressible Navier–Stokes equations. To simplify the problem, we consider only pulsations with a single base frequency $\varOmega$. Any periodic pulsations can be analysed in a similar fashion by considering the corresponding Fourier series. The pulsating streamwise velocity component relative to the base frequency $\varOmega$ can be expressed as the sum of a parabolic steady component (superscript 0) and purely oscillatory component (superscript osc)

(4.1) \begin{equation} U_x(\kern0.07pt y,t) = U_x^{(0)}(\kern0.7pt y) + U_x^{osc}(\kern0.07pt y,t). \end{equation}

This velocity profile is generated by a spatially uniform and time-periodic streamwise pressure gradient $-G_x(t)$ and corresponds to the mass-flow rate $Q(t)$ of similar form

(4.2)\begin{equation} Q_x(t) = Q_x^{(0)} + Q_x^{osc}(\kern0.07pt y,t). \end{equation}

Some details on the formulae for the general case can be found in Von Kerczek (1982) and Pier & Schmid (2017).

The steady and unsteady velocity components are respectively given by

(4.3)\begin{equation} \left.\begin{gathered} U_x^{(0)} (\kern0.07pt y) = \frac{3}{2} \frac{Q_x^{(0)}}{2h} \left( 1 - \left( \frac{y}{h} \right)^{2} \right) ,\\ U_x^{osc} (\kern0.07pt y,t) = \frac{ Q_x^{osc} (\kern0.07pt y,t)}{2h} = \frac{Q_x^{(1)}}{2h} \mathrm{Re} \left( \mathcal{W} \left( \frac{y}{h},Wo \right) \textrm{e}^{\textrm{i}\varOmega t} \right), \end{gathered}\right\} \end{equation}

where $Q_x^{(1)}$ is the amplitude of the oscillating component, Wo is the Womersley number relating the channel half-height $h$ to the thickness of the oscillating boundary layer $\delta =\sqrt {\nu /\varOmega }$ defined as

(4.4)\begin{equation} Wo = \frac{h}{\delta} = h \sqrt{\varOmega/\nu}, \end{equation}

and $\mathcal {W}(\xi ,Wo)$ is a function defining the profile of the unsteady velocity component given by

(4.5)\begin{equation} \mathcal{W}(\xi,Wo) = \left( \frac{\cosh(\sqrt{i}\xi Wo)}{\cosh(\sqrt{i} Wo)} - 1\right)/ \left( \frac{\tanh(\sqrt{i} Wo)}{\sqrt{i} Wo} - 1 \right), \end{equation}

where $\xi = y/h$ and we have used $\sqrt {i}= (i+1)/\sqrt {2}$.

Some authors use $\beta = Wo/\sqrt {2}$ as the frequency scale but we will follow Pier & Schmid (2017) to facilitate comparison. The pulsation frequency can be computed as $\varOmega = Wo^{2}/\mbox {Re}$ which, in turn, leads to a pulsation period of $T_0 = 2{\rm \pi} \mbox {Re}/Wo^{2}$.

Following Pier & Schmid (2017), we introduce $\tilde {Q}$, the normalised fluctuating part of the unsteady flow rate $Q(t)$, given by

(4.6)\begin{equation} \tilde{Q} = \frac{Q_x^{(1)}}{Q_x^{(0)}} \quad \mbox{such that}\ Q(t) = Q_x^{(0)} ( 1 + \tilde{Q} \cos (\varOmega t) ), \end{equation}

where we choose $Q_x^{(n)} \in \mathbb {R}$ for $n=-1, 0, 1$ without loss of generality. Note that $Q_x^{(1)}$ is twice the corresponding value in Pier & Schmid (2017) to yield an identical definition of $\tilde {Q}$.

From the equations above it is clear that setting $Q_x^{(1)} = 0$ recovers the solution for steady plane Poiseuille flow with

(4.7)\begin{equation} U_x(\kern0.07pt y) = U_c \left( 1 - \left( \frac{y}{h} \right)^{2} \right), \end{equation}

where we have used that $Q_x^{(0)}= \frac {2}{3}U_c \cdot 2h$. The associated steady pressure gradient is

(4.8)\begin{equation} G_x^{(0)} = \frac{2U_c\nu}{h^{2}} . \end{equation}

Figure 1 shows an example of the base flow variation over one oscillation period including the characteristic boundary layer thickness $\delta$ for reference. Note the complex interplay of reverse flow and inflection points in the local profiles. The stability of the inflection points is based on the local Fjørtoft criterion (Schmid & Henningson 2001). Note that, while reverse flow only appears at higher pulsation amplitudes, inflection points will occur along the profile for virtually all pulsating cases.

Figure 1. Schematic representation of the components of pulsating Poiseuille flow for $\tilde {Q}=1.0$ and $Wo=10$ over one period ($T_0 = 471.2$). The complete profile (c) is plane Poiseuille flow (a) superimposed with an oscillating flat Stokes layer (b). (d) Shows a schematic representation of one pulsation cycle including the regions of reverse flow (shaded in blue), the Stokes layer thickness (dashed black lines) as well as the local base flow profiles (thick black lines, the number indicates the corresponding time $t/T_0$). The red lines indicate the location of inflection points in the base flow profile and their stability (full) or instability (dashed) based on the local Fjørtoft criterion.

We have chosen to normalise the velocity components with the centreline velocity of plane Poiseuille flow $U_c$ and the channel half-height $h$ which differs from the normalisation adopted in Pier & Schmid (2017). To facilitate comparison, we briefly summarise the relation to the normalisation adopted in this work for a few central parameters in table 1.

Table 1. Conversion table between different choices of non-dimensionalisation in Pier & Schmid (2017), Von Kerczek (1982) and the present work.

4.2. Local stability analysis

Parallel and geometrically homogeneous shear flows such as plane Poiseuille flow can be analysed by defining a fundamental pair of real wavenumbers $(\alpha ,\beta )$ in the stream- and spanwise directions, respectively, which uniquely describe the problem.

We consider wave-like perturbations of the form

(4.9)\begin{equation} \begin{pmatrix} u(x,y,z,t) \\ v(x,y,z,t) \\ w(x,y,z,t) \\ p(x,y,z,t) \end{pmatrix} = \begin{pmatrix} \tilde{u}(\kern0.07pt y) \\ \tilde{v}(\kern0.07pt y) \\ \tilde{w}(\kern0.07pt y) \\ \tilde{p}(\kern0.07pt y) \end{pmatrix}_{\alpha,\beta} \exp({\textrm{i}\alpha x + \textrm{i}\beta z - \omega t}), \end{equation}

where $u,v,w$ are the velocity components, $p$ is the pressure field, $\tilde {u},\tilde {v},\tilde {w}, \tilde {p} \in \mathbb {C}^{n}$ are their respective Fourier transforms in the homogeneous directions and $\omega \in \mathbb {C}$ is the complex growth rate; $\mathrm {Re}(\omega ) = \omega _r$ is the (real) temporal growth rate and $\mathrm {Im}(\omega ) = \omega _i$ is the angular frequency. Instead of $\omega _i$, the phase speed computed as $c= \omega _i/\alpha$ is often used in this work.

Introducing this ansatz into (3.2) in the absence of external forcing and removing the pressure dependence by combining the two equations using the wall-normal vorticity $\eta$, the linearised Navier–Stokes equations for a fixed base flow profile can be identically recast as an eigenvalue problem for the $\omega$ as

(4.10)\begin{equation} - \omega \tilde{q} = \mathcal{L}_{{OS}} \tilde{q}, \end{equation}

where $\tilde {q}$ is the state vector containing the Fourier transforms of the wall-normal velocity and vorticity

(4.11)\begin{equation} \tilde{q} = \begin{pmatrix} \tilde{v} \\ \tilde{\eta} \end{pmatrix}_{\alpha,\beta}. \end{equation}

This is the standard form of the well-known Orr–Sommerfeld (OS)/Squire (SQ) equations. The full derivation of the OS/SQ equations can be found in Schmid & Henningson (2001).

The reference spectra for plane Poiseuille flow are obtained for the linear operator in (4.10) using a Chebyshev collocation discretisation with a wall-normal resolution of 196 points solved in MATLAB.

The eigenvalue problem (4.10) that is at the heart of the local analysis, by construction, does not allow for time dependence of the base flow. The validation of an implementation of the OTD framework against solutions of the OS/SQ equations is therefore necessarily restricted to steady cases in which both methods yield the same results. Furthermore, care needs to be taken to ensure that both methods solve the same local stability problem in order to be able to compare results from Nek5000 with results of (4.10) or Pier & Schmid (2017) who also use the normal mode ansatz in their analysis. Because of the finite extent of computational domain in the homogeneous direction(s) needed for Nek5000, we cannot a priori set the stream- and spanwise wavenumbers $\alpha$ and $\beta$ to define the local problem. Hence, numerical noise may lead to energy growth in other wavenumbers that fit the computational box and contaminate the results. Therefore, these spurious modes are projected out at every timestep.

Note also that (4.10) is expressed in (complex) Fourier space whereas Nek5000 works in physical space. This means that, since the OS/SQ equations do not have steady eigenmodes, each complex OS/SQ mode will correspond to two complex conjugate modes computed with Nek5000. For this reason, all computations in this work are performed with an even number of modes and we will count the complex conjugates as a single mode for clarity.

5. Numerical set-up in Nek5000

The present implementation is done in the high-order spectral element code Nek5000 Fischer et al. (2008) solving the incompressible Navier–Stokes equations. In the spectral element method (SEM), the computational domain is divided into a conformal grid of deformable, hexahedral subdomains (quadrilateral in two dimensions). On each element, the velocity solution is expanded in terms of Lagrange polynomial basis functions of degree $N$ in each spatial direction at Gauss–Lobatto–Legendre quadrature points (for the pressure, polynomials of degree $N-2$ on Gauss–Lobatto points are used), which is referred to as the $\mathbb {P}_N-\mathbb {P}_{N-2}$ formulation. Note that the choice of $N$ affects both the number of degrees of freedom and the rate of spatial convergence of the solution. Generally, $N\geq 7$ is recommended and often suffices, but convergence tests are performed with $N=15$. The SEM approach combines the high accuracy of spectral methods with the geometric flexibility of finite element methods and is therefore easily extensible to more complex geometries. Albeit, a simple geometry has been chosen in the present work for illustrative purposes.

Timestepping is performed via the third-order accurate backward difference scheme for both linear and nonlinear simulations where the viscous terms are treated implicitly and the nonlinear and forcing terms are computed explicitly using a third-order accurate extrapolation scheme with over-integration. To ensure an accurate solution, the timestep was set such that the Courant–Friedrichs–Lewy number was at most 0.4 relative to the base flow. The timestep carries over directly to the linear simulations.

To simplify comparisons with the results from the OS/SQ equations, a domain length $L$ of $2{\rm \pi}$ was chosen in the homogeneous directions. The set of wavenumbers that can be supported by the system is therefore $\alpha ,\beta \in \mathbb {N}$. The structure of the mesh used in the stability analysis of two-dimensional (2-D) plane Poiseuille flow is shown in figure 2. The corresponding mesh for the 3-D case used to validate the code has the same resolution in both homogeneous directions (see the Appendix in § A). Since the solutions are periodic, they are in principle independent of the resolution in the homogeneous direction. Nonetheless, it was found that reasonable resolution was necessary for numerical stability of the linearised solver and 8 elements were placed in these directions. The wall-normal resolution is considerably higher, especially towards the walls where the base flow gradient is largest. The simulations were run with polynomial order $N=9$ with solver tolerances set to $10^{-12}$ for both velocity and pressure equations. No-slip boundary conditions are imposed at the wall and periodic conditions in the homogeneous direction(s), as indicated in figure 2. Since the base flow is analytically known for both steady and unsteady configurations it is not computed explicitly but instead supplied directly to save computational effort.

Figure 2. Two-dimensional spectral element mesh for Nek5000. The base flow profile for plane Poiseuille flow is shown in black. The colours indicate the boundary conditions. The top and bottom walls are no slip (red) and periodic conditions are imposed on the streamwise boundaries (blue).

Further details on some numerical aspects of the solution of the OTD equations in Nek5000 and the validation for plane Poiseuille flow can be found in the Appendix.

6. Numerical results

6.1. OTD modes in pulsating Poiseuille flow

The addition of the pulsatile component to the parabolic base flow profile of plane Poiseuille flow changes the dynamics of the system dramatically. The pulsation frequency $\varOmega$, with which the system is ultimately synchronised, introduces another time scale into the system, the pulsation period $T_0$, separating the long-term intercyclic behaviour from the intracyclic growth and decay within a pulsation cycle. Note that in a time-periodic flow the limit cycle once all initial transients have passed will have the same temporal periodicity as the base flow. Since we consider mainly growth rates over time, the corresponding traces are purely periodic in the limit cycle and we therefore refer to it as the periodic regime. The OTD modes instantaneously span the subspace of most unstable directions and therefore allow us to study the initial phase of transient growth, the intracyclic growth over a base flow period as well as the intercyclic behaviour dictated by the Floquet exponents.

6.1.1. Long-term evolution and Floquet exponents

Floquet analysis is usually performed separately from classical linear stability analysis since it requires the equations of motion to be recast to remove the intracyclic variations. Using OTD modes, the computation of FTLEs and Floquet exponents is a byproduct of the method without significant overhead. Using this framework, the first Floquet exponent was computed for two-dimensional cases with $\mbox {Re}=7500$, $\alpha =1$ for different combinations of Womersley number $Wo$ and relative amplitude of the oscillating flow rate $\tilde {Q}$. Since only the leading Floquet exponent is computed, these simulations were run with only one complex mode. The corresponding linear temporal growth rates are shown in figure 3. The values for $Wo$ and $\tilde {Q}$ were chosen such as to cover the parameter range used in Pier & Schmid (2017) including unstable and transiently stable cases and also extend the range to higher relative mass-flow rates. The range of Womersley numbers is related to experimental values reported for physiological flows (Pier & Schmid 2017) and lies in the parameter range for pulsating Poiseuille flow where the time scales of the pulsation and the Tollmien–Schlichting (TS) waves associated with the steady base flow are of similar magnitude, i.e. the pulsations are neither fast nor slow (Davis 1976) so that they are expected to noticeably interact with the instability waves. It is also close to the strongest pulsatile stabilisation at moderate pulsation amplitudes ($Wo\approx 28$, Singer, Ferziger & Reed 1989). Note that plane Poiseuille flow is linearly unstable for the chosen parameter combination ($\textrm {Re}= 7500$). Three Womersley numbers were chosen exhibiting very different stability behaviour for increasing values of $\tilde {Q}$. These results are in very good agreement with reference data provided by B. Pier (private communication).

Figure 3. Temporal intercyclic growth rates for different values of $Wo$ and $\tilde {Q}$. Full lines adapted from Pier & Schmid (2017) (data computed up to $\tilde {Q}\leq 0.6$ received in private communication), filled circles computed in Nek5000. The black symbol is the steady configuration. Note that the reference values are scaled to match the normalisation used in this work. (a) Is a close-up of the region of small pulsation amplitudes with quadratic dependence on $\tilde {Q}$.

While low pulsation frequencies ($Wo=10$) are destabilising in comparison with the steady case at the same Reynolds number for all values of $\tilde {Q}$ investigated, we observe that the destabilisation seems to saturate beyond relative pulsation mass-flow rates of $\tilde {Q}=0.5$ and even slightly decrease again at the highest values of $\tilde {Q}$. For higher frequencies, the behaviour is very different. For intermediate pulsation frequencies ($Wo=18$), an increase in pulsation amplitude leads to a nearly monotonic increase in intercyclic stabilisation up to a maximum cycle-to-cycle decay rate of $\mathrm {Re}(\mu _1) = -0.038$ for $\tilde {Q}=1.00$. For high pulsation frequencies ($Wo=25$), although pulsations are stabilising overall, the relationship between pulsation amplitude and intercyclic growth rate is highly non-monotonic. For this parameter combination, we observe a maximum stabilisation of $\mathrm {Re}(\mu _1) = -0.031$ at around $\tilde {Q}=0.38$ that subsequently reduces considerably to the end of the considered amplitude range at $\tilde {Q}=1.0$ where similar values as for $\tilde {Q}=0.1$ are observed.

6.1.2. Intracyclic growth rate modulation over a pulsation period

While Floquet exponents shed light on the long-term fate of the system, they do not allow for the analysis of the intracyclic variations that ultimately dictate the cycle-to-cycle growth or decay of disturbances. Using the OTD framework, we can trace the eigenvalues of the reduced operator $L_r$ and the associated instantaneous OTD modes during the pulsation. In the following, we will repeatedly refer to steady values with which we mean the corresponding values for the steady case at the same Reynolds number, i.e. $\tilde {Q}=0$, without distinction of whether these are computed using Nek5000 or MATLAB.

The evolution of the real part of the eigenvalues as well as the numerical abscissa of $L_r$ for the pulsating case is shown in figure 4(a) for $r=6$ for $Wo=25$ and moderate pulsation amplitudes $\tilde {Q}=0.2$. Comparing the signal with the most unstable part of the spectrum of the corresponding OS operator (dashed lines), we see that, while the growth rates of all eigenvalues are affected by the pulsation during the transient response, the periodic regime is characterised for each mode by either an alignment of the eigenvalue modulation with the pulsation frequency or the convergence to the respective steady value. Doubling the subspace size to $r=12$, shown in figure 4(b), we observe a faster convergence of the least stable modes to the periodic regime because the more stable modes shield them from the dynamics outside of the subspace that may lead to eigenvalue crossings. These are general features that are largely independent of $Wo$, $\tilde {Q}$ and the subspace size for moderate values of $r$. The situation changes for very large subspace sizes; see § 6.1.6 for details.

Figure 4. Instantaneous real eigenvalues and numerical abscissa of the reduced operator $L_r$ for 2-D pulsating Poiseuille flow with $\mbox {Re}=7500$, $\alpha =1$, $\tilde {Q}=0.2$ and $Wo=25$; (a) $r=6$, (b) $r=12$. At the points marked $t_1,t_2$ and $t_3$ in (b), the instantaneous structure of the corresponding mode is shown in figure 5.

Figure 5 shows the velocity components of several OTD modes at different instants of the pulsation cycle in the periodic regime, labelled according to the order of the eigenvalues of their steady counterparts. Mode $u_1$ is reminiscent of the TS wave in steady Poiseuille flow that is modulated in the near-wall region, especially at the moment of maximum amplification (figure 5a). When the corresponding mode is most damped (figure 5b), the spatial structure is strongly affected by the pulsation close to the wall. Note that, since there is non-normal growth potential at all times, the most unstable direction will be a superposition of different modes, as can be seen in the region close to the centreline for mode $u_1$. Figure 5(c) shows one of the centre modes (mode $u_3$) that is very similar to its steady counterpart. Although the spatial structure is somewhat altered due to non-normality, the spatial structure of the OS mode is clearly visible and the growth rate is unaffected. Finally, figure 5(d) shows the spatial structure of mode $u_4$ when it is most amplified, resembling the second least stable mode of the T-branch of the OS spectrum.

Figure 5. Instantaneous $u$ and $v$ components of $u_i$ at different times $t_j$ along the pulsation cycle in the periodic regime as marked in figure 4(b); (a) $u_1$ at maximum amplification ($t_1 = 428$), (b) $u_1$ at maximum damping ($t_2 = 473$), (c) $u_3$ at $t=t_1$, (d) $u_4$ at maximum amplification ($t_3 = 432$). The parameters are the same as in figure 4(b). Videos of the evolution of $u_1$ and $u_4$ over the pulsation cycle are available in the supplementary material as Movie 1 and Movie 2 available at https://doi.org/10.1017/jfm.2021.743, respectively.

In order to characterise the different modes more precisely, it is instructive to consider also the phase speed i.e. the imaginary part of the eigenvalue since $\alpha =1$. Combining the traces of the instantaneous real and imaginary parts of the eigenvalue in the periodic regime and comparing them for different values of $\tilde {Q}$ exhibits a clear pattern shown in figure 6 for the two dominant oscillating eigenvalues for $Wo= 10$ (a,b) and $Wo= 25$ (c,d). We observe that the eigenvalue excursions are more pronounced for large pulsation amplitudes, especially at low frequency, and are generally more erratic for higher pulsation frequencies. For $Wo=10$, low pulsation amplitudes lead mostly to a comparatively small variation of the growth rate at relatively constant phase speed and then abruptly begin with very large excursions of both growth rate and phase speed for amplitudes beyond $\tilde {Q}=0.1 - 0.15$. Nevertheless, as $\tilde {Q}$ approaches 0 (steady state), the periodic orbits contract around the most unstable eigenvalues of the A-branch, confirming the identification based on the spatial structure of the modes in figure 5. Comparing the eigenvalue orbits for the same pulsation frequency, we see that both oscillating modes exhibit similar trends throughout the tested amplitude range, with mode 4 being more stable throughout than mode 1, as in the steady case.

Figure 6. Periodic orbits in the complex plane of the instantaneous eigenvalues $\lambda _1$ and $\lambda _4$ of $L_r$ for different pulsation amplitudes $\tilde {Q}$ at $\mbox {Re}=7500, \alpha =1, \beta = 0$ compared with the corresponding OS spectrum. The shaded area indicates linear instability. The rows correspond to the two different values of $Wo$ and every column corresponds to a different mode. For $Wo = 10$: $\lambda _1$ (a) and $\lambda _4$ (b). For $Wo = 25$: $\lambda _1$ (c) and $\lambda _4$ (d).

The two distinct behaviours of the modes in the pulsating case in comparison with the steady counterpart is linked to the spatial localisation of the modes in relation to the base flow modulation. The pulsations superimposed onto the parabolic profile exhibit the most variation close to the walls leading to the eigenvalue excursions of modes 1 and 4 that have similar behaviour. Close to the centreline of the channel, the pulsating velocity component is flat so that centre modes see very little variation of the base flow curvature in comparison with the steady case. Their growth rates are largely unaffected. It is a general feature that all modes on the A-branch of the OS spectrum are more strongly affected by the pulsations than the modes on the P-branch, especially the closer they are to the centreline.

Also in terms of the phase speed the eigenvalues have a much less erratic behaviour as the pulsation amplitude increases. While the phase speeds of all eigenvalues oscillate with the mean flow amplitude, only the eigenvalues from the A-branch noticeably deviate from a purely sinusoidal variation in phase speed. Towards the upper limit of the amplitude range ($\tilde {Q} = 1.0$) and low frequencies, where we observe the largest excursions of the modulated eigenvalue traces, we also begin to see the pulsations slightly affect the growth rates of the centre modes. This is because, at these extreme amplitudes, the modifications of the local base flow profiles due to the pulsations are no longer confined to the near-wall regions. This effect is more pronounced for lower frequencies since the associated oscillating boundary layer is thicker and thus reaches further into the channel.

The periodic orbits give a good overview of the eigenvalue excursions in the complex plane but are ill suited to compare different eigenvalue traces at specific instants along the pulsation cycle since the time dependence is implicit. In this context, it is more useful to consider the intracyclic modulation of the real and imaginary parts of mode 1 (modulated TS wave) in the periodic regime separately over time to link the instability characteristics to flow features. The eigenvalue traces are shown in figure 7 at $Wo=10$ 7(a,b) and $Wo=25$ 7(c,d) for the same values of $\tilde {Q}$ as in figure 6.

Figure 7. Intracyclic modulation of the most unstable eigenvalue ($\lambda _1$) in the periodic regime for different pulsation amplitudes $\tilde {Q}$ at $\mbox {Re}=7500, \alpha =1, \beta = 0$. The horizontal dashed line indicates the leading eigenvalue of the corresponding OS spectrum; (a,b) $Wo=10$, (c,d) $Wo=25$.

Comparing the eigenvalue traces with the corresponding Floquet exponents (cf. figure 3), it is evident that the growth rates over the pulsation cycle are considerably larger than the net cycle-to-cycle growth rate. This becomes particularly clear when comparing them with the temporal growth rate of the corresponding steady case (dashed black line). The steady growth rate is of roughly the same order of magnitude as the leading Floquet exponent at moderate pulsation amplitudes which themselves are 1–2 orders of magnitude smaller than the peak growth and decay rates of even very low pulsation amplitudes ($\tilde {Q}=0.05$). These large intracyclic growth rates that are sustained for a large fraction of the period lead to the well-documented large perturbation amplitudes during a pulsation cycle (Pier & Schmid 2017).

Furthermore, it can be noted that, while the growth and decay rates tend to increase with the pulsation amplitude for all values of $Wo$, the peak growth rates saturate at around 0.11 and damping rates rarely exceed $-0.15$ when the pulsation amplitudes exceed values of roughly $\tilde {Q} = 0.5$. The maximum growth rates achievable with pulsations are therefore lower than the peak values exploiting transient growth mechanisms. Due to the fact that, unlike transient growth mechanisms that act over short time scales, the high intracyclic growth rates are maintained over long periods of time, they lead to higher energy growth overall. Since the extremal growth rates saturate, the continued variation of the Floquet exponents for higher pulsation amplitudes is mainly due to variations of the distribution of amplification and damping throughout the cycle.

Due to the large parameter space it is prohibitively expensive to give an exhaustive description of the effect of pulsation here and we will focus on two combinations $(Wo,\tilde {Q})$, namely $(10,0.5)$ and $(25,0.5)$, to illustrate the range of possible intracyclic variation and show some typical features of the stability of pulsatile Poiseuille flow. The streamwise velocity component of the OTD mode corresponding to $\lambda _1$ for these cases is shown in figures 8(b)–8(a) and 9(b)–9(a), respectively, at selected instants during the pulsation cycle showcasing prominent features. The simulations for this section were run with $r=12$. We only show mode 1 since mode 4 shows very similar behaviour and the remaining modes are centre modes that exhibit little variation over the pulsation cycle.

Figure 8. Instantaneous streamwise velocity component of $u_1$ at different instants over a pulsation cycle for $Wo=10$ and $\tilde {Q}=0.5$ in the periodic regime (a) compared with the base flow profiles, the location of inflection points and the boundary layer thickness (b) as well as the instantaneous growth rate of the mode (c). The velocity fields and base flow profiles are shown for $t/T_0 = 0.06, 0.15, 0.59, 0.75$ marked as $t_1, \ldots ,t_4$, respectively. The distinction between stable and unstable inflection points is as in figure 1. The computations were performed with $r=12$. A video of the evolution of $u_1$ over the pulsation cycle is available in the supplementary material as Movie 3.

Figure 9. Same as figure 8 but for $Wo=25$ and $\tilde {Q}=0.5$. The velocity fields and base flow profiles are shown for $t/T_0 = 0.06, 0.56, 0.62, 0.68$ marked as $t_1, \ldots ,t_4$, respectively. A video of the evolution of $u_1$ over the pulsation cycle is available in the supplementary material as Movie 4.

6.1.3. Stability characteristics at low Womersley numbers

For low frequency pulsations ($Wo = 10$) and low pulsation amplitudes, the growth rate variations increase in amplitude with $\tilde {Q}$. The phases of growth and decay happen at approximately the same time along the cycle, positive amplification rates existing in the range $t/T_0 \in [0.14, 0.6]$. As the pulsation amplitude is further increased, these phases start lagging behind, aligning with the appearance of backflow for large pulsation amplitudes. The amplification phase is well correlated with the existence of an unstable inflection point outside of the oscillating boundary layer and the snapshots of the modes during this phase (figure 8(a), $t_2$ to $t_4$) clearly show the localisation of the mode at the wall-normal height of the inflection points. In fact, as a second (stable) inflection point appears at the beginning of the acceleration phase, it is also followed by the mode ($t_3$). The structure of the amplified mode is largely unchanged over the considered range of pulsation amplitudes. The damping phase initiates as the inflection points collide and vanish ($t_4$). The observation that the mode follows the base flow inflection points only once they move out of the oscillating Stokes layer is likely linked to the fact that oscillating Stokes layers are known to have a stabilising effect on disturbances, decaying faster than in a quiescent fluid due to viscous dissipation (Davis 1976). The fact that regions very close to the wall do not promote instability has also been linked to low values of the wall-normal velocity ($v$) prevailing there that plays a key role in linear growth mechanisms (Obremski et al. 1969).

As $\tilde {Q}$ is increased past approximately 0.1, oscillations of the eigenvalue start appearing during the damping phase. As the amplitude is further increased, these oscillations intensify and finally coalesce, leading to strong damping at the beginning of the cycle that abruptly transitions to amplification. Snapshots of the $u$-component of mode 1 contrasting the spatial structure during the damping and the beginning of the amplification phase after are shown in figure 8(a) from $t_1$ to $t_2$, respectively. Instead of TS-like structures close to the wall that are not visible at all, the mode consists of vortices located half-way between the wall and the centre of the channel that are strongly tilted in flow direction. Their structure is comparable to Orr waves that, once they have tilted with the base flow and can no longer extract energy from it, also exhibit similar decay rates. The shift towards the centre of the channel is also reflected in the phase speed that increases dramatically when the mechanism picked up by the mode changes. The oscillations prior to the coalescence are related to intermittent changes of the growth mechanism from TS-like structures to Orr-like structures and back before the Orr waves shown in figure 8(a) at $t_1$ become dominant throughout the damping phase.

An interesting note is that the appearance of the oscillations that herald the shift of the mode as well as the increase of the maximal damping coincide with the limit between the cruising and ballistic regimes identified by Pier & Schmid (2017) in the saturated nonlinear regime. The latter is characterised by the complete disappearance of perturbations of the base flow over parts of the pulsation cycle. When damping is strongest, the leading OTD mode has negligible amplitude close to the wall. In fact, none of the leading 6 OTD modes pick up on TS-like structures close to the wall during the damping phase, indicating that these structures are even more damped. At these lower pulsation frequencies the pulsation period is considerably longer than the viscous time scale, giving the perturbations time to decay completely from the nonlinear amplitudes with the major energy supply from the linear mechanisms cut off. The effect is similar to the intermittent transition observed in direct numerical simulations of pulsating pipe flow by Xu & Avila (2018), in which turbulent puffs completely decay at low Womersley numbers.

6.1.4. Stability characteristics at high Womersley numbers

In comparison with the low frequency case that seems to fall into one of two categories of intracyclic behaviour, the variations of growth rate and phase speed are much more complicated at higher pulsation frequencies. This reflects on the intercyclic growth that is highly non-monotonic. The large changes in intracyclic behaviour for small changes in pulsation amplitude suggest that several competing mechanisms exist simultaneously.

Due to the higher frequency the oscillating boundary layer is very thin compared with the channel diameter and the phase lag between the boundary layer and the flow rate variation is more pronounced and leads to a more inflectional base flow profile. In fact, already for $\tilde {Q}=0.35$ the base flow profile exhibits inflection points at all times. Figure 9(b) shows the location of inflection points at $\tilde {Q} = 0.5$. Due to the phase lag of the boundary layer, the amplification and decay phases are shifted compared with the lower pulsation frequencies. A similar shift in the timing of the maximum decay is seen as $\tilde {Q}$ is increased, but no oscillations in the eigenvalue can be observed. The shape of the dominant mode never completely changes as is observed in figure 8(a) from $t_1$ to $t_2$, but remains similar to distorted waves localised at the wall, as shown in figure 9(a).

Also here, the amplification phase is generally correlated with the existence of two inflection points outside of the Stokes layer and the modes localise on them. The abrupt beginning of the amplification phase is related to the emergence of the second inflection point outside of the oscillating boundary layer that leads to the breakup of the strongly sheared vortex ($t_2$) into two separate vortices ($t_3$), each localising on one of the inflection points. The phase of the two vortices quickly aligns as the growth rate peaks to reestablish a single near-wall wave ($t_4$). The wall waves stay localised on the inflection points during the entire amplification phase. The beginning of the damping phase is again closely related with the collision and disappearance of the inflection points ($t_1$). Although the fact that pulsation and viscous time scales are similar would suggest that local properties of the base flow profile play a less important role for the stability characteristics, we see that the existence of inflection points outside of the oscillating boundary layer is intimately linked to the amplification phases even at high pulsation frequencies.

6.1.5. Pulsations and non-normality

The linearised Navier–Stokes operator is highly non-normal, which can lead to considerable transient growth due to the non-orthogonality of its eigendirections. In order to analyse the impact of non-normality on the stability characteristics, we need a measure for the non-normal growth potential, i.e. the proportion of numerical abscissa $\sigma _{max}$, including both modal and non-modal contributions, that cannot be explained by modal mechanisms alone. Since $\sigma _{max}$ always equals or exceeds the growth rate of the least stable eigenvalue $\lambda _1$, equality implying normality of the operator for which the non-normal growth potential is zero, we choose the difference $\sigma _{max}-\lambda _1$ as a measure for this potential (although other choices are possible). Obviously, the perturbations that will experience the growth rates given by $\sigma _{max}$ and $\lambda _1$, respectively, will not be the same in the non-normal case. Therefore, their difference is not a measure for the non-normal growth potential of a specific perturbation but instead a global indicator of the non-normal growth potential over all possible perturbations due to non-normality.

Comparing again figures 4(a) and 4(b), this time focusing on the numerical abscissa, we see that it follows the modulation of the oscillating eigenvalues. The non-normal growth potential is much less variable than the eigenvalues themselves and seems to change only little as the subspace size is doubled to $r=12$, even though we observe both large variations of the intracyclic growth rates and highly damped Floquet exponents compared with the steady case (see figure 3). Moreover, the fact that the general structure of the eigenvectors does not change dramatically over the pulsation cycle makes an incidental correlation of $\sigma _{max}$ and $\lambda _i$ improbable. This might indicate that the numerical abscissa is already saturated with few OTD modes. This is in fact not the case, as can be seen in figure 10(a), showing the evolution of the numerical abscissa as the OTD subspace dimension is increased up to $r=100$.

Figure 10. Variation of the numerical abscissa $\sigma _{max}$ of the reduced operator $L_r$ for different choices of $r$ at $Wo=25$ and $\tilde {Q}=0.2$ in the limit cycle. Blue to red: $r = 6, 12, 16, 20, 30, 40, 70, 100$ (a). The base flow acceleration/deceleration is indicated by the black line representing $Q(t)$ schematically (not to scale). (b) A close-up around the numerical abscissa of the corresponding OS spectrum (dashed line) and only the curves for $r=70,100$ are visible.

We observe that the maximum instantaneous growth in the subspace continues to increase as $r$ is increased, only converging slowly as we reach $r=100$. It is known for plane Poiseuille flow that an accurate prediction of transient growth using an eigenmode expansion requires a large number of modes up to the S-branch, although these are highly damped (Reddy & Henningson 1993). The present results indicate that this is true also when we look for transient growth potential within the instantaneously most unstable subspace in pulsating Poiseuille flow, i.e. that growth and decay due to the base flow pulsations do not exploit non-normality but instead coexist with it. This conclusion is in line with the observation that the maximum instantaneous growth rate in the pulsating flow case for $r=100$ is oscillating symmetrically within less than 1 % of the steady value (figure 10b). To cross-check, the maximum non-normal growth rate was computed for the steady case using 100 modes, matching with the value computed from the OS operator to within 0.17 % (not shown), proving that this subspace size is sufficient to capture almost the full non-normal growth potential. The results show that the small fluctuation of the numerical abscissa relative to the steady value over a pulsation cycle is nearly in phase with the acceleration/deceleration of the base flow measured by the sign of ${\textrm {d}}Q/{\textrm {d}}t$. Nevertheless, further studies are necessary to confirm that the growth indeed saturates and to analyse the behaviour for other frequencies and pulsation amplitudes.

The evidence shows that the pulsations mainly affect certain modal growth mechanisms (i.e. the exponential growth rates of the eigenvalues on the A-branch) without noticeably influencing the overall potential for non-normal growth compared with the steady case. To relate this result to earlier studies of optimal initial conditions for Poiseuille flow, these results practically mean that the known transient growth mechanisms in Poiseuille flow (e.g. the Orr mechanism for 2-D channel flow) are unaltered by the pulsations, i.e. the initial growth curves for these perturbations will look similar in both pulsating and steady cases. This does not mean, that there cannot be other transient phenomena that crucially rely on pulsations to emerge but these must be very different from the steady 2-D transient mechanisms, probably act over more than a pulsation period and are heavily restricted in terms of time scales if they are to be relevant compared with modal growth.

6.1.6. Subharmonic perturbation cycles

As we have seen in the previous sections, the variation of the instantaneous eigenvalues of the reduced operator $L_r$ in the limit cycle are characterised either by convergence to the steady values in the case of centre modes or, in the case of eigenvalues from the A-branch, an alignment with the base pulsation frequency. While this is true for most eigenvalue traces, a third alternative exists for the evolution of the eigenvalues. It was found that when the subspace size is sufficiently large, two or more eigenvalue traces merge to form interlaced subharmonic eigenvalue traces. This phenomenon is exemplified in figure 11, showing the eigenvalue traces of $L_r$ over three consecutive pulsation cycles for the same parameters as in figure 10(a) with $r=40$ modes. Note that we use the term ‘subharmonic’ as a general term denoting periodicity with integer multiples of the base period $T_0$.

Figure 11. Traces of selected eigenvalues of $L_r$ corresponding to the A-branch of the spectrum for $\mbox {Re}=7500$, $\alpha =1$, $\beta =0$, $Wo=25$ and $\tilde {Q} = 0.2$ with $r=40$ in the limit cycle. Evolution of the growth rate in time for three consecutive pulsation cycles (a) and orbits in the complex plane over a single pulsation cycle (b). The three eigenvalues that are subharmonic are highlighted in both panels showing their period is $T = 3T_0$. All other eigenvalue traces, including the centre modes (not shown) have the same period as the base flow ($T_0$). The black dots in (b) show the corresponding OS spectrum for reference.

In this situation, the orbit of the dominant eigenvalue has merged with two other orbits and the three eigenvalues follow the same path staggered in time such that the instantaneous spectrum of the reduced operator is periodic with the base frequency. Note that the changes in the periodicity of the individual eigenvalue traces does not affect the intercyclic growth rates nor the instantaneous stability characteristics since the subharmonic eigenvalues simply exchange places from one pulsation cycle to the next. It is interesting to note that period tripling, as observed in this example, is not the only subharmonic that can appear and other combinations have been found for different parameter values, e.g. eigenvalue traces of period $2T_0$. It is likely that such subharmonics are the norm, especially for higher pulsation amplitudes that lead to larger eigenvalue excursions overall.

A comparison of the eigenvalue traces in figures 4(a,b) and 11(a) indicates that the reason the subharmonic cycles are absent for smaller subspace sizes ($r$) is the intense damping of the modes during a large part of the period. If a particular eigenvalue is too damped relative to the others it drops out of the subspace spanned by the OTD basis, which invariably follows the most unstable directions. In fact, the comparatively large fluctuations of the eigenvalue when it is most damped compared with the rest of the cycle suggest that there are more eigenvalue crossings with more damped modes. This could mean that an even larger subspace size is required to trace the full path of the dominant eigenvalue that might, in fact, cover more than just three periods when completely mapped.

A detailed analysis of the subharmonic orbits is outside of the scope of this work, not least because the accurate computation of highly damped eigenvalues leads to unwieldy requirements on wall-normal resolution. A more promising approach for the study of the subharmonic orbits is a fully spectral approach. These analyses can then also accurately trace the damped part of the eigenvalue orbits that are possibly underresolved in this work. While the resolution affects the details of the eigenvalue orbits, the existence of subharmonic perturbation cycles is a persistent feature of the flow case since the eigenvalue crossings that herald them lie in the upper part of the spectrum that is well resolved. In fact, the first 6 modes (i.e. $r\leq 12$) are accurately computed, as can be seen from the match between the growth rates of the leading centre modes and the OS spectrum in figures 4(a,b) and 11(a).

7. Discussion and conclusion

The OTD modes are a recent method for constructing an orthonormal basis of a finite-dimensional subspace spanning and optimally following the instantaneously dominant directions of the tangent space in time along a reference trajectory proposed by Babaee & Sapsis (2016). Apart from the modal stability characteristics, they also capture non-modal growth within the most unstable subspace around a given trajectory. In this work, the potential of the OTD framework for the linear stability analysis of time-dependent flows is explored.

The OTD framework has been implemented in the spectral element code Nek5000 and applied to pulsating Poiseuille flow, allowing for a rigorous validation against the large amount of data that are available for this canonical flow configuration and proof of the accuracy of the implementation. In particular, the results for the leading Floquet exponents from Pier & Schmid (2017) were accurately reproduced and new results obtained for a larger range of pulsation amplitudes than previously reported. Using the framework, the instantaneous variation of the dominant OTD mode of linearly (Floquet) unstable two-dimensional pulsating Poiseuille flow was tracked in the limit cycle and related to the well-known OS modes. This procedure allows for clear characterisation of the stability mechanisms and their interaction that lead to the non-uniform effect of pulsations on the stability of pulsating Poiseuille flow.

In the following, we discuss the main findings of this work. The first point focuses on general considerations concerning the OTD framework and guidelines for its application, while the remaining points pertain to the linear stability characteristics of pulsating Poiseuille flow obtained using the OTD modes.

(i) The OTD modes will recover the $r$-dimensional most unstable subspace of the tangent space. Unfortunately, a general rule for the choice of the subspace size $r$ is not available since it is highly dependent on the system and the goal of the stability analysis. The existing guidelines are either based on the unstable part of the full operator and its symmetric part, which is exhaustive but can be prohibitively large in many systems, or on heuristics. Once $r$ is chosen, care must be taken in the analysis of the OTD modes since the transient phase between initialisation and convergence of the OTD modes to the instantaneously most unstable directions of the tangent space lacks a precise physical meaning. Although the alignment is exponentially fast, the exact time it takes in practice is highly dependent on the system as well as the initial conditions. Once converged, however, the OTD modes depend only on the trajectory pointwise and are memoryless (Blanchard & Sapsis 2019b).

The OTD modes capture both modal and non-modal instabilities within the OTD subspace. This does not imply that the full non-normal growth potential of the complete system is captured since non-normal growth mechanisms may, and often do, involve the interaction of one or more very stable directions that do not appear in finite OTD subspaces of practical size.

The OTD framework can be used to compute reduced-order FTLEs along the trajectory and using the Blanchard & Sapsis formulation makes their computation particularly efficient at negligible extra cost (Babaee et al. 2017; Blanchard & Sapsis 2019a). When considering time-periodic flows, it consequently yields the Floquet exponents with the practical advantage over classical Floquet methods that the framework can be applied directly to the unsteady problem requiring very little extra machinery beyond the standard routines for the solution of the base flow and the linearised Navier–Stokes operator. Moreover, the fact that the OTD modes capture the instantaneously dominant directions in phase space opens many opportunities for reduced-order modelling and control.

When applied to steady base flows, the OTD modes can be readily analysed and completely understood with the existing framework of linear stability analysis. It is therefore in the realm of unsteady flows, in particular flows exhibiting spatio-temporal variation of the base flow that the strengths of the OTD framework are most salient since these flows are inaccessible to global mode analysis. Examples of such flow cases include the analysis of streak instability in bypass transition (Andersson et al. 2001), chaotic edge states (Duguet, Willis & Kerswell 2008; Beneitez et al. 2019) or the minimal seed of turbulence (Pringle & Kerswell 2010; Duguet et al. 2013; Vavaliaris, Beneitez & Henningson 2020) or fully non-autonomous flows dictated by time-dependent external forcing or boundary conditions. While the OTD modes have shown to be successful in a statistical framework (Farazmand & Sapsis 2016) or for flow control (Blanchard et al. 2018), most applications only use them as a stable numerical tool since the physical interpretation of a particular instance of the instantaneous OTD modes on top of a fully transient, nonlinear base flow is far from trivial but might be the key to the analysis of non-autonomous flow cases.

(ii) The OTD framework applied to pulsating Poiseuille flow in the periodic regime generates modes that correspond to the pulsating equivalent of the OS modes. The amount of modulation these modes experience is determined by their spatial localisation in the channel. The wall modes of the A-branch, including the ubiquitous TS waves, are strongly modulated over the pulsation cycle by the changes in the base flow profile close to the oscillating Stokes layer and exhibit considerable fluctuations in both exponential growth rate and phase speed. On the other hand, the growth rates of the modes of the P-branch or centre modes are largely unaffected by the pulsation. It is only for extreme pulsation amplitudes at low frequency that the base flow modulation reaches the centreline, inducing a noticeable intracyclic variation in the leading centre modes.

(iii) The base flow profiles for pulsating Poiseuille flow are highly inflectional and it was found that the amplification phase of the pulsation cycle is well correlated with the existence of two inflection points of the base flow profile outside of the oscillating Stokes layer close to the wall. This result was consistent over the entire range of Womersley number considered. The transition from amplification to damping is generally linked to the collision and disappearance of these inflection points. For low frequencies, where the base flow profile is much less inflectional, the amplification begins already as the first inflection point moves out of the boundary layer. Intuitively, the amplification phase is associated with a W- or $\varLambda$-shaped base flow profile that exhibits two separate shear layers and low flow rates while perturbations decay when the bulk flow is forward.

(iv) The leading OTD mode generally has the structure of a travelling wave close to the wall, similar to the TS waves in steady flow for most of the pulsation cycle. The eigenfunctions are localised in the wall-normal direction at the height of the inflection points in the base flow profile, in particular during the amplification phase. This is in line with the nonlinear simulations of Pier & Schmid (2017) that observe the modulated exponential growth of similar instability waves. As the pulsation amplitude is increased, the vortices are heavily sheared by the mean flow but generally maintain their spatial structure and stay localised close to the wall indicating that, while modulated by the pulsation, the basic mechanism governing the stability is preserved over most of the parameter range. This behaviour prevails also in regions of the parameter space that are Floquet stable.

The only exception to this rule is found during the damping phase for low pulsation frequencies and pulsation amplitudes beyond approximately $\tilde {Q}=0.1$. For these configurations, the strongly damped leading mode moves away from the wall entirely and its spatial structure resembles that of decaying Orr waves that lean with the mean shear. The corresponding nonlinear simulations by Pier & Schmid (2017) also exhibit a different behaviour in this parameter range, termed the ‘ballistic regime’, in which the perturbations completely vanish over part of the pulsation cycle. Considering that characteristic TS-like instabilities were not found in the six-dimensional most unstable subspace during this phase, indicating that they are very damped, it is likely that the change of linear instability mechanism is intimately linked to the appearance of the ballistic regime. This view is supported by the fact that, for the parameter choices investigated, the modulated TS waves (as well as similar wall structures from the A branch) are the only modes that are unstable during the pulsation cycle so that, during their decay phase, there are no linear growth mechanisms available to the flow. At the same time, TS-like structures similar to the saturated nonlinear waves are heavily damped over a long period of time due to the low pulsation frequency, which could explain their complete disappearance and the fact that no other perturbations are excited during part of the oscillation cycle. Direct numerical simulations of pulsating pipe flow by Xu & Avila (2018) show full relaminarisation of turbulent puffs at low pulsation amplitudes that may hinge on a similar mechanism.

(v) By increasing the size of the OTD subspace considerably to $r=100$ for a specific case in order to capture the full non-normal growth potential, it was shown that the pulsations have little effect on the non-normality of the operator compared with the steady case. The intracyclic variations of the numerical abscissa for a moderate pulsation amplitude ($\tilde {Q}=0.2$) were found to lie within 1 % of the steady value. This result indicates that the changes in the instability characteristics due to pulsation do not rely on non-normality but are instead due to the modulation of modal mechanisms that coexist independently. This result would suggest that, while large transient growth responses can be obtained from an optimal initial condition for pulsating flows as shown by Xu et al. (2021) and Tsigklifis & Lucey (2017) in similar geometries, it is likely difficult for the transient mechanisms to fuel their growth efficiently using the base flow pulsations.

(vi) For small and moderate subspace sizes, all OTD modes are periodic with the base flow pulsation frequency in the limit cycle. When the subspace is large enough, the orbit of the dominant eigenvalue of the reduced operator merges with more damped orbits to create subharmonic oscillations that have not been documented for pulsating Poiseuille flow before. They indicate that the perturbations have a more intricate behaviour than previously assumed and that a restriction to perturbations at the base frequency will hide relevant subharmonic structures. The appearance of the subharmonic orbits is likely linked to the transient proximity of the eigenvalues in the complex plane during the pulsation cycle, leading them to switch place. On the other hand, a bifurcation cascade that is one of the routes to chaos (Feigenbaum 1980) is unlikely their cause since the subharmonic orbits do not exhibit the characteristic period doubling and do not seem to emerge at an accelerating rate. Further research is necessary to analyse and understand the appearance and behaviour of the subharmonic perturbation cycles. Given that the dominant eigenvalues are very damped during part of the subharmonic cycle, such an endeavour will need to have a high wall-normal resolution and should therefore be carried out with a spectral approach.

It may be noted that, irrespective of whether a spectral or pseudo-spectral method is employed for the study of perturbations on a periodic base flow, the OTD framework is ideally suited for the task because it allows the perturbations to be followed even if they are highly damped without assumptions about their temporal periodicity.

To highlight the physical insights gained from applying the OTD framework to the case of pulsating Poiseuille flow, we summarise the main aspects of the discussion.

The main driver of instability in pulsating Poiseuille flow just beyond to the critical Reynolds number for linear instability is a modulated version of the TS wave known from steady shear flows that is prominent across the parameter space considered in this work. This result is in qualitative agreement with DNS results in the literature. Moreover, at low pulsation frequencies, the damping phase is associated with a different spatial structure reminiscent of Orr waves leaning with the mean shear and the disappearance of perturbations akin to TS waves which may explain the emergence of the ballistic regime in nonlinear simulations by Pier & Schmid (2017). The amplification phase is marked by the co-existence of two inflection points in the base flow profile outside of the boundary layer. In all cases damping begins close to the collision and disappearance of the inflection points closest to the centreline. For low pulsation frequencies, where the base flow profile only has a single inflection point for most of the cycle, the amplification begins when it moves out of the boundary layer. The dominant perturbations are shown to localise at the inflection points, highlighting their importance even in transient flow stability at high pulsation frequencies. The computations spanning almost the full non-normal growth potential of the linear operator indicate that pulsations have little effect on known 2-D transient growth mechanisms such as the Orr mechanism. Finally, as the size of the OTD subspace was increased considerably, some of the eigenvalue orbits were seen to merge into subharmonic cycles. This interesting phenomenon is not yet fully understood but is likely crucial for a complete picture of the transient stability of pulsating Poiseuille flow.