Nonlinear and Non-Modal Stability of Structures Evolving in Shear Flows Conor A. Daly Queens’ College, Cambridge. A Dissertation submitted for the degree of Doctor of Philosophy at the University of Cambridge September 2013 Declaration This thesis describes research carried out in the Department of Applied Mathematics and Theoretical Physics of Cambridge University. The dissertation is the result of my own work and includes nothing which is the outcome of work done in collaboration, except where specifically indicated in the text. The numerical calculations presented were com- puted using my own MATLAB code, with the exception of the nonlinear simulations of the Navier-Stokes equations and computation of nonlinear solutions in §3. These calcula- tions were made using Navier-Stokes solvers SIMSON and Channelflow with permission from Dr P. Schlatter and Dr T. Schneider respectively. No part of the work contained herein has been submitted to any other university or place of learning for any degree, diploma or other qualification. Conor A. Daly Cambridge, September 2013 Acknowledgements First and foremost I would like to thank my supervisor Prof. Nigel Peake for his perpet- ual support and limitless patience throughout my time in Cambridge. Nigel’s willingness to listen and readiness to provide direction have been invaluable in shaping my thesis, and guiding my entry into the world of research. This, alongside regular discussions on the current state of LFC and the occasional whipping on the squash court, have allowed the past four years to pass very quickly. A great deal of thanks must also go to my office- partner Dr. Ed Brambley, without whose generous insights, my research would not be where it is today. I have been extremely fortunate to share a space with someone whose passions for mathematics, computing and freshly ground Arabica are congruent with my own. I would also like to thank past and present members of the Waves Group: Rosie, He´le`ne, Justin and Lorna for their feedback on my work and for generously passing on their expertise in LaTeX. I would like to further thank all the friends I have made in DAMTP and the CMS, who have made my time here enjoyable and unforgettable. It has been a great pleasure to spend three months of my postgraduate study working at KTH in the beautiful city of Stockholm. I would like to thank Prof. Henrik Alfredsson, whose presentation during his visit to DAMTP directly inspired two of the chapters in my thesis; Dr. Philipp Schlatter, who taught me a great deal of computational fluid dy- namics and under whose influence I took strides forward in my research; Dr. Alexandre Suryadi, with whom I had many insightful discussions on experimental fluid dynamics; and Prof. Dan Henningson, whose passion for science was simply inspirational and whose garden house provided much more than a roof over my head! I would also like to thank the many friends I made at the Mechanics department in KTH, who made me feel at home during my visit. I would like to thank Dr. Tobias Schneider for his very generous invitation to the Max Planck Institute, which allowed me to develop my work on nonlinear structures. My short time in the charming town of Go¨ttingen was extremely fruitful, and for this I owe much to the boundless enthusiasm of Dr. Schneider and the computational prowess of Dr. Hecke Schrobsdorff. I would not be where I am today were it not for the support and encouragement of my parents and family. Mum, who has always pushed me to achieve without pressuring me to succeed, and Dad, who eventually learned not to ask “Is your course finished yet?”, but is unrelenting in his ability to administer the occasional whipping on the squash court. I would like to thank the many friends I have made here in Cambridge and in Queens’, all of whom have made my years here very special. There are too many to mention everyone by name, but huge thanks go to my housemates Marion, Brendan, Liz and Stuart, and to Kavita who has given me encouragement and support beyond what I deserve. Finally I would like to thank the EPSRC, DAMTP and St. Malachy’s College for pro- viding the finanical backing which has made my research possible. Abstract This thesis explores a range of stability techniques and their application to structures in various constant density fluid flows. In particular, the stability of nonlinear struc- tures which develop in rotating plane Couette flow is analyzed using Floquet theory, which allows the global stability of an important secondary nonlinear structure called a Taylor vortex to be determined. From this the distinct tertiary states which emerge as Taylor vortices break down are characterized and their bifurcation behaviour is studied. Non-modal stability analyses are conducted in rotating plane Couette flow and annular Poiseuille-Couette flow. In each flow the growth mechanisms and the form of the pertur- bations responsible for the maximum linear energy amplification are discussed. Finally, the non-modal behaviour of the Papkovitch-Fadle operator is treated with reference to spatially developing disturbances in Stokes channel flow. The mechanisms and the rates of convergence of the linear spatial energy amplification are investigated and contrasted with temporal energy amplification. Contents 1 Introduction 5 2 Non-modal growth in annular Poiseuille-Couette flow 11 2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1 Alternative non-dimensionalizations . . . . . . . . . . . . . . . . . 14 2.1.2 Thread-annular flow . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2 Linear stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Non-modal growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Optimal growth in APCF . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.1 Optimal perturbations and growth mechanism . . . . . . . . . . . 31 2.4.2 Axisymmetric growth . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4.3 Radius ratio dependence of G . . . . . . . . . . . . . . . . . . . . 38 2.5 Comparison to experimental data . . . . . . . . . . . . . . . . . . . . . . 40 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3 Stability and bifurcations in rotating plane Couette flow 45 3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Linear stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.1 Flowfield visualization . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Nonlinear states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.4 Floquet stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.1 Steady states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.2 Time-periodic states . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4.3 Global stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.5 Secondary states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 1 2 Contents 3.5.1 Taylor vortex flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.5.2 Second Taylor vortex flow . . . . . . . . . . . . . . . . . . . . . . 62 3.5.3 Oblique vortex flow . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.5.4 Summary of secondary states . . . . . . . . . . . . . . . . . . . . 67 3.6 Stability of Taylor vortex flow . . . . . . . . . . . . . . . . . . . . . . . . 68 3.6.1 Streamwise instability maps . . . . . . . . . . . . . . . . . . . . . 70 3.6.2 Eckhaus instability maps . . . . . . . . . . . . . . . . . . . . . . . 79 3.6.3 Summary of TVF stability properties . . . . . . . . . . . . . . . . 82 3.7 Bifurcations of tertiary states . . . . . . . . . . . . . . . . . . . . . . . . 83 3.7.1 Wavy vortex flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.7.2 Twist and wavy twist vortex flow . . . . . . . . . . . . . . . . . . 84 3.7.3 Oscillatory wavy vortex flow . . . . . . . . . . . . . . . . . . . . . 85 3.7.4 Skewed vortex flow . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.7.5 Bifurcation of twist, wavy twist and wavy vortex flows . . . . . . 88 3.7.6 Bifurcation of wavy and oscillatory wavy vortex flows . . . . . . . 94 3.7.7 Bifurcation of skewed vortex flow . . . . . . . . . . . . . . . . . . 98 3.8 Transition modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.8.1 Numerical procedure . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.8.2 Model prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.8.3 Time-periodic flow: Ω = 2.5 . . . . . . . . . . . . . . . . . . . . . 106 3.8.4 Wavy vortex flow: Ω = 10 . . . . . . . . . . . . . . . . . . . . . . 109 3.8.5 Taylor vortex flow: Ω = 20 . . . . . . . . . . . . . . . . . . . . . . 111 3.8.6 Taylor vortex flow: Ω = 50 . . . . . . . . . . . . . . . . . . . . . . 112 3.8.7 Wavy twist vortex flow: Ω = 90 . . . . . . . . . . . . . . . . . . . 115 3.8.8 Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4 Non-modal growth in rotating plane Couette flow 123 4.1 Non-modal growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.1.1 Maximum optimal growth . . . . . . . . . . . . . . . . . . . . . . 126 4.2 Optimal growth in RPCF . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.2.1 Cyclonic rotation and the non-rotating case . . . . . . . . . . . . 130 4.2.2 Subcritical anti-cyclonic rotation . . . . . . . . . . . . . . . . . . 137 4.2.3 Rotational independence of the Orr mechanism . . . . . . . . . . 140 CONTENTS 3 4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5 Nonnormality and optimal growth of the Papkovitch-Fadle operator 145 5.1 Governing equations and modal solution . . . . . . . . . . . . . . . . . . 145 5.2 Nonnormality and optimal energy growth . . . . . . . . . . . . . . . . . . 148 5.2.1 Optimal growth and convergence . . . . . . . . . . . . . . . . . . 150 5.2.2 Optimal growth with bounded wall shear . . . . . . . . . . . . . . 157 5.3 Short wavelength asymptotics . . . . . . . . . . . . . . . . . . . . . . . . 159 5.3.1 Normalization factor . . . . . . . . . . . . . . . . . . . . . . . . . 161 5.3.2 Interaction of short wavelength modes . . . . . . . . . . . . . . . 162 5.3.3 Interaction of long and short wavelength modes . . . . . . . . . . 163 5.3.4 Short wavelength growth mechanism . . . . . . . . . . . . . . . . 164 5.4 Spatial Poiseuille flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 5.4.1 Spatial spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 5.4.2 Optimal growth and convergence . . . . . . . . . . . . . . . . . . 170 5.4.3 Optimal disturbances . . . . . . . . . . . . . . . . . . . . . . . . . 174 5.5 Temporal plane Couette flow . . . . . . . . . . . . . . . . . . . . . . . . . 176 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 6 Conclusions & further work 185 6.1 Annular Poiseuille-Couette flow . . . . . . . . . . . . . . . . . . . . . . . 185 6.2 Rotating plane Couette flow . . . . . . . . . . . . . . . . . . . . . . . . . 186 6.3 The Papkovitch-Fadle operator and spatial non-modal growth . . . . . . 189 Bibliography 191 4 Contents Chapter 1 Introduction Hydrodynamic stability and the transition to turbulence of fluid flow have been un- der investigation and scientific scrutiny since the late 19th century, when the studies of Rayleigh (1880) and Reynolds (1883) provided early insights into the mathematics and physics involved. With the equations governing the dynamics of a viscous fluid already established by Navier (1822) and Stokes (1845), Rayleigh’s contribution, though it is one of many in the field, was to derive an equation governing the evolution of infinitesimal, inviscid waves which perturb a finite mean flow. The ideas underpinning the formulation of this equation, known as the Rayleigh equation, are a major influence on our modern understanding of stability in its mathematical sense. Reynolds’ experimental investi- gations on transition to turbulence in pipe flow heralded the discovery that transition could be characterized in terms of a non-dimensional quantity which takes account of the spatio-temporal scales of the system and the viscosity of the fluid. The Reynolds number, defined Re = Uhν for flow in a channel where U is the mean flow velocity, h is the channel half-width and ν is the kinematic viscosity, has become one of the most important numbers in fluid dynamics. Its discovery has precipitated the identification of a range of non-dimensional numbers which describe various aspects of fluid flow. The inconsistency between experimental observations of the onset of turbulence and theoret- ical predictions of when transition should occur remains an outstanding problem in 21st century physics, and motivation for current researchers around the world. The possibility of transiently growing structures in shear flows has long been recognized in fluid mechanics. Among the first to predict the existence of such disturbances was Orr (1907) who, on analyzing the dynamics of a two-dimensional inviscid fluid and under- scoring the existence of a transient instability mechanism, wrote, “It accordingly appears that, in this simple case, although the disturbance, if sufficiently small, must ultimately decrease indefinitely, yet, before doing so, it may be very much increased”, Orr (1907, p. 32). Transient, non-modal, energy growth of linear perturbations has now become a 5 6 1. Introduction cornerstone of hydrodynamic stability analysis, alongside classical eigenvalue, or modal, analysis. Researchers such as Ellingsen & Palm (1974) and Landahl (1980) added to the literature on inviscid, algebraic (as opposed to expontential) growth mechanisms extant in linearly stable fluids. Gustavsson (1991), Butler & Farrell (1992), Reddy & Henningson (1993) and Schmid & Henningson (1994) showed that transient energy growth of linear disturbances over many orders of magnitude is possible in viscous flows, with growth driven by the inviscid, algebraic instabilities and subsequently cut off by viscous dissipa- tion after long times. It was recognized that this transient energy growth is a consequence of the underlying nonnormality of the linearized Navier-Stokes operator (Trefethen et al. (1993), Trefethen (1997), Trefethen & Embree (2005)), which renders eigenmodes of the system non-orthogonal with respect to an appropriate norm. Hence we term the study of transient linear energy growth non-modal analysis, as a combination of many modes are a necessary requirement for growth. Non-modal growth has been postulated to play a role in the onset of turbulent dynamics in subcritical shear flows, where experimental ob- servations of turbulence are very much at odds with the predicitons of the exponentially growing linear modes of classical linear stability theory. From linear stability theory, we have that plane Poiseuille flow has critical Reynolds number Recr = 5772 (Orszag (1971)) and plane Couette flow is stable for all Reynolds numbers (Romanov (1973)). Numerical calculations suggest that Hagen-Poiseuille flow is linearly stable, though a formal proof of this continues to elude mathematicians. However, in experimental and numerical in- vestigations, transition is found to occur as low as Re = 1000 for plane Poiseuille flow, Re = 360 in plane Couette flow and Re = 1020 in Hagen-Poiseuille flow. Recent studies have proposed fully nonlinear mechanisms governing subcritical shear flow transition. These include the self-sustaining process proposed by Waleffe (1997, 1998) and its asymptotic analogue of vortex-wave interaction theory featured in the works of Hall & Smith (1991); Hall & Sherwin (2010), Deguchi et al. (2013) and Blackman et al. (2013). The self-sustaining process and vortex-wave interaction theory are reliant on streak instabilities; a streaky structure forms and loses stability to an inviscid wave; the neutral inviscid wave drives streamwise oriented vortices; and the streamwise vortices in turn drive the streaks. Non-modal growth, in the form of the lift-up mechanism (Lan- dahl (1980)), has been proposed as the cause for the streamwise vortices being driven into streaks. In another approach to the problem of transition, optimization over the fully nonlinear Navier-Stokes equations for the minimal seed required for transition away from laminar flow has been undertaken by Pringle et al. (2012), Rabin et al. (2012) and Duguet et al. (2013); non-modal linear energy amplification plays a role in the nonlinear optimals discovered by these authors. The mechanism proposed by Orr (1907) drives the early growth states of the nonlinear optimal and the lift-up mechanism is once again behind the emergence of streaks. In this light, we consider a non-modal stability analysis 7to be an important contribution to the understanding of stability and transition in fluid flows, and this is the motivation for the non-modal stability analyses in §§2, 4 and 5 of this thesis. In §2 we study the non-modal stability of annular Poiseuille-Couette flow (APCF), a cylindrical flow where the fluid is driven by the combination of an axial pressure-gradient and a sliding inner cylinder. APCF has received attention recently as it incorporates thread-annular flow which is relevant to the medical procedure of thread-injection, a minimally-invasive surgical technique whereby a thread is inserted into a body using a spool and fluid-filled syringe apparatus (see Frei et al. (2000) for more details). Walton (2003, 2004, 2005) built upon the earlier linear stability analyses of Sadeghi & Higgins (1991) and Gittler (1993) to study the nonlinear stability and axisymmetric linear stabil- ity of APCF. Sadeghi & Higgins (1991) computed non-axisymmetric neutral curves and, following the approach of Cowley & Smith (1985) for planar Poiseuille-Couette flow, de- rived asymptotic cut-off inner cylinder velocities beyond which no linear instability is supported. Gittler (1993) showed that disconnected axisymmetric neutral curves exist for APCF, a finding later confirmed by Walton (2004) in a comprehensive axisymmetric stability analysis. Walton (2003, 2005) also found non-axisymmetric, large Re nonlin- ear waves to be supported over a wide parameter range and explored their relevance to thread-annular flow. Experimental studies of APCF are found in papers by Frei et al. (2000), who focus on thread-annulus applications, and a short paper by Shands et al. (1980), where transition to turbulence and fully turbulent flow are studied. In particular, Shands et al. (1980) show that transition occurs far below the linear stability threshold and, furthermore, that an increase of inner cylinder velocity de-stabilizes the flow which is at odds with the trend predicted by linear stability theory. Most recently, Wong & Wal- ton (2012) computed finite-amplitude axisymmetric travelling wave solutions and showed that they can be continued below the linear neutral curve to subcritical parameter space. Non-modal growth analysis has been applied to related flows, such as Hagen-Poiseuille flow (Schmid & Henningson (1994)) and sliding Couette flow (Liu & Liu (2012)), but not APCF. Therefore, our aim is to complement the existing works concerning transition in APCF with a non-modal stability analysis. The seminal work of Taylor (1923) mixed experimental and theoretical methods to de- scribe the emergence of vortex structures in a differentially rotated, concentric-cylinder apparatus. These vortex structures are now called Taylor vortices and are considered key structures in transition of cylindrical and curved fluid flow (see Koschmieder (1993) for a full review). The Taylor vortices were shown to develop due to a linear inviscid instabil- ity of the base flow, putting experiment and theory in good agreement with one another 8 1. Introduction and providing some hope that the phenomena of transition and turbulence may be the- oretically understood. Taylor’s experiments sparked a flurry of research activity, which included: theoretical studies of Taylor vortex instability such as Davey et al. (1968) and Eagles (1971), who used weakly nonlinear analysis to determine the instabilities which affect Taylor vortices; experimental papers such as Coles (1965), Andereck et al. (1986) and Hegseth et al. (1996), who mapped in parameter space the different flow regimes observed in the experiments; and numerical stability analyses wherein finite-amplitude Taylor vortices are calculated numerically in addition to the higher-order structures they bifurcate towards, such as Nagata (1988), Weisshaar et al. (1991) and Antonijoan & Sa´nchez (2000). Recent experimental attention has been focused on the flow of a dif- ferentially rotated fluid through a linear shear layer, known as rotating plane Couette flow (RPCF). Alfredsson & Tillmark (2005), Hiwatashi et al. (2007), Tsukahara et al. (2010) and Suryadi et al. (2013) carried out experimental investigations of transition, with Tsukahara et al. (2010) following Andereck et al. (1986) by making a map in pa- rameter space demarcating the different flow regimes which exist. Though RPCF is a more challenging flow to model experimentally, it is perhaps more amenable to a the- oretical or numerical analysis due to the Cartesian geometry and the ease with which rotation can be added to the governing equations. The equations governing RPCF can be interpreted as a local approximation to Taylor-Couette flow in the co-rotating, narrow- gap limit, and indeed many of the aforementioned theoretical and numerical studies are based in the Cartesian framework. In §3 we attempt to further the theoretical understanding of transition in supercritical RPCF. We begin with a characterization of the linear instabilities which de-stabilize the basic, primary flow and the finite-amplitude structures, often called secondary states in this context, which emerge from each unstable mode. The most important secondary structures are the Taylor vortices discussed in the preceeding paragraph, as the least stable perturbations of the primary flow bifurcate towards them. We perform a thorough investigation of the stability properties of Taylor vortices using Floquet theory to solve for perturbation growth rates. We study the stability of Taylor vortices with a range of spanwise wavenumbers β, to extend previous works such as Nagata (1998) who focused on Taylor vortices in RPCF with the critical primary wavenumber β = 1.558. Our sta- bility analysis is then used to investigate tertiary bifurcations in RPCF, with a study of the structures that Taylor vortices develop into as they lose stability. Many of these structures have an analogue in the cylindrical Taylor-Couette flow, such as wavy vortex flow (Davey et al. (1968), Nagata (1988)) and twist vortices (Weisshaar et al. (1991), An- tonijoan & Sa´nchez (2000)), but our approach uncovers a new streamwise-independent structure and a basic oscillatory state. Our analysis of the bifurcation sequence in RPCF is subsequently used to study the flow regimes observed in the experiments of Tsukahara 9et al. (2010) and Suryadi et al. (2013). We propose a simplified transition model, the aim of which is to provide insight into the complex dynamics of transition, which we use in a comparison between experimental observation and theoretical and numerical modelling. §4 continues our study of RPCF with a non-modal stability analysis on its subcritical parameter subspace. RPCF has two rotational regimes: cyclonic rotation, for which the mean vorticity is aligned with the system rotation and anti-cyclonic rotation, for which the mean vorticity and the system rotation counteract one another. The anti-cyclonic rotation regime has been of particular interest in the astrophysics community, where it is used as a local model for large-scale astrophysical accretion disks. The possibility of hydrodynamic turbulence in accretion disks is a subject of lively debate amongst as- trophysicists. It is postulated that hydrodynamic turbulence could be the driving force behind an inward transport of angular momentum, which is necessary for the disks to accrete matter (Narayan et al. (1987), Longaretti (2002), Tevzadze et al. (2003), Yecko (2004), Shartman et al. (2012)). In this vein, Yecko (2004) made calculations of the non-modal growth for anti-cyclonic rotation. Using the rotation number Ro = 2ΩhU as a non-dimensional measure of rotation, where Ω is the rotation rate of the system, h is the channel half-width and U is the speed of each opposingly moving plate, Yecko found strong growth in the case Ro = 1, with this rotation rate related to a constant angular momentum rotating disk. This growth is driven by the anti lift-up effect (Antkowiak & Brancher (2007)), a mechanism with opposing characteristics to the lift-up mechanism of non-rotating flow. Yecko reported weaker two-dimensional growth for Ro = 43 , which is a local model for a Keplerian accretion disk, with this growth attributed to the mechanism originally proposed by Orr. Yecko subsequently concluded that anti-cyclonic hydrody- namic transition would be possible for high Reynolds numbers. We do not address the specific route in which a disturbance reaches turbulence in our work, rather, we identify the dominant linear energy growth mechanisms which we suspect will be relevant to the transition process. Finally, in §5 we study the nonnormality of the Papkovitch-Fadle operator. The Papkovitch-Fadle operator is derived from the mathematical system which independently describes physical problems in both fluid mechanics and solid mechanics. In fluid mechanics, the problem corresponds to the Stokes flow of a viscous fluid within a semi-infinite two-dimensional strip. In solid mechanics, modelling the elastic deformation of a clamped semi-infinite plate requires the solution of the Papkovitch-Fadle operator. The operator is named after Papkovitch (1940) and Fadle (1941) who independently worked on the problem of elastic deformation in a finite strip, or rectangle, and in so doing gave expressions for the eigenfunctions relating to the semi-infinite domain. Much 10 1. Introduction of the previous literature on the problem is concerned with the suitability of expressing arbitrary edge data as a sum of these eigenfunctions. Smith (1952) examined the solid mechanics problem and developed an algorithm to determine the bending displacement of a semi-infinite elastic strip, under arbitrary displacements and stresses at the free edge. Smith’s solution takes the form of a series of Papkovitch-Fadle eigenfunctions and conditions are established to guarantee convergence of that series. Noting that these conditions are too restrictive for most applications, Joseph (1977); Joseph & Sturges (1978), Gregory (1979, 1980a,b) and Spence (1982) showed that convergence could be guaranteed beyond Smith’s conditions, and used their results to solve a wider range of Stokes flow and plate elasticity problems than previously possible. Trefethen (1997); Trefethen & Embree (2005) highlighted the nonnormality of the Papkovitch-Fadle operator and computed its pseudospectra and plots of the optimal amplification of arbitrary edge conditions, along the semi-infinite dimension. The opti- mal disturbances were recognized to be Moffatt eddies (Moffatt (1963)), which are known to develop in a wall bounded corner of a viscous fluid. Optimal growth of spatially devel- oping disturbances has also been treated in the context of boundary layers by Andersson et al. (1999), Luchini (2000) and Tempelmann et al. (2010) where it has been used as a key component in prediction models for transition to turbulence. A spatial optimal growth framework was also developed by Reshotko & Tumin (2001) in Hagen-Poiseuille flow, and for planar, wall-bounded parallel flows by Schmid & Henningson (2001). These works on spatial non-modal growth complement the much larger literature on non-modal growth in temporal flows, and together they form the basis of our understanding of the spatio-temporal development of structures within a fluid. In §5 we return to the question of the maximum spatial amplification that can be achieved by an arbitrary edge condi- tion expressed as a finite series of Papkovitch-Fadle eigenfunctions. Framing the problem in a Stokes flow context, we investigate the convergence of the optimal edge condition with modal truncation. We then assess the extent to which our findings for Stokes flow are relevant to flows with spatial convection by studying two-dimensional spatial plane Poiseuille flow. To conclude, we compare our findings for spatial flows with an asymp- totic analysis of the convergence of optimals in temporal plane Couette flow. Chapter 2 Non-modal growth in annular Poiseuille-Couette flow In this chapter we consider the temporal transient energy amplification of non-modal disturbances in annular Poiseuille-Couette flow (APCF); the flow between concentric cylinders with an axially translating inner cylinder and a pressure gradient driving the fluid across the annular gap. The outline of the chapter is as follows: in §2.1 we introduce the governing equations and derive the axially and azimuthally independent laminar base profile. In §2.2 we recap previous authors’ linear stability results and in so doing verify our own calculations. We further maximize over wavenumbers to produce neutral curves in (Re, V ) parameter space, with Re the Reynolds number based on the pressure gradi- ent and V the inner cylinder velocity. In §2.3 we formulate the optimal energy growth problem and define the maxima by which we analyse the transient energy growth. In §2.4 we present the results of our non-modal growth calculations and assess the growth mechanisms driving the linear energy amplification. Finally, in §2.5 we compare our findings with the experimental study of Shands et al. (1980). 2.1 Governing equations We consider annular Poiseuille-Couette flow (APCF); the flow of an incompressible, New- tonian fluid bounded between two rigid, co-axial cylinders with impermeable walls. A constant axial pressure gradient drives the fluid in the gap between the cylinders and, in addition, the inner cylinder is free to move with a prescribed velocity, in the axial direction. We work in cylindrical polar coordinates (x˜, r˜, θ˜), aligned along the common axis of the cylinders. The fluid velocity u˜ and pressure p˜ are governed by the cylindrical 11 12 2. Non-modal growth in annular Poiseuille-Couette flow Navier-Stokes equations ∂u˜ ∂t˜ + u˜ · ∇˜u˜ = − 1 ρ ∇˜p˜+ ν∇˜2u˜, (2.1a) ∇˜ · u˜ = 0, (2.1b) where ρ and ν are the constant fluid density and kinematic viscosity, respectively. Given an inner cylinder radius A and outer cylinder radius B, the no-slip boundary condition is applied on the cylinder walls r = A,B. For axial pressure drop ∆p over pipe length L and an inner cylinder velocity of V˜ , we have the following axisymmetric and axially independent exact solution of (2.1), which serves as a laminar flow profile U˜(r˜) = U˜p(r˜) + U˜c(r˜), (2.2a) with U˜p(r˜) = ∆p 4µL log BA ( r˜2 log B A − (B2 − A2) log r˜ +B2 logA− A2 logB ) , (2.2b) U˜c(r˜) = V˜ log Br˜ log BA , (2.2c) where µ = ρν is the dynamic viscosity and A < r˜ < B. U˜p(r˜) and U˜c(r˜) denote the Poiseuille and Couette components, respectively. The geometry is depicted in Figure 2.1. To non-dimensionalize the flow we use the radial gap length scale B−A, the fluid density ρ and the mean Poiseuille velocity U . For cross-sectional area Ac, U is defined as U = 1 Ac ∫ Ac U˜p dS = 2 B2 − A2 ∫ B A U˜p(r) rdr, (2.3a) = ∆p 8µL log BA ( B2 − A2 − (B2 + A2) log B A ) . (2.3b) The flow can then be described in terms of non-dimensional quantities; the radius ratio η of the cylinders, the Reynolds number Re based on the pressure gradient and the non-dimensionalized velocity V of the translating inner cylinder η = A B , Re = U(B − A) ν , V = V˜ U . (2.4) The non-dimensionalized laminar profile is U(r) = 2 (1− r2(1− η)2) log η − (1− η2) log(r(1− η)) 1− η2 + (1 + η2) log η + V log(r(1− η)) log η , (2.5) 2.1 Governing equations 13 for η1−η < r < 1 1−η and with 0 < η < 1. In Figure 2.2 we plot the laminar profile at various values of η and V . Assuming an infinitesimal perturbation of the form [u¯, v¯, w¯, p¯] = [u, v, w, p](r)ei(kx+mθ)+ωt, (2.6) with k ∈ R, m ∈ Z and ω ∈ C, we linearize the Navier-Stokes equations about the laminar profile U(r) to get ωu+ ikUu+ vU ′ = −ikp+ 1 Re { u′′ + u′ r − ( m2 r2 + k2 ) u } , (2.7a) ωv + ikUv = −p′ + 1 Re { v′′ + v′ r − ( 1 +m2 r2 + k2 ) v + 2imw r2 } , (2.7b) ωw + ikUw = − imp r + 1 Re { w′′ + w′ r − ( 1 +m2 r2 + k2 ) w + 2imv r2 } , (2.7c) iku+ v′ + v r + imw r = 0, (2.7d) where ′ = ddr . The no-slip boundary condition holds on each cylinder surface, giving u = v = w = 0 at r = η 1− η , 1 1− η . (2.8) Figure 2.1: Flow geometry with dimensional quantities labelled. 14 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.2: Base flow profiles in r for the radius ratios η = 142 , 0.2, 0.5, 0.8. U(r) with V = −15,−10,−5,−2,−1, 0, 1, 2, 5, 10, 15 is plotted for each radius ratio. The black dashed line in each Figure marks the inner cylinder boundary r = η1−η . 2.1.1 Alternative non-dimensionalizations Our choice of non-dimensionalization is not unique; various alternative characteristic length scales and velocities have been used by previous authors studying APCF. In the experiments of Shands et al. (1980) and Frei et al. (2000), lengths and velocities are non-dimensionalized by the hydraulic diameter, D = 2(B −A), and the total mean flow velocity Û = 1 Ac ∫ Ac U˜p + U˜c dS = U + V , (2.9a) = ∆p 8µL log BA ( B2 − A2 − (B2 + A2) log B A ) + V˜ ( B2 − A2 − 2A2 log BA 2(B2 − A2) log BA ) . (2.9b) 2.1 Governing equations 15 The Reynolds number R̂e and non-dimensional inner cylinder velocity V̂ , defined in terms of D and Û , are R̂e = ÛD ν , V̂ = V˜ Û , (2.10) where ν is the kinematic viscosity of the fluid as before. In order to convert from the experimental quantities to ours, we first note that, V = V˜ ( η2 − 1− 2η2 log η 2(1− η2) log η ) . (2.11) After some straightforward manipulations we derive the following relations for Re and V in terms of R̂e and V̂ Re = R̂e 2 { 1− V̂ ( η2 − 1− 2η2 log η 2(1− η2) log η )} , (2.12a) V = 2V̂ (1− η2) log η 2(1− η2) log η − (η2 − 1− 2η2 log η)V̂ . (2.12b) In the theoretical works of Walton (2003, 2004, 2005) the outer cylinder radius, B, is used as the length scale and velocities are non-dimensionalized using the pressure driven velocity Uw = ∆pB2 4µL . Hence Walton’s parameters, Rew Vw and kw are related to ours through Re = 1− η f(η) Rew, (2.13a) V = f(η)Vw, (2.13b) k = (1− η)kw, (2.13c) where we define f(η) = 2 log η 1− η2 + (1 + η2) log η . (2.13d) Further conversion is needed for comparison with the results of Sadeghi & Higgins (1991), who chose the maximum pressure driven velocity, maxr˜ U˜p, and cylinder half-gap, B−A2 , as their characteristic quantities Re = 2 f(η)g(η) Resh, (2.14a) V = f(η)g(η)Vsh, (2.14b) k = 2ksh, (2.14c) 16 2. Non-modal growth in annular Poiseuille-Couette flow where f(η) is defined as before, and g(η) = η2 + 1− η2 2 log η ( 1− log ( η2 − 1 2η2 log η )) . (2.14d) 2.1.2 Thread-annular flow A pertinent physical application of APCF is in thread-annular flow, which is a model for a medical injection procedure involving insertion of an implant into a body in a minimally invasive manner. In the procedure, a thread is surrounded by a fluid and stored on a spool within a container with a long cylindrical outlet. It is then inserted into the desired location by pressurizing the fluid such that the thread travels along the cylindrical section of the container before entering the body. More details on the technique can be found in Frei et al. (2000). A crude mathematical model for the dynamics of the procedure is a system of pressure-driven fluid constrained to move in the gap between concentric cylinders with a sliding inner cylinder. Thread-annular flow can be considered a specific case of APCF, where the inner cylinder is not independently driven, but free to slide. A maximum inner cylinder velocity can be derived for the thread-annulus problem, by balacing the pressure-driven forces on the inner cylinder with the resistive action of the shear. Assuming the fluid to be in the laminar solution (2.5), the non-dimensionalized resistive shear force Fτ acting on the inner cylinder over a length ` is Fτ = 2piη`Re 1− η ( dU dr )∣ ∣ ∣ ∣ r= η1−η = 2pi`Re { V log η − ( 4η2 log η + 2(1− η2) 1− η2 + (1 + η2) log η )} . (2.15) The non-dimensionalized force on the inner cylinder due to the pressure gradient across the annular gap is Fp = −pi` ( η 1− η )2dP dx = pi`Re ( 8η2 log η 1− η2 + (1 + η2) log η ) (2.16) The total resistive force on the inner cylinder is then Fr = Fp +Fτ . When Fr is zero, the thread reaches its maximum velocity Vmax = 2(1− η2) log η 1− η2 + (1 + η2) log η . (2.17) 2.2 Linear stability 17 2.2 Linear stability By modification of the Rayleigh inflection point criterion to account for cylindrical ge- ometry, we investigate whether APCF is linearly stable to inviscid perturbations. We begin with the Rayleigh equation in cylindrical coordinates (U − c) ( φ r − ( φ′ κr )′) − ( U ′ κr )′ φ = 0, (2.18) where φ = rv, κ = k2 + m 2 r2 and c = −iω k in our notation. Multiplication of (2.18) by the factor ( φ ∗ U−c), where φ ∗ is the complex conjugate of the potential φ, and integrating over the annular gap gives ∫ 1 1−η η 1−η φφ∗ r − ( φ′ κr )′ φ∗dr − ∫ 1 1−η η 1−η ( U ′ κr )′ φφ∗ U − c dr = 0. (2.19) Integration by parts of the integral on the left gives ∫ 1 1−η η 1−η |φ|2 r + |φ′|2 κr dr − ∫ 1 1−η η 1−η ( U ′ κr )′ |φ|2 U − c dr = 0. (2.20) The integral on the left is strictly positive and real. The second integral is complex valued, and with ci denoting the imaginary part of the complex wavespeed c, we express the imaginary part of the second integral as ci ∫ 1 1−η η 1−η ( U ′ κr )′ |φ|2 |U − c|2 dr = 0. (2.21) Since |φ| 2 |U−c|2 is strictly positive, in order for the integral to be satisfied for ci > 0 the factor (U ′ κr ) ′ must change sign somewhere in the interval η1−η ≤ r ≤ 1 1−η . However, for the base flow (2.5), we have ( U ′ κr )′ = 2(2Pm2 −Qk2)r (k2r2 +m2)2 , (2.22) with the constants P and Q given by P = − 2(1− η)2 log η 1− η2 + (1 + η2) log η (2.23) Q = V log η − 2(1− η2) 1− η2 + (1 + η2) log η . (2.24) Provided 2Pm2 − Qk2 6= 0, (U ′ κr ) ′ cannot change sign in the interval η1−η ≤ r ≤ 1 1−η , implying ci = 0 in order for expression (2.21) to be satisfied. Therefore, the Rayleigh stability criterion only tells us that the base flow is linearly stable to inviscid waves 18 2. Non-modal growth in annular Poiseuille-Couette flow whenever 2Pm2 − Qk2 6= 0. However, since the Rayleigh argument is only a necessary condition for instability, we cannot deduce that the flow is unstable if 2Pm2 −Qk2 = 0. In order for 2Pm2 − Qk2 = 0 to be satisfied, P and Q must have opposite sign. Since 0 < η < 1, P is strictly positive. Therefore we require Q < 0, which implies V > V ∗ = 2(1− η2) log η 1− η2 + (1 + η2) log η . (2.25) Here V ∗ is precisely the maximum inner cylinder velocity Vmax achieveable in thread- annular flow, meaning that thread-annular flow is always stable to inviscid waves. Viscous waves can de-stabilize APCF and we numerically determine these instabilities by solving the linear system (2.7) for ω ∈ C, using the Chebyshev collocation tech- nique outlined in Schmid & Henningson (2001). Typically, the collocation technique uses Chebyshev polynomials of the first kind, Tn, to approximate functions of the non-periodic coordinate ξ, say, at discrete locations ξj called Gauss-Lobatto grid points ξj = cos ( (j − 1)pi N − 1 ) , f(ξj) = N−1∑ n=1 anTn−1(ξj), j = 1, 2, . . . , N. (2.26) Since the Gauss-Lobatto points lie between −1 ≤ ξj ≤ 1, we map from the annular gap r ∈ [ η1−η , 1 1−η ] to ξ ∈ [−1, 1] using rj = 1 2 ( 1 + η 1− η − ξj ) , η 1− η ≤ rj ≤ 1 1− η . (2.27) With derivatives of the dependent variables represented by Chebyshev differentiation matrices, this discretization allows us to solve the generalized 4N×4N eigenvalue problem Lqˆ` = ω`Mqˆ`, (2.28) for ω` and qˆ` = (u`, v`, w`, p`), where the operators A and B are given by L =       χ+ 1Re 1 r2 − dU dr 0 −ik 0 χ − 1Re 2im r2 − d dr 0 1Re 2im r2 χ − im r ik ddr + 1 r im r 0       , (2.29) M =       1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0       , (2.30) 2.2 Linear stability 19 and we have defined χ = 1 Re ( d2 dr2 + 1 r d dr − ( 1 +m2 r2 + k2 )) − ikU. (2.31) Our solutions to the linear stability equations are found by numerical computation of the full spectrum using the eigenproblem routines in MATLAB. A number of spurious modes are found by solving the system in this way. We guard against spurious modes entering our calculations in the following way: the rows of L which are modified to enforce the no-slip boundary condition on the cylinder walls are multiplied by an arbitrary constant A. Therefore the equation for u at the inner wall, for example, becomes Au`|r= η1−η = ω`u`|r= η 1−η . (2.32) This forces spurious eigenvalues to take the value ωs = A, and by choosing an appropriate A we can ensure spurious modes are kept out of our calculations. Having arrived at a set of modes, we list the eigenvalues in terms of decreasing real part, ω1, ω2, . . . , ω`, . . ., such that the primary mode corresponds to the least stable. In Figure 2.3 we plot the eigenvalue spectra for the case (k,m, η,Re) = (1, 1, 0.5, 2000), at various inner cylinder velocities V . We use the computations of Walton (2005) to validate our numerical routines, as their pa- per includes calculations of neutral curves for radius ratio η = 0.4, over various azimuthal wavenumbers and inner cylinder velocities. Using a bisection technique, we calculate the corrsponding neutral curves using our non-dimensionalization, and compare them with Walton’s results. In Figure 2.4 we compare our calculations with those of Walton (2005) by re-printing Figure 4(a) from Walton’s paper alongside our own computations re-scaled to the variables used by Walton. We find good agreement between our results and those of Walton, with little noticeable difference between the curves. 20 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.3: Complex eigenvalues ω` of the linear stability equations for (k,m, η,Re) = (1, 1, 0.5, 2000) at the inner cylinder velocities V = 0, 0.1, 1,−0.1 as indicated. Figure 2.4: On the left we re-print Walton’s neutral curves for η = 0.4, m = 1 and V = 0 (solid), V = 0.015 (dotted), V = 0.03 (dashed) and V = 0.1 (dot-dashed). The axial wavenumber in Walton’s paper is denoted α. On the right we plot our calculations re-scaled to Walton’s quantities, with Vw = 0 (solid) and Vw = 0.03 (dashed). 2.2 Linear stability 21 Further to the (Re, k) curves, we compute neutral curves in the (Re, V )-parameter space by maximizing unstable growth rates over axial wavenumber k for fixed azimuthal wavenumber m. In particular we plot neutral curves for the radius ratios considered in §2.4: η = 0.2, 0.5, 0.8. We could not find any linearly unstable modes for Re ≤ 105 when η = 142 . Our results are shown in Figures 2.5, 2.6 and 2.7. In each Figure, we plot neutral curves at the least stable azimuthal orders for Re < 105. Across all radius ratios, |V | & 1 is stablizing, as the critical Re is pushed higher than its V = 0 value. V small and neg- ative is de-stabilizing in the cases η = 0.2 and η = 0.5, as shown in Figures 2.5 and 2.6. From Figure 2.7 we see that small positive V is de-stabilizing for η = 0.8. We find the lowest critical (Re, V ) combination when η = 0.2 for (Re, V ) = (6441,−0.15). Walton (2005) reports disconnected neutral axisymmetric neutral curves; our calculations are in agreement with this, though our results do not include Re large enough to support dis- connected axisymmetric neutral curves. Walton additionally mentions that they do not find disconnected non-axisymmetric neutral curves. Our calculations support this con- clusion, however, we find that the m = 1 neutral curve for η = 0.5, and all neutral curves for η = 0.8 are scalloped, with three local turning points in (Re, V ), for V < 0, V ≈ 0 and V > 0. Sadeghi & Higgins (1991) found asymptotic cutoff velocities for 0.5 < η < 1 for axisymmetric and non-axisymmetric modes and positive inner cylinder velocities, adapt- ing the methods used by Cowley & Smith (1985) for planar Poiseuille-Couette flow. In Figure 2.6 we plot the cutoff velocities for η = 0.5, as listed in their paper, alongside our neutral curves. Our computations are consisent with their asymptotic findings. Figure 2.5: Left: Neutral curves in (Re, V ) at azimuthal order m = 1 and radius ratio η = 0.2. Right: Critical axial wavenumbers k corresponding to the neutral modes in the left panel. 22 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.6: Left: Neutral curves in (Re, V ) at azimuthal orders m = 1, 2 and radius ratio η = 0.5. The dashed lines indicate the cutoff velocities derived by Sadeghi & Higgins (1991). Right: Critical axial wavenumbers k corresponding to the neutral modes in the left panel. Figure 2.7: Left: Neutral curves in (Re, V ) at azimuthal orders m = 0, 1, 2, 3, 4 and radius ratio η = 0.8. Right: Critical axial wavenumbers k corresponding to the neutral modes in the left panel. 2.3 Non-modal growth 23 2.3 Non-modal growth Our solutions to the linear stability equations are used to compute the non-modal energy growth in APCF. We solve for ω` ∈ C and vector qˆ` = (uˆ`, pˆ`) satisfying the linear system (2.28) Lqˆ` = ω`Mqˆ`. (2.33) An arbitrary velocity disturbance can be expanded in terms of the first L modes of the stability equations u(x, r, θ, t) = ei(kx+mθ) L∑ `=1 a`uˆ`(r)eω`t, (2.34) for a` ∈ C. We then introduce the norm E(t) = ‖u(t)‖2 = ∫ 1 1−η η 1−η uHu rdr, (2.35) which is related to the physical kinetic energy of the fluid, E˜, in the periodic domain [0, 2pik ]× [ η 1−η , 1 1−η ]× [0, 2pi m ] by the relation E = mkE˜ 2ρpi2 , (2.36) where ρ is the constant fluid density. For aˆ = (a1, a2, . . . , aL) and Λ = diag(ω1, ω2, . . . , ωL), E(t) can be written E(t) = aˆHeΛ H tMeΛtaˆ, (2.37) where eΛt is the matrix exponential of Λt, and M is the L× L matrix Mjk = ∫ 1 1−η η 1−η uˆHk uˆj rdr. (2.38) M is both Hermitian (M = MH) and positive definite. Since M is Hermitian, its singular value decomposition can be written M = UMΣMU H M , (2.39) where ΣM = diag(σ1, σ2, . . . , σL) is a diagonal matrix of the singular values of M in descending order σ1 > σ2 > . . . > σL, and UM is a unitary matrix. For √ ΣM = diag( √ σ1, √ σ2, . . . , √ σL) we have F = UM √ ΣM , (2.40) 24 2. Non-modal growth in annular Poiseuille-Couette flow which allows us to factor M into M = FHF . Hence, we have ‖u(t)‖2 = aˆHeΛ H tFHF eΛtaˆ = ‖F eΛtaˆ‖22, (2.41) where ‖.‖2 is the Euclidean 2-norm. The maximum energy amplification G is simply G(t) = max u0 6=0 E(t) E(0) = max u0 6=0 ‖u(t)‖2 ‖u0‖2 , (2.42) = max aˆ6=0 ‖F eΛtaˆ‖22 ‖F aˆ‖22 , (2.43) = max aˆ6=0 ‖F eΛtF−1F aˆ‖22 ‖F aˆ‖22 , (2.44) = ‖F eΛtF−1‖22. (2.45) G(0) = 1 by definition. Linear non-modal energy growth is supported in a flow if G(t) > 1 for any non-zero t. G(t) represents the maximum possible energy amplification at time t of an initial perturbation, though different initial perturbations will be responsible for maximum growth at distinct times. Suppose, having found the maximum energy amplification achieveable at time tm, we wish to find the form of the optimal initial perturbation and also its form at tm. Let u0 denote the optimal perturbation, with u0 = ei(kx+mθ) L∑ `=1 a`uˆ`, and ‖u0‖ = 1. (2.46) Then, since u0 will have energy ‖F eΛtmF−1‖22 at tm, we have Baˆ = ‖F eΛtmF−1‖2aˆe Λtm = µaˆeΛtm , (2.47) where B = F eΛtmF−1. The constant µ = ‖F eΛtmF−1‖2 is simply the principal singular value of B. The singular value decomposition of B allows us to write BVB = UBΣB, (2.48) for unitary matrices VB and UB, and the matrix of singular values of B, ΣB = diag(σ1, σ2, . . . , σL). If we denote the principal eigenvectors of VB and UB as v1 and u1, respectively, then since µ = σ1 we have Bv1 = σ1u1, (2.49) from which we can identify aˆ = v1. The linear optimal u0 can then be found by ex- panding aˆ over the basis vectors uˆ1 . . . uˆL. The optimal at time tm is found by simply expanding aˆ over the vectors uˆ1eω1tm . . . uˆLeωLtm 2.3 Non-modal growth 25 Given that the gain G = G(t, k,m, η, Re, V ) is a six parameter function, we introduce two maxima which allow us to represent the level of non-modal growth possible for a perturbation at a fixed Reynolds number, inner cylinder velocity and radius ratio. First, we maximize G over time Gmax(k,m, η,Re, V ) = max t G(t, k,m, η, Re, V ) = G(tmax, k,m, η, Re, V ). (2.50) Gmax represents the greatest linear energy amplification that a perturbation with pre- scribed wavenumbers can experience. We further maximize over the wavenumbers k and m to define G(η,Re, V ) = max k,m Gmax(k,m, η,Re, V ) = G(Tmax, kmax,mmax, η, Re, V ), (2.51) such that G denotes the maximum linear energy growth possible for a given Re, V and η combination. kmax, mmax and Tmax correspond to the axial wavenumber, azimuthal wavenumber and the amplification time of the perturbation responsible for the ampli- fication G. If G = 1, then all perturbations are non-modally stable. We shall refer to the perturbation which achieves the maximum, optimal, growth G as the linear optimal. We use G as a measure of the non-modal stability of the flow, with large values of G interpreted as being de-stabilizing. The MATLAB golden mean search fminbnd is used to find Gmax and G. Finding Gmax requires a straightforward optimization over time; however, G may have local maxima, so we vary the time horizon at each step to increase the likelihood of finding Gmax. Since the azimuthal wavenumber m is integer valued, we consider a range of m in turn, and use the golden mean search to maximize over axial wavenumber k to find G. This is computationally more expensive than finding Gmax, as a new spectrum of modes must be computed for each wavenumber combination (m, k). In order to reduce the number of wavenumber searches required to find G, we assume that the optimal wavenumber mmax varies smoothly with Re, V and η. To check against this assumption, we perform random spotchecks of Gmax over a wider range of wavenumbers than typically considered with Re, V and η fixed. 26 2. Non-modal growth in annular Poiseuille-Couette flow 2.4 Optimal growth in APCF We begin our analysis of non-modal growth by plotting curves of Gmax in the (Re, k)- plane, to investigate how transients interact with linearly unstable modes. The results in Figure 2.8 show that, while changes in V have a significant effect on the linearly unstable region, there is little qualitative effect on the transient behaviour. The curves of Gmax approach the unstable boundary smoothly, as illustrated in Figure 2.9. Isolines of G in (Re, V )-space are presented in Figure 2.10 at the radius ratios η = 1 42 , 0.2, 0.5, 0.8. Each plot is generated from a grid of 31 × 115 grid-points in (Re, V ). The radius ratio η = 142 ≈ 0.0238 was chosen to coincide with the ratio considered by Shands et al. (1980), for which there is experimental data on the transitional values of Re and V . We will make a comparison to this data in §2.5. The isolines are asymmetrical with respect to postive or negative inner cylinder velocity, with a turning point at slightly negative V . However, the asymmetry in V becomes less pronounced as η increases. The qualitative trend is that non-zero inner cylinder velocity has a de-stabilizing influence on the base flow as, aside from the small stabilizing region of V < 0, G is found to increase with |V |. This is in contrast to the modal linear stability results from §2.2, where a non- zero inner cylinder velocity has a stabilizing effect. The de-stabilizing trend is robust across all the radius ratios considered, with little change in the maximum energy growth achieved by each optimal. The thick red curves in Figure 2.10 are the energy stability boundaries, defined for V fixed to be the maximum Re for which G = 1. We compute the energy stability boundary using a bisection method, for fixed V , to calculate the energy- critical Re to an accuracy of  = 10−6. To the left of the energy stability curve, the energy of all disturbances is strictly monotically decaying. In parallel with results from the canonical shear flows, the Reynolds numbers on the energy stability boundary, ReE are many orders of magnitude smaller than the critical Re from linear stability theory. In Figure 2.11 we plot the optimal azimuthal wavenumbers mmax at each point in the (Re, V )-plane. These plots invariably look less well resolved than the contours of G, but this is a consequence of the integer nature of mmax. The same grid is used in the calculation of both Figures. Typically mmax varies in unit steps, however for η = 142 when Re > 103.5, mmax jumps from one to three as V increases from below to V ≈ −0.63. mmax then falls back to two as V is increased further. This sequence is depicted for Re = 104 in Figure 2.12. We find that the values of mmax increase as the radius ratio η approaches unity. For η = 142 we have 1 ≤ mmax ≤ 3, whereas for η = 0.8 we find 13 ≤ mmax ≤ 19. Lower azimuthal orders are selected when |V | becomes large. The behaviour of kmax also changes for |V | large: the optimals are axially independent for 2.4 Optimal growth in APCF 27 small inner cylinder velocities, but weakly oblique disturbances are optimal for large |V |, as is shown in Figure 2.13. The scenario is reminiscient of planar flows, in that streamwise independent disturbances maximize growth for plane Poiseuille flow, whereas slightly oblique disturbances maximize the transient growth in plane Couette flow (Butler & Farrell (1992), Reddy & Henningson (1993)). We find that when kmax is non-zero, kmax ∝ Re−1 at fixed V . Examples of this relationship are displayed in Figure 2.14. Figure 2.8: Contours of Gmax. Top row: η = 0.2, m = 1. Bottom row: η = 0.5, m = 2. The red curves in each panel mark the energy stability boundary. The black curves are contours of Gmax, from log10Gmax = 0.25, 0.5, . . . , 4.5. The black shaded region is linearly unstable. Figure 2.9: Gain curves showing the approach to the neutral curve. For k > kc G is unbounded in the limit t→∞. (Re, V, η,m) = (104.5, 0, 0.5, 2), kc = 0.688. 28 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.10: Contours of maximum linear optimal growth G in the (Re, V )-plane for radius ratios η = 142 , 0.2, 0.5, 0.8. The thick red curve demarcates the energy stability boundary. The blue curves are contours of log10 G = 0.25, 0.5, 0.75, . . . , 6. 2.4 Optimal growth in APCF 29 Figure 2.11: Colour-plots of azimuthal wavenumbers mmax corresponding to G in the (Re, V )-plane for radius ratios η = 142 , 0.2, 0.5, 0.8. For each η, mmax increases for |V | small. Figure 2.12: G for (Re, η) = (104, 142) is the supremum of the intersection of the depicted curves. A non-unit change in mmax from one to three can be seen for V ≈ −0.63. 30 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.13: Optimal wavenumbers kmax and mmax against V , for Re = 102 (blue), Re = 103 (red) and Re = 104 (black), at the radius ratios η = 0.2 and η = 0.5. Non-zero axial wavenumbers are preferred for |V | large, and mmax increases for |V | small. Figure 2.14: Where kmax is non-zero, it follows kmax ∝ Re−1 with V fixed. 2.4 Optimal growth in APCF 31 2.4.1 Optimal perturbations and growth mechanism We find that both the axially independent (kmax = 0) and the weakly oblique (kmax ∝ Re−1) optimals have energy gains that scale with the square of the Reynolds number G ∝ Re2. In addition, we find that the time taken to reach the optimal, Tmax, scales with the Reynolds number as Tmax ∝ Re. A typical example is plotted in Figure 2.15. These scalings are suggestive of energy amplification via the lift-up effect, a well-documented inviscid, algebraic growth mechanism (Landahl (1980)), that has been found to drive op- timal growth in viscous planar flows (Gustavsson (1991), Butler & Farrell (1992), Reddy & Henningson (1993), Schmid & Henningson (2001), Chapman (2002)) and in Hagen- Poiseuille flow (Bergstro¨m (1993), Schmid & Henningson (1994)). Physically, the lift-up mechanism refers to extraction of energy from the mean flow as low-energy vortices are ‘lifted-up’ into high-energy streaks in the direction of the mean flow. Perturbations must be three-dimensional to experience this growth, as the mechanism involves wall-normal velocity interacting with the base flow in such a way that wall-normal vorticity is ampli- fied. To confirm that the lift-up mechanism is responsible for the energy growth of the opti- mals in APCF, we analyse the evolution of the energy fraction of an arbitrary optimal disturbance. The energy fractions are defined Eu = ∫ 1 1−η η 1−η u∗u rdr ∫ 1 1−η η 1−η uHu rdr , Ev = ∫ 1 1−η η 1−η v∗v rdr ∫ 1 1−η η 1−η uHu rdr , Ew = ∫ 1 1−η η 1−η w∗w rdr ∫ 1 1−η η 1−η uHu rdr , (2.52) and can be calculated at an arbitrary time for a given velocity perturbation. We plot the evolution of the energy fractions for two optimals in Figure 2.16: one axially indepen- dent and one weakly oblique. We find that in each case energy is transferred to the axial velocity. We take this as confirmation of the lift-up effect being the amplification mecha- nism driving the transient growth. Combined with the results of Figure 2.10, a non-zero inner cylinder velocity can be interpreted as enhancing the lift-up mechanism by causing greater linear energy growth. In Figures 2.17, 2.18, 2.19 and 2.20 we plot cross-sectional velocity fields of a selection of optimal perturbations at t = 0 and t = Tmax. Each Figure consists of a vector plot of the (v, w) velocity field and a colour-plot of the axial velocity u, with positive and negative velocities coloured red and blue respectively. In Figure 2.17 we observe a localization feature of the optimals whereby the streaks at Tmax have moved towards the inner cylinder. We find this feature only to be present for small radius radio and in cases with large positive or negative inner cylinder velocity, The optimal in Figure 2.18 has an inner and outer streak layered structure, with the inner ring of streaks 32 2. Non-modal growth in annular Poiseuille-Couette flow strengthening between t = 0 and Tmax. Our calculations suggest this structure is typical of the optimals with fixed or low velocity inner cylinder. As the radius ratio increases in Figures 2.19 and 2.20, the increase in azimuthal order is evident in the greater number of vortices and streaks in the optimal. Figure 2.15: G and Tmax against Re for radius ratio η = 0.2 and V = 0 (left) and V = 10 (right). The red dashed lines represent the asymptotic scalings G ∝ Re2 and Tmax ∝ Re. Figure 2.16: Energy fractions showing the emergence of streaks in the optimals for both V = 0 and V = 10. The optimal for V = 0 is axially independent with (kmax,mmax, Tmax) = (0, 3, 21.84), and the optimal for V = 10 is weakly oblique with (kmax,mmax, Tmax) = (0.0312, 2, 25.69). (Re, η) = (103, 0.2) in both cases. 2.4 Optimal growth in APCF 33 Figure 2.17: Optimal perturbations at t = 0 and t = Tmax = 19.15 for (Re, V, η, kmax,mmax) = (103,−5, 142 , 0.2666, 1). Figure 2.18: Optimal perturbations at t = 0 and t = Tmax = 21.84 for (Re, V, η, kmax,mmax) = (103, 0, 0.2, 0, 3). 34 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.19: Optimal perturbations at t = 0 and t = Tmax = 34.22 for (Re, V, η, kmax,mmax) = (103, 5, 0.5, 0, 5). Figure 2.20: Optimal perturbations at t = 0 and t = Tmax = 30.26 for (Re, V, η, kmax,mmax) = (103, 10, 0.8, 0.0214, 15). 2.4 Optimal growth in APCF 35 2.4.2 Axisymmetric growth The azimuthal velocity in the linear system (2.7) decouples from the axial and radial ve- locities for axisymmetric perturbations (m = 0). Therefore, axisymmetric perturbations cannot be amplified by the three-dimensional lift-up mechanism since momentum can only be exchanged between the axial and radial velocities. We find instead that optimals grow via the Orr mechanism (Orr (1907)), typified by perturbations which are initially inclined against the mean flow and gain energy by evolving to align themselves with the mean shear. The Orr mechanism is a much weaker effect than the lift-up mechanism, with characteristic Re scalings G ∝ Re 2 3 and Tmax ∝ Re 1 3 . Optimal growth of oblique modes in plane Couette and plane Poiseuille flows was shown to have this scaling by Chapman (2002) and Heaton & Peake (2007), who additionally showed that the same scalings are recovered in a cylindrical geometry. Axisymmetric modes in APCF recover the aforementioned Re scalings; we demonstrate this in Figure 2.22. Since the axisymmetric modes grow on a different scale to non- axisymmetric optimals, we repeat the contour plots of Figure 2.10 to explore how ax- isymmetric transients behave in (Re, V )-space. Gm=0 is calculated as described in §2.3 but with wavenumber optimization over k only, and m fixed at zero. Contours of Gm=0 are plotted in Figure 2.21. The effect of inner cylinder velocity on Gm=0 is generally the same as for G, in that large |V | creates an increase in the energy amplification. However, each contour has a maximum Re for Vm with Vm < 0. Vm decreases as contour size in- creases, and the effect is more pronounced at smaller radius ratios. This means that there is a region of negative V which has a stabilizing influence on the transient behaviour. In particular, for η = 142 and Gm=0 > 10, the region 0 < V < −15 is stabilizing with respect to the stationary inner cylinder since the contours of Gm=0 > 10 do not bend back to their value at V = 0, as can be seen in Figure 2.21. In all cases, V > 0 has a de-stablizing effect on the transients. The axial wavenumbers kmax from our calculations of Gm=0 are always O(1). The wavenumbers experience a sharp turning point when plotted against V , at V ≈ 0 for fixed Re, and we display some examples of this in Figure 2.23. Additionally, the curves in Figure 2.23 appear to flatten out for |V | large, which is perhaps an indication of Couette shear dominating over pressure-driven forces. 36 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.21: Contours of maximum axisymmetric growth Gm=0 in the (Re, V )-plane for radius ratios η = 142 , 0.2, 0.5, 0.8. The thick red curve demarcates the axisymmetric energy stability boundary. The blue curves are contours of log10 Gm=0 = 0.25, 0.5, 0.75, . . . , 2. 2.4 Optimal growth in APCF 37 Figure 2.22: G and Tmax against Re for radius ratio η = 0.5 and V = −10 (left) and V = 15 (right). The red dashed lines represent the asymptotic scalings G ∝ Re 2 3 and Tmax ∝ Re 1 3 . Figure 2.23: Wavenumbers kmax for the axisymmetric optimals Gm=0. 38 2. Non-modal growth in annular Poiseuille-Couette flow 2.4.3 Radius ratio dependence of G We now investigate how G changes as the radius ratio is varied over the range 0 < η < 1. From Figures 2.10 and 2.21, we saw that the contours of G and Gm=0 were relatively unchanged over the radius ratios considered in the Figures, though the wavenumbers responsible for the maximum growth changed. This is confirmed in Figure 2.24, where G is plotted against η at various inner cylinder velocities. The curves have a scalloped appearance, owing to a succession of azimuthal orders contributing to the optimal. In each plot, in spite of quantitative changes in G as η varies, the order of magnitude of G is unchanged, meaning that the maximum linear energy amplification possible is not strongly dependent on the specific geometry of the annulus considered. Figure 2.24: G against η at various inner-cylinder velocities. Re = 2000. Dashed lines indicate the locations where mmax changes value. Furthermore, we can determine the optimal radial-azimuthal aspect ratio arθ from our computational data. We define arθ to be the ratio of the wall-normal gap and the average azimuthal periodic length of the optimal disturbance arθ = B − A 2pi mmax (B+A2 ) = mmax pi ( B − A B + A ) = mmax pi ( 1− η 1 + η ) . (2.53) 2.4 Optimal growth in APCF 39 Generally arθ converges on a constant value as η → 1, indicating that the optimals prefer a certain parameter dependent aspect ratio as they approach the narrow gap limit. A particularly interesting case is when V = 0, where we find that as the radius ratio approaches unity, arθ approaches the wall-normal-spanwise aspect ratio ayz of the optimal perturbation in plane Poiseuille flow. This is in keeping with the notion that annular Poiseuille flow approaches plane Poiseuille flow in the narrow-gap limit. We plot our findings in Figure 2.25, where we have calculated ayz for plane Poiseuille disturbances based on an optimal spanwise wavenumber β ≈ 2, which gives ayz = 2 ( 2pi β ) = 2 ( 2pi 2 ) = 2 pi ≈ 0.6366. (2.54) Figure 2.25: The optimal radial-azimuthal aspect ratio arθ is plotted against the radius ratio η (blue dots), for (Re, V ) = (2000, 0). It can be seen that arθ converges towards a constant as η → 1. The red dashed line denotes the value of the wall-normal-spanwise aspect ratio ayz = 2pi for optimal disturbances in plane Poiseuille flow. We infer that in the narrow gap limit the optimal disturbances of annular Poiseuille share the characteristics of optimal disturbances in plane Poiseuille flow. 40 2. Non-modal growth in annular Poiseuille-Couette flow 2.5 Comparison to experimental data Experimental data relating to APCF is relatively scarce, with Shands et al. (1980) and Frei et al. (2000) being the only authors conducting experiments in APCF that we are aware of. Frei et al. (2000) are primarily concerned with thread-annular flow, and explicit data on the parameters at which the flow transitions is not given in their work. In the Shands et al. (1980) experiments they do report specific (Re, V ) values at which the lam- inar flow was observed to break down; we therefore focus on data from their experiments. Shands et al. (1980) gauge the transition to turbulence of the fluid in their experiments by measuring the friction factor λ = 2∆pD ρLÛ2 . The quantities ∆p, D, ρ, L and Û are defined in §2.1. From equation (2.9b) we have that under laminar flow, the total mean velocity, Û , must satisfy Û = ∆p 8µL log BA ( B2 − A2 − (B2 + A2) log B A ) + V˜ ( B2 − A2 − 2A2 log BA 2(B2 − A2) log BA ) . (2.55) Dividing through by Û gives 1 = ∆pB2 8ρνLÛ ( 1− η2 + (1 + η2) log η log η ) + V˜ Û ( η2 − 1− 2η2 log η 2(1− η2) log η ) (2.56) 1 = 1 32 ( ∆pD ρLÛ2 )( ÛD ν )( 1− η2 + (1 + η2) log η (1− η)2 log η ) + V˜ Û ( η2 − 1− 2η2 log η 2(1− η2) log η ) . (2.57) Hence for laminar flow λ, R̂e and V̂ have the following relationship 1 = λR̂e 64 ( 1− η2 + (1 + η2) log η (1− η)2 log η ) + V̂ ( η2 − 1− 2η2 log η 2(1− η2) log η ) . (2.58) At fixed R̂e, Shands et al. (1980) measure the V̂ > 0 for which (2.58) first breaks down, thereby determining the point at which flow transitions away from laminar flow. In Fig- ure 2.26 we re-print Fig. 3 from Shands et al. (1980), which shows the friction factor λ plotted against V̂ for R̂e = 1100, 1900 and gap ratio η = 142 . In the Figure, mea- surements of the friction factor fall along equation (2.58) until transition occurs when V̂ = 3.3, 2 for R̂e = 1100, 1900, respectively. Additionally, Shands et al. (1980) report transition when V̂ = 2.7 and R̂e = 1300, and transition in the fixed inner cylinder case for R̂e = 2000. In our non-dimensionalization, these values correspond to Re = 308, 416, 697 and V = 5.89, 4.22, 2.73. Without any explicit comment on how linear optimal perturbations may contribute to 2.5 Comparison to experimental data 41 transition away from the laminar base flow, we plot the values at which Shands et al. (1980) report transition for η = 142 , alongside our curves of G in Figure 2.27. The transi- tion values lie within 102 < G < 102.5 and, suggestively, greater inner cylinder velocity is de-stabilizing in both the experiments and in the sense that it increases G. Since linearly unstable modes appear to be, in general, stabilized by inner-cylinder velocity, we draw the tentative conclusion that the energy amplification mechanisms driving the linear op- timal growth play a greater role in the break down of the laminar flow than the linearly unstable modes. Furthermore, Shands et al. (1980) remark that, having repeated their experiments for different radius ratios, transition is found to occur at similar parameter values. This is in keeping with our finding in §2.4.3 that the maximum energy amplifi- cation G is not strongly sensitive to changes in the radius ratio. 42 2. Non-modal growth in annular Poiseuille-Couette flow Figure 2.26: Re-print of Fig. 3 from Shands et al. (1980). The friction factor λ is plotted against inner cylinder velocity V̂ = V˜ Û for R̂e = 1100 (squares) and R̂e = 1900 (triangles). The two straight lines indicate the laminar relationship (2.58), which breaks down when V̂ = 3.3, 2 and R̂e = 1100, 1900 respectively. Figure 2.27: Contours of G in (Re, V )-space for η = 142 are plotted alongside tran- sition values observed in Shands et al. (1980). The contour values are log10 G = 0.25, 0.5, . . . , 5.25 with contours log10 G = 1, 2, 3, 4, 5 labelled. The black triangles mark the points of transition. The red curve demarcates the energy stability boundary. 2.6 Conclusion 43 2.6 Conclusion This chapter began with a re-cap of the linear stability properties of APCF, which we augmented by calculating neutral curves in (Re, V )-space. We find that non-zero inner cylinder velocity V tends to have a stabilizing influence on the modal stability of pertur- bations to the base flow. We then conducted a non-modal stability analysis, using the maximum energy amplifi- cation over all wavenumber combinations as a measure of the non-modal instability for an arbitrary (Re, V ) pair. We found transient energy amplification driven by the lift-up mechanism is possible for non-axisymmetric waves. The order of magnitude of the max- imum growth is not found to be strongly correlated to the radius ratio of the cylinders. Axisymmetric waves experience energy amplification via the Orr mechanism, though this amplification is not as strong as the lift-up mechanism. However, the maximum energy amplification associated with both axisymmetric and non-axisymmetric modes is found to increase with increasing inner cylinder velocity magnitude, which implies that non- zero inner cylinder velocity has a de-stabilizing influence on perturbations to the mean flow, in contrast with the modal stability results. Finally we compare our results with the study of Shands et al. (1980), who studied APCF in an experimental set up, and indicated the (Re, V ) parameters at which they observed transition to turbulence. The transitional values lie within a few contours of maximum growth, suggesting a relationship between the linear energy growth mechanisms and a break down in the laminar flow. 44 2. Non-modal growth in annular Poiseuille-Couette flow Chapter 3 Stability and bifurcations in rotating plane Couette flow In this chapter we focus on the stability and bifurcation behaviour of structures in super- critical rotating plane Couette flow (RPCF). The governing equations and the geometry of the problem are introduced in §3.1. In §3.2 we review the linear stability properties of the basic flow. In §3.3 an account is given of the methods we use to find nonlinear nu- merical secondary and tertiary solutions to the governing equations. In §3.4 we describe the Floquet theory techniques which we use to determine the stability of spatially and temporally periodic states. §3.5 is a discussion of the secondary states which bifurcate from the base flow. In §3.6 we perform a Floquet stability analysis of Taylor vortex flow, an important secondary structure in RPCF. Bifurcations of tertiary solutions, which emerge following the break-down of Taylor vortices, are presented in §3.7. Finally, in §3.8 we investigate the relevance of the secondary and tertiary states to the transition phenomena observed in physical experiments of RPCF, and propose a simple model with which the transition process can be analyzed. 3.1 Governing equations Rotating plane Couette flow (RPCF) is the flow of an incompressible, Newtonian fluid under linear shear and Coriolis rotation. The fluid is constrained between moving im- permeable walls; with velocity difference 2U˜ between the upper and lower walls and channel width 2h˜. We define unit vectors i, j and k pointing in the streamwise, wall- normal and spanwise directions respectively, as shown in Figure 3.1. Working towards a non-dimensionalized system of equations, we first consider the dimensional momentum balance in a rotating frame of reference ∂u˜ ∂t˜ + u˜ · ∇˜u˜ = − 1 ρ˜ ∇˜p˜+ ν˜∇˜2u˜+ 2 ( u˜× Ω˜ ) − Ω˜× Ω˜× x˜− dΩ˜ dt˜ × x˜, (3.1) 45 46 3. Stability and bifurcations in rotating plane Couette flow where ρ˜ is the fluid density, ν˜ is the kinematic viscosity and x˜ is a position vector pointing from the axis of rotation. For steady rotation about the z-axis we have Ω˜ = Ω˜k, which allows us to drop the time-dependent Euler force term dΩ˜ dt˜ × x˜. Since the centrifugal force term, Ω˜× Ω˜× x˜, contains no dynamic variables, we introduce the pressure p¯ p¯ = p˜− ρ˜Ω˜2 2 x˜2, (3.2) where x˜2 = x˜ · x˜. This allows the momentum equations to be expressed as ∂u˜ ∂t˜ + u˜ · ∇˜u˜ = − 1 ρ˜ ∇˜p¯+ ν˜∇˜2u˜+ 2 ( u˜× Ω˜ ) . (3.3) We now introduce three non-dimensional parameters: the Reynolds number Re, the rotation number Ro and a second rotation number Ω Re = U˜ h˜ ν˜ , Ro = 2Ω˜h˜ U˜ , Ω = 2Ω˜h˜2 ν˜ = Re Ro. (3.4) Hence, after non-dimensionalization and inclusion of the continuity equation, the velocity and pressure (u, p) satisfy ∂u ∂t + u · ∇u = −∇p+ 1 Re ∇2u+Ro ( u× k ) , (3.5a) ∇ · u = 0, (3.5b) with the no-slip boundary condition imposed on the each of the walls. The linear shear profile U = yi, is a solution to this system with pressure distribution P = P0 − 12Ro y 2, where P0 is the ambient pressure. Of the two rotation numbers Ro and Ω , Ro is useful since it appears as a single coefficient controlling the Coriolis force in (3.5), whereas Ω is independent of the wall velocity in terms of its non-dimensionalization meaning that it describes the effects of rotation independent of Re. For the majority of this chapter we describe rotation in terms of Ro, however, in §3.8 we use Ω , in accordance with experimental preference. Figure 3.1: Sketch of the flow geometry. 3.2 Linear stability 47 3.2 Linear stability It has been shown that RPCF is linearly unstable by Lezius & Johnston (1976), who drew on the equivalence of the linearized system (3.5) to two-dimensional Rayleigh- Be´nard convection. The flow loses stability in the rotation parameter range 0 < Ro < 1, to spanwise periodic disturbances with wavenumber β ≈ 1.55. The neutral curve takes the form Re2 = 107 Ro(1−Ro) . (3.6) More generally we can solve directly for the eigenfrequencies and eigenmodes of the linearized system using a pseudo-spectral technique. We Fourier transform the velocity vectors in the directions x and z and assume exponential dependence in time, with Fourier coefficients α, β and ω respectively. The fluid vectors become [u, v, w, p](x, y, z, t) = [u, v, w, p](y)ei(αx+βz)+ωt. (3.7) With base flow, or primary state, U = U(y)i = yi, the linear stability equations are ωu+ U iαu+ vU ′ = −iαp+ 1 Re ( u′′ − (α2 + β2)u ) +Rov (3.8a) ωv + U iαv = −p′ + 1 Re ( v′′ − (α2 + β2)v ) −Rou (3.8b) ωw + U iαw = −iβp+ 1 Re ( w′′ − (α2 + β2)w ) (3.8c) iαu+ v′ + iβw = 0. (3.8d) The system can be cast into a generalized eigenvalue problem Aq = ωBq, (3.9) for q = (u, v, w, p) and matrix operators A and B given by A =       A Ro− dUdy 0 −iα −Ro A 0 − ddy 0 0 A −iβ iα ddy iβ 0       , (3.10) B =       1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0       . (3.11) 48 3. Stability and bifurcations in rotating plane Couette flow where we have defined A = 1 Re ( d2 dy2 − (α2 + β2) ) − iαU. (3.12) We solve for ω by expressing q as a truncated series of Chebyshev polynomials and dis- cretizing the system accordingly. The technique is outlined in greater detail in Schmid & Henningson (2001). Above the critical Reynolds number Recr ≈ 20.7 there is a finite range of Ro or, equiva- lently, Ω for which the flow is linearly unstable. Figure 3.2: Left: Neutral curve. Right: ωr denotes the real part of the leading eigenvalue for Re = 100, Ro = 0.1. Perturbations with α = 0 are least stable. Our calculations show that linearly unstable perturbations take the form of streamwise oriented vortices, with wavenumbers α = 0 and β non-zero or obliquely oriented vortices with both α and β non-zero. We find that, for fixed β, Re and Ro, the leading streamwise vortex (SV) eigenmode is generically more unstable than the leading oblique vortex (OV) eigenmode: ωr(0, β) > ωr(α, β). (3.13) 3.2 Linear stability 49 Figure 3.3: Velocity fields u of (a) SV and (b) OV perturbations in the (y, z)-plane (top) with x = 0 and ((c) and (d)) the (x, z)-plane (bottom) with y = 0 for (Re,Ro) = (50, 0.1) and (α, β) = (0.5, 2). 3.2.1 Flowfield visualization It is often useful to visualize a fluid structure, whether it is a stability eigenmode or nonlinear solution to the Navier-Stokes equations. However, it is not straightforward to display a fully three-dimensional flowfield in such a way that none of its features are lost. Throughout this work we use two-dimensional projections of flowfields to concisely give a visual sense of their character. In each Figure we have a colour-plot of the streamwise velocity u, and a vector plot of the wall-normal and spanwise velocities (v, w). To display the streamwise variation of a flowfield we fix y = 0 and plot the (x, z)-plane. To view the vortical structure of a flowfield we fix x = 0 and plot the (y, z)-plane. Each Figure has a true aspect ratio. 50 3. Stability and bifurcations in rotating plane Couette flow 3.3 Nonlinear states We refer to numerical solutions of the nonlinear Navier-Stokes equations as nonlinear states. For a given initial guess, a Newton-Krylov-hookstep algorithm is used to find a steady or time-periodic state. Typically, our initial guesses come from predictions given by linear stability theory. Calculations are made using two direct numerical simulators. SIMSON (Chevalier et al. (2007)) is used for numerical integration of the Navier-Stokes equations, and Channelflow (Gibson (2008)) is used to integrate the Navier-Stokes equa- tions, find equilibrium and time-periodic solutions, and to continue the solutions in pa- rameter space. Both codes use spectral discretization in the periodic directions, and discretization of the wall-normal direction in Chebyshev polynomials. Time integration is performed using a finite-differencing technique. Both codes have been widely validated by other authors, and we find excellent agreement between the two in our simulations. The Newton-Krylov-hookstep algorithm outlined in Viswanath (2007), and implemented within Channelflow, is used to determine equilibrium and time-periodic solutions to the Navier-Stokes equations. Let ∂u ∂t = LNSu, F tNS(u) = u+ ∫ t 0 LNSu(τ) dτ (3.14) be the Navier-Stokes equations and their time-t forward map. The equations are solved in perturbation formulation about the primary flow U = yi, so the operator LNS is defined as LNS = 1 Re ∇2u−∇p+Ro ( u× k ) −U · ∇u− u · ∇U − u · ∇u. (3.15) The Newton-Krylov-hookstep algorithm finds solutions by minimizing the residual G(u) G(u) = F TNS(u)− u. (3.16) The algorithm was used by other authors to find solutions in plane Couette flow, where the absence of linear instabilities makes finding an initial condition which converges towards a solution difficult (Viswanath (2007), Gibson et al. (2009)). Though the basins of attraction of solutions in our case will generally be greater than in non-rotating plane Couette flow due to supercriticality, the algorithm is still useful for its computational efficiency. We will give a brief outline of how the algorithm works. Writing the solution vector as x = (u, αs, βs, ωs), where αs and βs are the wavenumbers in the periodic directions and ωs is the frequency of a time-periodic state, the Newton step requires 3.3 Nonlinear states 51 solution of DG(xn)(xn+1 − xn) = −G(xn), (3.17) where G(u) has been appended to account for the invariance of the solution under trans- lations in its periodic directions x, z or t and guard against solution updates which are simply translations, and DG(x) is the Jacobian of the system. Equation (3.17) is solved for the updated solution guess xn+1 using the GMRES Krylov subspace method (Saad & Schultz (1986), Sa´nchez et al. (2004)). If the Newton step fails to decrease the residual below a prescribed tolerance, a hookstep algorithm is employed to find a solution update which does decrease the residual sufficiently. The hookstep algorithm uses the Lagrange multipliers technique to minimize ‖DG(xn)(xn+1 − xn) +G(xn)‖2, (3.18) subject to ‖xn+1 − xn‖2 < δ2, (3.19) so that, in effect, we are taking a reduced Newton step. In summary, in the Newton- Krylov-hookstep algorithm the classic Newton iteration scheme is augmented by a Krylov subspace technique to enhance computational efficiency and a hookstep algorithm to in- crease the accuracy of each solution update. More detail on the full algorithm can be found in Viswanath (2007) and the references cited above. All solutions found in this thesis have an accuracy of ‖G(u)‖ < 10−14. Having found a state with velocity field u = (u, v, w) we use the cross-flow energy Ecf as a measure of the nonlinear solution, Ecf (αs, βs, Re,Ro) = 1 2DV ∫ D v2 + w2dV. (3.20) where D is the x and z periodic domain [0, Lx]×[−1, 1]×[0, Lz] = [0, 2piαs ]×[−1, 1]×[0, 2pi βs ], and DV is its volume. Ecf is defined in this way so that it is independent of the domain. For example, no distinction is made between a given solution and the same solution on a doubled domain. The four parameter dependence of Ecf invites us to fix three parame- ters and investigate trajectories of solutions on a low-dimensional subspace. As such we typically fix the geometric parameters (αs, βs) and investigate bifurcations in Ro, with the remaining dynamic parameter Re fixed. 52 3. Stability and bifurcations in rotating plane Couette flow 3.4 Floquet stability analysis Suppose we have found a nonlinear Navier-Stokes solution, that is periodic in any of the directions x and z or time t. Its stability properties can be determined through Floquet analysis. 3.4.1 Steady states First considering steady solutions, we utilize the spatial periodicity of the state by ex- pressing it as a Fourier series in the periodic variables. Given wavenumbers αs in x, βs in z; we have U s(x, y, z) = ∞∑ j=−∞ ∞∑ k=−∞ U jk(y)ei(jαsx+kβsz), (3.21) where U jk(y) = αsβs 4pi2 ∫ 2pi αs 0 ∫ 2pi βs 0 U s(x, y, z)e−i(jαsx+kβsz) dxdz. (3.22) We insert the Fourier series expansion of U s into the linearized Navier-Stokes equations in the form ∂u ∂t +U s · ∇u+ u · ∇U s = −∇p+ 1 Re ∇2u+Ro ( u× k ) , (3.23a) ∇ · u = 0. (3.23b) Following Floquet’s theorem, we write the perturbation vector q = (u, v, w, p) as q(x, y, z, t) = ei(αx+βz)+σt ∞∑ m=−∞ ∞∑ n=−∞ qˆmn(y)e i(mαsx+nβsz), (3.24) with the parameters α, β and σ free to be tuned away from the fundamental wavenumbers and frequencies. Using the expressions for q and U s we can match exponents in the perturbation equations and derive an infinite system of coupled equations of the form Amnqˆmn + ∑ j,k F jkmnqˆm−j,n−k = σBqˆmn, ∀m,n ∈ Z. (3.25) The operators A, B and F are defined as Amn =       A dU0dy −Ro 0 i(α +mαs) Ro A+ dV0dy 0 d dy 0 dW0dy A i(β + nβs) i(α +mαs) ddy i(β + nβs) 0       , (3.26) 3.4 Floquet stability analysis 53 F jkmn =       ijαsUjk + F dUjk dy ikβsUjk 0 ijαsVjk dVjk dy + F ikβsVjk 0 ijαsWjk dWjk dy ikβsWjk + F 0 0 0 0 0       , (3.27) B = −       1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0       , (3.28) where we have used A = i(α+mαs)U0 +V0 d dy +i(β+nβs)W0− 1 Re ( d2 dy2 − (α+mαs) 2− (β+nβs) 2 ) , (3.29) F = i(α + (m− j)αs)Ujk + Vjk d dy + i(β + (n− k)βs)Wjk, (3.30) and U jk = (Ujk, Vjk,Wjk), U 0 = U 00. (3.31) We can then prescribe wavenumbers α and β and, provided |U jk| < 1 to guarantee convergence, solve by truncation and discretization for σ. We truncate the system in the periodic directions with −Tm < m < Tm, −Tn < n < Tn. (3.32) The wall-normal direction, y, is discretized using the standard pseudo-spectral technique with N Chebyshev polynomials. This procedure reduces the system to an M ×M eigen- value problem with M = 4N(2Tm + 1)(2Tn + 1). (3.33) All eigenvalue computations are carried out using the Arnoldi routine as implemented in MATLAB, since only the first few eigenvalues and eigenmodes are relevant in linear stability analyses. 3.4.2 Time-periodic states We now consider the stability of a state which is periodic in time with period T U s(x, y, z, t) = U s(x, y, z, t+ T ). (3.34) The equations governing an infinitesimal perturbation to the state are, as before, equa- tions (3.23). However, with a view to transforming (3.23) into a closed evolution system 54 3. Stability and bifurcations in rotating plane Couette flow and reducing computational memory requirements, we switch from velocity and pressure variables to poloidal-toroidal potentials. Since the velocity is divergence free, we express the velocity as the sum of the poloidal potential φ and the toroidal potential ψ u = ∇×∇× φj +∇× ψj, (3.35) where ∇× denotes the curl operator. The no-slip boundary conditions impose φ = φ′ = ψ = 0, (3.36) on the walls y = ±1. We apply the operators j · ∇ × ∇× and j · ∇× to the linearized Navier-Stokes equation to form two independent equations for φ and ψ ∂t∇ 2∆2φ = ( 1 Re ∇2−U∂x ) ∇2∆2φ+U ′′∂x∆2φ−Ro∂z∆2ψ−j· ( ∇×∇× ( u·∇U s+U s·∇u ) ) , (3.37a) ∂t∆2ψ = ( 1 Re ∇2−U∂x ) ∆2ψ−U ′∂z∆2φ+Ro∂z∆2φ− j · ( ∇× ( u · ∇U s +U s · ∇u ) ) , (3.37b) where ∆2 = ∂2x + ∂ 2 z is the two-dimensional Laplacian, and we have assumed that the mean flow is purely streamwise U 0 = U(y, t)i. (3.38) The velocity perturbation u in system (3.37) is mapped to the poloidal-toroidal poten- tials, ξ = (φ, ψ), using (3.35) u =    ∂2xy −∂z −∆2 0 ∂2yz ∂x    ξ. (3.39) It is convenient to write the system (3.37) more succinctly using operator notation ∂tMξ = L(U s)ξ (3.40) =⇒ ∂tξ = L(U s)ξ, L(U s) = M−1L(U s). (3.41) Since U s is periodic in time so too is the operator L(U s). We now construct the Poincare´ map, P (U s), which is formally given by P (U s) = exp (∫ T 0 L(Us)dt ) . (3.42) 3.4 Floquet stability analysis 55 The Poincare´ map gives the action of system (3.41) on an arbitrary potential vector ξ over one period. i.e. ξ(x, t+ T ) = P (U s)ξ(x, t). (3.43) The eigenvalues µ of P (U s) are the Floquet multipliers over the period µ = eσT , (3.44) from which we determine the exponent σ σ = 1 T log µ. (3.45) En route to computing P (U s), by Floquet’s theorem we have that ξ need not have the same periodicity as U s, so we write ξ(x, y, z, t) = ei(αx+βz)+σt ∞∑ m=−∞ ∞∑ n=−∞ ξˆmn(y, t)e i(mαsx+nβsz) (3.46) as we had for q = (u, p) in equation (3.24). The modes ξˆmn(y, t) are assumed to be periodic in time with the same period as U s. Therefore we can frame (3.41) as a system of infinitely many coupled equations ∂tξˆmn = Amnξˆmn + ∑ j,k F jkmnξˆm−j,n−k, ∀m,n ∈ Z. (3.47) As in the case for steady states, we truncate the system in the periodic directions with −Tm < m < Tm, −Tn < n < Tn. (3.48) and discretize in the wall-normal direction with N Chebyshev polynomials. With the system suitably discretized and truncated, we use the second-order implicit trapezium rule method to advance the equations in time. For y′ = f(y, t), (3.49) the trapezium rule gives ys+1 = ys + h 2 ( f(ys, ts) + f(ys+1, ts+1) ) , (3.50) where h = ts+1 − ts. Thus we have ξs+1 = ξs + h 2 ( Lsξs + Ls+1ξs+1 ) , (3.51) 56 3. Stability and bifurcations in rotating plane Couette flow ( I − h 2 Ls+1 ) ξs+1 = ( I + h 2 Ls ) ξs, (3.52) ξs+1 = ( I − h 2 Ls+1 )−1( I + h 2 Ls ) ξs, (3.53) where I is the identity operator. The trapezium rule is chosen because we found that explicit methods were generally unstable and higher-order methods have greater com- putational memory requirements at each time step. The action of the Poincare´ map is found by time-stepping the equations over one period T . At each time step the base flow U s must be updated. We do this by Fourier series approximation of the base flow at the desired time point. From DNS of U s over one period we save the flowfield at 100 time points, and from these we numerically approximate the temporal Fourier transform U˜ ` = ∫ T 0 U s(t)e i2pi`t T dt, (3.54) so that the flowfield at an arbitrary time t can be expressed by the Fourier series U s(t) = ∞∑ `=−∞ U˜ `e− i2pi`t T . (3.55) A summation truncation of −30 ≤ ` ≤ 30 is used. Comparison of the Fourier series approximation of the flowfield with a previously saved field has an error O(10−11). The advantages of using the poloidal-toroidal potentials rather than primitive variables are that we do not need an extra calculation to update the pressure at each time step, and the problem size is now reduced by half such that the operators are approximated by square matrices of order M = 2N(2Tm + 1)(2Tn + 1). This eases the memory requirements for each stability calculation and allows for higher truncations to be reached. Time-stepping of the linear perturbation equations and the final eigenvalue calculations are performed in MATLAB, while DNS of U s required for the temporal Fourier transform is carried out by Channelflow (Gibson (2008)). 3.4.3 Global stability Having now outlined how to find the stability of a spatially-periodic stateU s(αs, βs, Re,Ro) to an infinitesimal perturbation, we can thus produce a stability map by finding the ex- ponents σ(α, β, αs, βs, Re,Ro) over all de-tuning parameters α and β. In this way, we can determine whether a state is globally stable. Owing to the periodicity of U s there is some saving in the range of α and β at which σ must be computed to find the global stability properties. In what follows, we discuss the relevant (α, β) stability domain for 3.4 Floquet stability analysis 57 a steady base flow, however the argument for unsteady states is analogous. Recall (3.25) Amnqˆmn + ∑ j,k F jkmnqˆm−j,n−k = σBqˆmn, ∀m,n ∈ Z. (3.56) We can consider the operators Amn and F jkmn as functions of α and β with αs, βs, Re and Ro fixed. Clearly Amn(α + αs, β + βs) = Am+1,n+1(α, β), (3.57) and F jkmn(α + αs, β + βs) = F jk m+1,n+1(α, β). (3.58) Since the system (3.25) is completed by summation over m,n ∈ Z, we have that σ(α + αs, β + βs) = σ(α, β). (3.59) By induction we can generalize this to σ(α + pαs, β + qβs) = σ(α, β), p, q ∈ Z. (3.60) Furthermore, we have that Amn(−α,−β) = A−m,−n(α, β) ∗, (3.61) and F jkmn(−α,−β) = F −j,−k −m,−n(α, β) ∗, (3.62) where .∗ denotes the (element-wise) complex conjugate, and we have used U jk(y)∗ = ∫ 2pi α 0 ∫ 2pi β 0 U s(x, y, z)ei(jαsx+kβsz) dxdz = U−j,−k. (3.63) Suppose we have found σ = σ(α, β) which satisfies Amn(α, β)qˆmn + ∑ j,k F jkmn(α, β)qˆm−j,n−k = σBqˆmn, ∀m,n ∈ Z (3.64) and we wish to find σ˜ = σ(−α,−β) Amn(−α,−β)q˜mn + ∑ j,k F jkmn(−α,−β)q˜m−j,n−k = σ˜Bq˜mn, ∀m,n ∈ Z. (3.65) Then using (3.61) and (3.62) we have A−m,−n(α, β) ∗q˜mn + ∑ j,k F−j,−k−m,−n(α, β) ∗q˜m−j,n−k = σ˜Bq˜mn, ∀m,n ∈ Z. (3.66) 58 3. Stability and bifurcations in rotating plane Couette flow Taking the conjugate of (3.64) gives Amn(α, β) ∗qˆ∗mn + ∑ j,k F jkmn(α, β) ∗qˆ∗m−j,n−k = σ ∗Bqˆ∗mn, ∀m,n ∈ Z (3.67) and since these equations must be satisfied for m,n, j, k ∈ Z we conclude that σ˜ = σ∗ =⇒ σ(−α,−β) = σ(α, β)∗. (3.68) Now, suppose we wish to find σ = σ(rαs, sβs), with 1 2 < r, s < 1. (3.69) Then σ = σ(−rαs,−sβs) ∗ = σ(−rαs + αs,−sβs + βs) ∗ (3.70) = σ((1− r)αs, (1− s)βs) ∗. (3.71) Hence the global stability properties of U s can be determined for de-tuning wavenumbers in the domain 0 ≤ α ≤ αs 2 , 0 ≤ β ≤ βs 2 . (3.72) If U s is periodic in only one direction and independent of the other, such as the Taylor vortex solution in §3.5.1, then the stability domain becomes the semi-infinite strip 0 ≤ α <∞, 0 ≤ β ≤ βs 2 . (3.73) 3.4.4 Convergence By variation of the truncation parameters Tm and Tn and the Chebyshev discretization N we can assess the convergence of the eigenvalues σ for a given solution U s. Our computational capability puts an upper limit on the truncation and discretization we can use. Running on a machine with 24GB RAM, we find that a matrix size M ≤ 17000 can be handled. For a Taylor vortex solution we can truncate to a high degree of accuracy, since the state is streamwise independent and truncation is required in the spanwise direction only. For fully three-dimensional flows we cannot truncate to the same degree, though we still find good convergence of the leading eigenmode. In Tables 3.1, 3.2, 3.3 and 3.4 we list the leading eigenvalues from stability calculations on a Taylor vortex, a wavy vortex (see WVF §3.7.1) and an oscillatory wavy vortex (see oWVF §3.7.3) at various truncations. 3.4 Floquet stability analysis 59 N 17 25 49 Tn 6 0.025574 0.025573 0.025573 7 0.025726 0.025725 0.025725 8 0.025687 0.025686 0.025686 9 0.025692 0.025691 0.025691 10 0.025691 0.025690 0.025690 11 0.025691 0.025690 0.025690 12 0.025691 0.025690 0.025690 13 0.025691 0.025690 0.025690 14 0.025691 0.025690 0.025690 15 0.025691 0.025690 0.025690 16 0.025691 0.025690 0.025690 17 0.025691 0.025690 0.025690 18 0.025691 0.025690 0.025690 19 0.025691 0.025690 0.025690 20 0.025691 0.025690 0.025690 Table 3.1: Values of the leading eigenvalue σ for a Taylor vortex solution with (α, β, βs, Re,Ro) = (0.5, 0, 0.5236, 2.513, 100, 0.1). Tn 4 5 6 7 8 9 Tm 4 0.023226 0.021836 0.021661 0.021647 0.021648 0.021648 5 0.022822 0.021309 0.021169 0.021164 0.021164 0.021164 6 0.022805 0.021334 0.021163 0.021163 0.021163 0.021163 7 0.022830 0.021351 0.021187 0.021184 0.021185 - 8 0.022836 0.021356 0.021191 0.021189 - - 9 0.022834 0.021356 0.021191 - - - Table 3.2: Values of the leading eigenvalue σ for a wavy vortex solution with N = 17 and (α, β, αs, βs, Re,Ro) = (0, 1.2565, 0.5236, 2.513, 100, 0.1). Tn 4 5 6 7 8 Tm 4 0.023229 0.021838 0.021663 0.021649 0.021650 5 0.022830 0.021317 0.021178 0.021172 - 6 0.022817 0.021345 0.021174 - - 7 0.022844 0.021364 - - - 8 0.022851 - - - - Table 3.3: Values of the leading eigenvalue σ for a wavy vortex solution with N = 25 and (α, β, αs, βs, Re,Ro) = (0, 1.2565, 0.5236, 2.513, 100, 0.1). 60 3. Stability and bifurcations in rotating plane Couette flow Tn 3 4 5 6 Tm 3 0.002903 0.006839 0.006993 0.006992 4 0.002937 0.006840 0.007024 0.007016 5 0.002933 0.006837 0.007021 0.007014 6 0.002933 0.006838 0.007021 0.007014 Table 3.4: Values of the leading exponent σ for oscillatory wavy vortex flow (α, β, αs, βs, Re,Ro) = (0.1, 0.375, 0.8, 3, 100, 0.55), N = 17, h = T500 = 0.03844. 3.4.5 Validation Predictions from the stability calculations are validated against the Navier-Stokes solver SIMSON (Chevalier et al. (2007)). We compare the growth rates predicted by the sta- bility analysis with growth rates of a divergence-free, random disturbance comprised of Stokes modes which is added to the solution and integrated forwards. A domain-doubling technique is used to test growth rates of a subharmonic perturbation; the solution is repeated in the periodic direction such that subharmonic disturbances can be superim- posed on to the original solution. In Figure 3.4 we plot the evolution of the perturbation energy-density E = E˜ ρ = 1 2 ∫ D u2 + v2 + w2dV, (3.74) of a given mode, where ρ is the fluid density, D is the periodic domain and E˜ is the physical kinetic energy. We can then compare the growth rate of the perturbation as observed from simulations with the prediction made by the stability analysis. As can be seen in Figure 3.4, excellent agreement is obtained in each case. 3.4 Floquet stability analysis 61 Figure 3.4: (a) TVF with (α, β, βs, Re,Ro) = (0.5, 0, 2, 100, 0.11), (b) WVF with (α, β, αs, βs, Re,Ro) = (0, 1, 0.5, 2, 100, 0.11) and (c) oWVF with (α, β, αs, βs, Re,Ro) = (0.1, 0.375, 0.8, 3, 100, 0.55). In each panel, the blue curve is the perturbation evolution from DNS of the Navier-Stokes equations and the red dashed line is the Floquet theory prediction based on the leading growth rate σr. 62 3. Stability and bifurcations in rotating plane Couette flow 3.5 Secondary states In §§3.5.1 to 3.5.3 we characterize the secondary states in supercritical RPCF. We inves- tigate the effect of rotation on the solutions by making bifurcation diagrams in (Ro,Ecf )- sapce for Re = 100. In each bifurcation diagram we include information on the harmonic stability of the state i.e. the stability of the solution to a perturbation with α = αs and β = βs as defined in (3.24). Unstable states are represented by a dashed line while stable states are indicated with an unbroken line. 3.5.1 Taylor vortex flow A bifurcation initiated by the leading streamwise vortex instability eigenmode results in Taylor vortex flow (TVF). TVF consists of a pair of streamwise independent counter- rotating vortices and counter-propagating streaks. The streamwise streak component of the velocity field changes appearance as the geometry of the vortices varies from long to short spanwise wavelength, with the streaks become more rounded as the wavelength is decreased. The ellipticity and orientation of the vortical part of the flowfield also changes as the wavenumber increases. These changes in the flowfields can be seen in Figure 3.5. We plot bifurcations of TVF from the primary state (Ecf = 0) in Figure 3.7. The Figure includes curves of the growth rate of the linear mode responsible for each bifurcation from which we can see that each bifurcation is transcritical, with TVF appearing only when the linear mode has a positive growth rate. TVF only becomes harmonically unstable for β = 4, in the range 0.084 < Ro < 0.91. The Taylor vortices possess mirror symmetry Z, where Z[u, v, w](x, y, z) = [u, v,−w](x, y,−z). (3.75) This symmetry feature is apparent from the flowfields in Figures 3.5, with the axis of reflection being the mid-point of the spanwise domain. 3.5.2 Second Taylor vortex flow As noted by Nagata (2013) the primary flow loses stability to a second streamwise vortex eigenmode, prompting a bifurcation to a second Taylor vortex flow (TVF2). TVF2 has a double layered vortical structure, with a pair of counter-rotating vortices aligned in the wall-normal direction, alongside another pair of counter-rotating vortices in the spanwise direction, as can be seen in Figure 3.6. The mirror symmetry Z holds in TVF2, with the axis of reflection the mid-point of the spanwise domain in our Figures. Bifurcations of TVF2 in Ro are plotted in Figure 3.8. The trajectories of the solutions are reminiscient of those of the first TVF in the (Ro,Ecf ) parameter space, and each bifurcation is similarly transcritial. TVF2 is unstable at each point in our bifurcation diagrams. 3.5 Secondary states 63 Figure 3.5: Taylor vortex flowfields at various spanwise wavelengths and rotation num- bers. Re = 100 in all cases. The orientation of the vortices and the structure of the streamwise streaks are subject to change as βs and Ro are varied. Figure 3.6: Flowfields of TVF2 at Re = 100 for βs = 2, Ro = 0.5 and β = 4, Ro = 0.3. 64 3. Stability and bifurcations in rotating plane Couette flow Figure 3.7: (a) Taylor vortex bifurcations for various spanwise wavenumbers. Re = 100 in each case. In each case the growth rate of the linear mode responsible for the bifurcation is plotted as a dashed line of the colour corresponding to the TVF trajectory. The linear growth rates are scaled ωrRe . (b) Blow-up of the low Ro region where the loss of stability of TVF with βs = 4, and the emergence of each TVF with the linear mode, can be clearly seen. Figure 3.8: Bifurcations of TVF2 at βs = 2 and β = 4, with Re = 100 in both cases. The linear mode responsible for the bifurcation is included in each panel. Both solutions are unstable along their entire trajectories. 3.5 Secondary states 65 3.5.3 Oblique vortex flow In addition to losing stability to streamwise independent perturbations, linear stability theory tells us that the primary state can also be de-stabilized by oblique waves. As α is increased from zero, we generically find that the two unstable oblique modes move closer together until they merge to become a pair of modes with non-zero frequency, i.e. with ωi 6= 0. An example of this is shown in Figure 3.9. The two zero frequency unstable modes develop into steady, nonlinear oblique vortex flow (OVF). OVF, like TVF, takes the form of pairs of counter-rotating vortices, but unlike TVF the counter-propagating streaks are oblique; inclined at an angle θ = tan αβ away from the streamwise direction. We investigate the trajectories of OVF in the (Ro,Ecf ) subspace as it bifurcates from the linear shear profile, for the case (β,Re) = (2, 100). We find that for α . 0.225 there are upper and lower solutions which bifurcate from each unstable eigenmode. We name these OVF1 and OVF2 respectively. In Figure 3.10 we see that each solution forms two distinct curves in (Ro,Ecf ). OVF1 is harmonically stable for a large portion of its bifurcation, whereas the OVF2 is always unstable. We note that OVF1 is unstable for Ro . 0.1, though we do not pursue tertiary bifurcations here. As α is increased the two solutions merge to form a closed trajectory in (Ro,Ecf ) and a separate, disjointed set of high and low rotation solutions. This is shown in Figure 3.11. A depiction of the flowfields of two OVF solutions is given in Figure 3.12. OVF has a Taylor-Couette flow analogue in spiral vortex flow (SPI), which has recently been studied by Hoffmann et al. (2009), Altmeyer et al. (2010) and Deguchi & Altmeyer (2013), though in these cases the spiral vortices are travelling waves, rather than steady flows as considered here. Figure 3.9: Linear stability modes (β,Re,Ro) = (2, 100, 0.5). The blue corresponds to the first unstable mode, red is the second mode and black is the merged mode with ωi 6= 0. 66 3. Stability and bifurcations in rotating plane Couette flow Figure 3.10: OVF1 and OVF2 for (αs, βs, Re) = (0.2, 2, 100) and (αs, βs, Re) = (0.225, 2, 100). The growth rates of the linear modes which initiate bifurcation from the primary state are indicated by dotted lines. Figure 3.11: Merged OVF solutions for (αs, βs, Re) = (0.23125, 2, 100) and (αs, βs, Re) = (0.2375, 2, 100). The curves are coloured differently to Figure 3.10 to indicate that merg- ing has taken place. The outer and inner solutions of each curve can still be considered OVF1 and OVF2, respectively. The growth rates of the linear modes which initiate bifurcation from the primary state are indicated by dotted lines. 3.5 Secondary states 67 Figure 3.12: Flowfields of the oblique vortex flows OVF1 and OVF2 for (αs, βs, Re,Ro) = (0.225, 2, 100, 0.5). 3.5.4 Summary of secondary states The secondary states which we have discussed fall into two categories: streamwise inde- pendent and streamwise periodic states. Two streamwise independent states, TVF and TVF2 have been introduced and their bifurcations in Ro have been calculated. Stream- wise periodic states OVF1 and OVF2 bifurcate as distinct solutions at low streamwise wavenumber αs before merging as αs is increased. It is to be expected that any transition from the primary Couette flow will pass through one of these secondary states. 68 3. Stability and bifurcations in rotating plane Couette flow 3.6 Stability of Taylor vortex flow In this section we determine the global stability of Taylor vortex flow. In particular we focus on the stability of the first Taylor vortex solution, labelled TVF, since the most unstable primary instability typically develops into this solution. Therefore it is useful to characterize the secondary instabilities that emerge when TVF loses stability. Let βs denote the spanwise wavenumber of a Taylor vortex pair. By varying the detuning parameters α and β from equation (3.24), we can determine the global stability of a Taylor vortex solution. It is typically found that the least stable perturbations are either streamwise periodic and tuned to the fundamental, i.e. (α, β) = (α, 0), or of Eckhaus type with (α, β) = (0, β) (Eckhaus (1965)). From §3.4.3 we have that the wavenumber domain which spans the full stability properties is the semi-infinite strip 0 ≤ α <∞, 0 ≤ β ≤ βs 2 . (3.76) Figure 3.13: Global stability of TVF with (a) (βs, Re,Ro) = (1.5, 100, 0.7) and (b) (βs, Re,Ro) = (4, 100, 0.3), with perturbation growth rates σr ≥ 0 coloured. In (a) the least stable perturbation is streamwise periodic and tuned to the fundamental spanwise wavenumber. In (b) the least stable perturbation is harmonic, with α = β = 0. In Figure 3.13 we plot the global stability properties of two typical TVF profiles. To clearly highlight the unstable parameter regions, we plot the intersection of the growth rate σr with the neutral plane σr = 0, such that all stable regions are coloured dark blue regardless of the value of the decay rate. In each of the cases considered we find that the least stable perturbation is tuned to the fundamental, i.e. has spanwise detuning wavenumber β = 0. As such, we choose three Reynolds numbers Re = 50, 100 and 250, and calculate the streamwise secondary instabilities at various βs. The chosen Reynolds numbers roughly correspond to the start, the middle and the end of the transition region, 3.6 Stability of Taylor vortex flow 69 as observed experimentally in Tsukahara et al. (2010). The spanwise wavenumbers of the Taylor vortices, βs, are chosen such that they are spread across the linearly unstable region of the primary flow, as shown in Figure 3.14. Figure 3.14: Colour-plots of the primary instability growth rate ωr at Re = 50, 100, 250. The thick black lines mark the locations of βs for the Taylor vortices considered at each Re. To create each stability map in this section, minimum grids of 100×200 points in (α,Ro) or (β,Ro) are used. The number of points used in each map varies, as we refine in the grid to better resolve cases with complicated instability regions. The colouring of the map is then interpolated to give a smooth visual effect. An example is given in Figure 3.15. 70 3. Stability and bifurcations in rotating plane Couette flow Figure 3.15: The colouring of the plot on the left is interpolated to give the plot on the right. 3.6.1 Streamwise instability maps In Figures 3.16, 3.21 and 3.26 we plot the streamwise instabilities which affect Taylor vortices of various spanwise wavelengths for Re = 50, Re = 100 and Re = 250. In Figure 3.16 where we have Re = 50, for βs = 1, 1.5, 2 and 2.5 we see a wedge shaped region of instability for 0 . Ro . 0.2 and 0 . α . 0.6. This is known as wavy vortex instability (Davey et al. (1968), Nagata (1986, 1998)), and we plot examples of its mode structure and eigenvalues of the stability operator in a case with wavy instability, in Figures 3.17 and 3.18. The wavy instability causes Taylor vortices to lose stability to structures elon- gated in the streamwise direction, with modulated high-velocity/low-velocity streaks. Throughout the wedge regions in Figure 3.16, the stability operator has one unstable eigenmode with σi = 0. As the spanwise wavenumber of the Taylor vortices increases to βs = 3 the wavy instability is no longer present, and the vortices become stable to all streamwise perturbations. For βs = 1 a region of instability occurs at high rotation parameters, 0.6 . Ro . 1 and 0.5 . α . 1.5. This is known as twist vortex instabil- ity (Weisshaar et al. (1991)). There are two unstable modes, the wavy twist and twist vortex modes, which compete to de-stabilize the Taylor vortices. We show examples of each mode and the stability eigenvalues in Figures 3.19 and 3.20. We find that the wavy twist mode has a larger growth rate compared to the twist mode for Re = 50. The twist instabilities are more sensitive to the spanwise wavelength of the TVF than the wavy instability, as they are not observed for βs ≥ 1.5. An increase in Re has the generic effect of increasing the range and strength of insta- bilities across the range of TVF, as can be seen for Re = 100 Figure 3.21. The wavy vortex instability maintains the wedge shape observed for Re = 50, though the region 3.6 Stability of Taylor vortex flow 71 stretches to α > 1 for βs = 2.6, 3, and is confined to a smaller Ro range in comparison to Re = 50. For 0.4 . Ro . 1 and βs = 1, 1.5 we see a region of twist instabilities. Whichever mode is least stable is dependent on α and Ro, as it switches between the wavy twist and twist mode in (α,Ro)-space. Figure 3.22 gives an example of this. A mid-Ro instability emerges for βs = 2.5 and strengthens as βs increases to 3.5. The stability operator has a pair of complex conjugate modes in this case so we term this the oscillatory wavy instability, since we expect unstable perturbations to oscillate in time as their amplitude grows. The instability first occurs in the range 0.5 . α . 1, giving the modes a streamwise elongated appearance similar to the low-Ro wavy instabil- ity modes. This instability was noted by Nagata (1988), whose attention was restricted to Taylor vortices with spanwise wavelength βs = 1.5585, in which case the instability first emerges for Re ≈ 137.5. We plot examples of the mode and the eigenvalues of the stability operator in Figures 3.23 and 3.25. As βs increases the range of α over which the oscillatory instability is present decreases until the region merges with a harmonic instability region in βs = 4. The harmonic instability persists for βs = 4.5. We find that the unstable mode in this case leans against the Taylor vortex streaks causing them to be tilted in the spanwise direction, therefore we call this the skew instability. The mode is depicted in Figure 3.24, and its stability operator eigenvalues are shown in Figure 3.25. Note in Figure 3.25 for the operators with skew instability, there is an eigenvalue close to zero, which we will call σ0. When calculating the stability of a periodic solution to harmonic perturbations, we expect zero eigenvalues to appear as an indication that the state is indeed a solution to the governing equations. For our calculations, we compute σ0r < 10 −6, which we interpret as an indication of the accuracy of our code. The stability maps for Re = 250 in Figure 3.26 show instabilities covering greater pro- portions of the (α,Ro) domain than was covered for Re = 50, 100. The instabilities themselves, however, are the same as those which appeared for lower Re. For βs = 1 a region of wavy twist instability dominates (α,Ro)-space, with a small region of wavy instability at low Ro. The wedge shaped region of wavy instability from Re = 50, 100 is present again for βs = 2, at low Ro, and a thin strip of oscillatory instability appears for Ro ≈ 0.02 and 1 . α . 3. For βs = 2 at higher Ro, there are three distinct regions of instability. Two oscillatory instabilities account for the two lower α regions, while a combination of steady wavy twist and twist instabilities dominate for the higher α region. For βs = 3 and 4 the oscillatory instability dominates over 0.1 . Ro . 0.9, with a wedge region of steady wavy instability persisting for Ro . 0.1. For βs ≥ 5 the harmonic skew instability dominates, and instabilities have increasingly longer streamwise wavelengths. Unlike the cases for Re = 50, 100, there are no regions other than close to the linear stability boundary where the Taylor vortices for Re = 250 are stable to streamwise per- turbations. 72 3. Stability and bifurcations in rotating plane Couette flow Figure 3.16: Streamwise instability for Re = 50. σr ≥ 0 is coloured in (α,Ro)-space. 3.6 Stability of Taylor vortex flow 73 Figure 3.17: Wavy instability mode for (α, βs, Re,Ro) = (0.4, 1, 50, 0.1). On the top left we have TVF and on the right is the wavy mode which de-stabilizes it. Beneath we plot a combination of the two modes. Figure 3.18: Eigenvalues of the stability operator with a wavy instability for (α, βs, Re,Ro) = (0.4, 1, 50, 0.1021) (left) and (α, βs, Re,Ro) = (0.5, 2, 50, 0.1261) (right). 74 3. Stability and bifurcations in rotating plane Couette flow Figure 3.19: Wavy twist and twist instability modes for (α, βs, Re,Ro) = (1, 1, 50, 0.7606). The top row have an unstable TVF and the twist perturbations which de-stablize the vortex. The bottom row depicts the composition of the TVF and each instability mode. Figure 3.20: Eigenvalues of the stability operator with a wavy instability for (α, βs, Re,Ro) = (1, 1, 50, 0.7606). The twist eigenvalue is circled in red so as to dis- tinguish it from the wavy twist eigenvalue, since the two are very close. On the right we plot the growth rates σr of the wavy twist (blue) and the twist (red) modes against α, for Ro = 0.7606. In this instance, the wavy twist mode dominates for all streamwise wavenumbers. 3.6 Stability of Taylor vortex flow 75 Figure 3.21: Streamwise instability for Re = 100. σr ≥ 0 is coloured in (α,Ro)-space. 76 3. Stability and bifurcations in rotating plane Couette flow Figure 3.22: Wavy twist (blue) and twist (red) instability modes competing in Ro. (α, βs, Re) = (1.5, 1.5, 100). Unlike the case for Re = 50 in Figure 3.20, the least stable modes changes from twist to wavy twist as Ro increases. Figure 3.23: From top to bottom we plot TVF; the oscillatory wavy perturbation mode; and the combination of the two flowfields. The parameters here are (α, βs, Re,Ro) = (0.5, 3, 100, 0.5). 3.6 Stability of Taylor vortex flow 77 Figure 3.24: From left to right we have TVF; the skew mode; and the combination of TVF and the skew mode. We have chosen the name “skew” instability as the combination flowfield resembles a Taylor vortex pair which is skewed in the (y, z)-plane. The parameter values are (βs, Ro,Re) = (4, 0.2, 100). Figure 3.25: Stability operators with unstable oscillatory wavy modes (top) and skew modes (bottom). Re = 100 in all cases. 78 3. Stability and bifurcations in rotating plane Couette flow Figure 3.26: Streamwise instability for Re = 250. σr ≥ 0 is coloured in (α,Ro)-space. 3.6 Stability of Taylor vortex flow 79 3.6.2 Eckhaus instability maps In Figures 3.27, 3.28 and 3.30 we plot the streamwise independent instabilities that af- fect Taylor vortices for Re = 50, 100 and 250. An instability caused by a perturbation which is periodic in the same direction as the base state is called Eckhaus instability (Eckhaus (1965)). For Re = 50, 100, Eckhaus instability is limited to TVF with higher βs wavenumber so in Figures 3.27 and 3.28 we include maps for the two largest βs consid- ered. In each Figure, low Ro and high Ro tongues of subharmonic instability appear for the Taylor vortices. These modes have σi = 0, and initiate a bifurcation towards TVF with lower βs as is outlined in Figure 3.29 for a case with Re = 50. In Figure 3.28, for Re = 100, these tongues are restricted to a smaller region of the parameter space, and the skew instability emerges between 0.1 . Ro . 0.9. As was the case with streamwise instabilities, for Re = 250 the Taylor vortices are in general much less stable than was the case for lower Reynolds numbers. Only βs = 1 is Eckhaus stable, all other vortices being de-stabilized by streamwise independent perturbations. We find that in most cases the maximum Eckhaus instability is in fact harmonic, being either the skew instability or an oscillatory version of it. Vortices with higher spanwise wavenumber have stronger Eckhaus instability than those with lower wavenumber, as evidenced by increase in max- imum growth rate as βs is increased in Figure 3.30. The low and high Ro tongues of instability present in the Re = 50, 100 Eckhaus instability maps appear only at βs = 7, 8 for Re = 250, with the tongues confined to Ro very close to the upper and lower linear stability boundaries. Figure 3.27: Eckhaus instability for Re = 50 σr ≥ 0 is coloured in (β,Ro)-space. 80 3. Stability and bifurcations in rotating plane Couette flow Figure 3.28: Eckhaus instability for Re = 100 σr ≥ 0 is coloured in (β,Ro)-space. Figure 3.29: Subharmonic Eckhaus modes at (β, βs, Re,Ro) = (1.75, 3.5, 50, 0.859). At the top we plot TVF with βs = 3 over the doubled wavelength 4pi3 . On the top right we plot uβ, the velocity perturbation about the Taylor vortex, corresponding to Eckhaus mode β = 1.75. u−β corresponds to the conjugate mode, plotted in the bottom left panel. In the final panel we plot the superposition of the two perturbations, which resemble a Taylor vortex flowfield with βs = 1.75. 3.6 Stability of Taylor vortex flow 81 Figure 3.30: Eckhaus instability for Re = 250. σr ≥ 0 is coloured in (β,Ro)-space. 82 3. Stability and bifurcations in rotating plane Couette flow 3.6.3 Summary of TVF stability properties The stability of TVF has been shown to be strongly dependent on the spanwise wavenum- ber βs of the Taylor vortices, the rotation number Ro and the Reynolds number Re. We have mapped the streamwise instabilities across all Ro for Re = 50, 100, 250 and across a set of βs for which there exist Taylor vortices, in Figures 3.16, 3.21 and 3.26, and the Eckhaus instabilities across the same parameters in Figures 3.27, 3.28 and 3.30. Care has been taken to analyze each stability mode, to give a characterization to the instabilities which cause TVF to break down at given points in parameter space. The instabilities cover a larger region of (α,Ro) space as Re is increased, offering some promise that exper- imentally observed transition phenomena can be explained via bifurcation of structures from TVF. Our analysis has focused on purely streamwise and purely spanwise insta- bilities, so it is of course possible that there exist oblique instabilities which de-stabilize the Taylor vortices. However, as shown in Figure 3.13, where we find oblique instabil- ities in our calculations, they are often not as strong as purely streamwise or spanwise modes. We take this as justification in studying the reduced range of instabilities, in or- der to give a clearer view of the parameter dependence of the stability properties of TVF. By giving an indication of where in the (α, β, βs, Re,Ro) parameter space Taylor vortices will break down, our stability maps provide an ideal starting point for investigating the tertiary bifurcations which branch from secondary TVF. 3.7 Bifurcations of tertiary states 83 3.7 Bifurcations of tertiary states We define a tertiary state to be a structure which bifurcates directly from a secondary state, but not the primary state. In this analysis we identify a range of the tertiary states which bifurcate from Taylor vortex flow at transitional Reynolds numbers. To this end, we restrict our attention to Re = 100, since §3.6 indicates a range of instabilities present for this Re. To assess the effects of rotation on the tertiary states we track their trajecto- ries in (Ro,Ecf )-space and create bifurcation diagrams in §§3.7.5, 3.7.6 and 3.7.7 We use the Newton-Krylov-hookstep algorithm outlined in §3.3 to find the tertiary states and an adaptive quadratic continuation technique to continue them in solution space. In each bifurcation diagram, the harmonic stability or instability of each structure is indicated by a solid or dashed curve respectively, as was our convention in §3.5. In addition we perform a global stability analysis on a representative set of the tertiary solutions which appear in our bifurcation diagrams. An extensive tertiary stability anal- ysis would require examination of each solution at all the wavenumbers at which it exists, as was done for TVF in §3.6, which would be computationally expensive. Therefore we simply select cases of each tertiary flow to analyze. The tertiary states on which we per- form global stability analyses are marked by dots where they appear in the bifurcation diagrams. Well converged eigenfunctions require a greater truncation than is needed to compute well converged eigenvalues, so we do not include plots of the eigenfunctions in this thesis. In order to produce coherent maps of the instabilities affecting each solution an extra layer of interpolation is added to the procedure for the maps in §3.6. For a steady solu- tion, stability values are calculated at 100 data points in the (α, β) domain. These data points are then interpolated, using a cubic spline algorithm, onto a grid of over 1000 points before we create the colour-plot. For time-periodic states, since stability values are more computationally expensive to compute, cubic spline interpolation is used over a grid of 20 data points. As in §3.6, we plot the growth rates σr with the intersection of the neutral plane, σr = 0, such that the regions of instability can be clearly seen. We begin with an introduction to each of the tertiary states in §§3.7.1, 3.7.2, 3.7.3 and 3.7.4. 84 3. Stability and bifurcations in rotating plane Couette flow 3.7.1 Wavy vortex flow Wavy vortex flow (WVF) is a steady structure which develops when Taylor vortices lose stability to wavy vortex modes. Wavy vortices are a well-known feature of Taylor- Couette flow (see for example Davey et al. (1968), Andereck et al. (1986), Nagata (1986), Koschmieder (1993)), and were discussed by Nagata (1998) in the context of rotating plane Couette flow. Wavy vortices lose the mirror symmetry Z of Taylor vortices, but they hold shift-reflect and shift-rotation symmetries S and Ω. An example of the streamwise and spanwise variation of the flowfield is given in Figure 3.31. The flowfields closely resemble the superposition of the linear wavy mode and TVF from Figure 3.17. S[u, v, w](x, y, z) = [u, v,−w](x+ pi α , y,−z + pi β ) (3.77a) Ω[u, v, w](x, y, z) = [−u,−v, w](−x,−y, z + pi β ) (3.77b) Figure 3.31: Wavy vortex flows (αs, βs, Re,Ro) = (0.7, 1.5, 100, 0.05) and (αs, βs, Re,Ro) = (0.8, 3, 100, 0.12). 3.7.2 Twist and wavy twist vortex flow Twist (TWI) and wavy twist (wTWI) vortex flows emerge from the twist and wavy twist instabilities of Taylor vortex flow (Weisshaar et al. (1991), Antonijoan & Sa´nchez (2000)). Examples of the flowfields are given in Figure 3.32 where, like WVF, the mid- plane flowfields of the twist flows closely resemble their linear counterparts from Figure 3.19. wTWI has the symmetries S and Ω loses mirror symmetry Z, whereas TWI retains the mirror symmetry of TVF. 3.7 Bifurcations of tertiary states 85 Figure 3.32: Twist (left) and wavy twist (right) vortices for (αs, βs, Re,Ro) = (1, 1, 100, 0.8). 3.7.3 Oscillatory wavy vortex flow The complex conjugate eigenvalues of the oscillatory wavy instability instigate a Hopf bifurcation from Taylor vortex flow towards a periodic orbit of the Navier-Stokes equa- tions. To our knowledge, calculations of the oscillatory flow bifurcating directly from TVF has not been previously reported in the literature for Taylor-Couette flow, circular Couette flow or rotating plane Couette flow, though the oscillatory instability has been known since Nagata (1988). However, Hopf bifurcations from tertiary states have been reported, such as those discussed in Nagata & Kawahara (2004), and oscillatory flows in plane Couette flow have also been found (Kawahara & Kida (2001), Viswanath (2007), Kreilos & Eckhardt (2012)). The structure has a remarkably simple sinusoidal oscillatory amplitude, which is manifest in its Ecf evolution in Figure 3.33. We have named it the oscillatory wavy vortex flow (oWVF) as its flowfield resembles that of a steady wavy vortex which oscillates in time, as can be seen in Figure 3.34. Its structure is much sim- pler than the periodic orbits found in non-rotating plane Couette flow by, for example, Kawahara & Kida (2001) and Viswanath (2007), though its Ecf evolution is reminiscent of the 1 orbit found by Kreilos & Eckhardt (2012). oWVF holds the same symmetries as WVF and wTWI: shift-reflection S and shift-rotation Ω. Figure 3.33: Perturbed TVF evolving towards oWVF for (αs, βs, Re,Ro) = (0.6, 2.4, 100, 0.5). 86 3. Stability and bifurcations in rotating plane Couette flow Figure 3.34: Oscillatory wavy vortex flow for (αs, βs, Re,Ro) = (0.6, 2.6, 100, 0.5), pic- tured at four time points throughout one period of its oscillation. The flowfield has a wavy structure and the symmetries S and Ω are satisfied throughout its orbit. 3.7 Bifurcations of tertiary states 87 3.7.4 Skewed vortex flow Skewed vortex flow (SVF) is the steady, streamwise independent flow which bifurcates from Taylor vortices via the skew instability. SVF loses mirror symmetry Z, like the previous tertiary states excluding TWI. The spanwise skewness of the vortices is readily apparent in velocity field pictured in Figure 3.35. The flowfield is quite different from the linear mode and TVF combination of Figure 3.24, suggesting that a greater number of spanwise harmonics are excited during nonlinear saturation of the linear mode. We are unaware of SVF being reported by other authors. Figure 3.35: Skewed vortex flow for (βs, Re,Ro) = (4, 100, 0.9). 88 3. Stability and bifurcations in rotating plane Couette flow 3.7.5 Bifurcation of twist, wavy twist and wavy vortex flows Having classified the instabilities that cause Taylor vortices to break down and the ter- tiary states that they develop into, we can trace their trajectories in the (Ro,Ecf ) sub- space, with αs, βs and Re fixed. In Figure 3.36 we show how WVF and wTWI bifurcate from TVF with βs = 1.5. WVF is born when TVF loses stability to the wavy vortex mode and when TVF regains stability WVF dies, thus bifurcation to WVF is considered tran- scritical in this case. The stability of this WVF solution changes twice, with instability to oscillatory perturbations in the regions 0.0154 < Ro < 0.02 and 0.0246 < Ro < 0.041. Global stability analysis on a harmonically stable and unstable WVF solution, shown in Figure 3.40, indicate that both states de-stabilize to a subharmonic instability with maximum growth rate at approximately (α, β) = (0.35, 0.135). wTWI emerges, in Figure 3.36, in a saddle-node bifurcation at Ro = 0.604 with an upper and lower branch of solutions in the range 0.604 < Ro < 0.756, and both solu- tions are harmonically unstable to complex conjugate modes. The lower branch becomes harmonically stable at Ro = 0.635 meaning that there are two stable attractors over 0.635 < Ro < 0.756 since TVF is also stable, while the upper branch solution is un- stable throughout the bifurcation. At Ro = 0.756, TVF loses stability and the lower branch wTWI becomes the only attractor, while the upper branch merges with TVF. The structure of upper branch wTWI differs from the structure of its lower counter- part, as highlighted in the (x, z)-projections of two example flowfields in Figure 3.38. The lower branch solution has a sinusoidal arrangement of vortex pairs embedded be- tween streamwise streaks, whereas the upper branch solution consists of three pairs of vortices embedded between each streak, along one streamwise wavelength. Each group of three, in the upper branch solution, contains two strong and one weak vortex pairs, with the weak and strong pairings aligning either side of the streak. In Figure 3.39 we plot the global stability properties of the harmonically stable lower branch solution and the harmonically unstable upper branch solution indicated by dots in Figure 3.36. The lower branch solution is globally stable, whilst the upper branch solution is unstable to a perturbation of any (α, β) wavenumber pairing. However, the upper branch solution is maximally unstable to harmonic perturbations. In Figure 3.37 we plot bifurcations for αs = 1.5. TVF remains stable for Ro < 0.55, above which stability is lost to a twist mode. We then see TWI and wTWI bifurcate transcrit- ically from TVF. TWI is harmonically stable in the range 0.55 < Ro < 0.618, beyond which it loses stability to a stationary mode. We investigate the global stability of the two solutions indicated in the bifurcation diagram in Figure 3.41. Harmonically stable TWI 3.7 Bifurcations of tertiary states 89 loses stability to a subharmonic stationary mode with maximum (α, β) = (0.37, 0.23), while the harmonically unstable TWI is globally unstable with maximum to an oscil- latory subharmonic mode at (α, β) = (0.72, 0.48). TVF becomes stable once more for Ro > 0.852 and the tertiary states collapse. In Figures 3.42 and 3.43 we investigate bifurcations in the geometric parameter αs at Ro = 0.7 and Ro = 0.8, to explore the relationship between the tertiary solutions at αs = 0.7 and αs = 1.5. For both rotation numbers, TWI bifurcates transcritically in αs. wTWI merges with TVF for large αs for both rotation numbers, but its behaviour at low αs is different in each case. For Ro = 0.7 in Figure 3.42, the upper and lower branch wTWI from Figure 3.36 appear, with an exchange of harmonic stability as wTWI increases in Ecf to re-join with TVF at the point where it de-stabilizes. For Ro = 0.8 in Figure 3.43 however, no upper and lower branches are found, and the solution decreases in Ecf as αs → 0, tending towards TVF with βs = 3. This suggests wTWI energes in a subcritical, subharmonic bifurcation from TVF with βs = 3. A long wavelength wTWI flowfield is shown in Figure 3.44, where it can be seen that the twists become increasingly localized as the state moves towards TVF. 90 3. Stability and bifurcations in rotating plane Couette flow Figure 3.36: Tertiary bifurcations for (αs, βs, Re) = (0.7, 1.5, 100). The global stability of the WVF solutions marked by red dots is calculated in Figure 3.40. The green dots correspond to the wTWI solutions plotted in Figure 3.38 and used in global stability calculations in Figure 3.39. Figure 3.37: Tertiary bifurcations for (αs, βs, Re) = (1.5, 1.5, 100). The purple dots correspond to the TWI solutions for global stability computations in Figure 3.41. 3.7 Bifurcations of tertiary states 91 Figure 3.38: Stable and unstable wTWI solutions at (αs, βs, Re,Ro) = (0.7, 1.5, 100, 0.7). Figure 3.39: Stability of harmonically stable and unstable wTWI solutions at (αs, βs, Re,Ro) = (0.7, 1.5, 100, 0.7). σr ≥ 0 is coloured. The harmonically stable so- lution is in fact globally stable, whilst the harmonically unstable solution is unstable to all perturbations. The least stable mode is the oscillatory harmonic mode. 92 3. Stability and bifurcations in rotating plane Couette flow Figure 3.40: Stability of harmonically stable and unstable WVF solutions at Ro = 0.05 and Ro = 0.032 for (αs, βs, Re) = (0.7, 1.5, 100). σr ≥ 0 is coloured. Both solutions lose stability to oscillatory subharmonic perturbations with the least stable mode at approximately (α, β) = (0.35, 0.135) for each state. The harmonic instability at Ro = 0.032 is also oscillatory. Figure 3.41: Stability of harmonically stable and unstable TWI solutions at Ro = 0.6 and Ro = 0.7 for (αs, βs, Re) = (1.5, 1.5, 100). σr ≥ 0 is coloured. Both solutions de-stabilize to stationary, subharmonic perturbations. The harmonic instability for the Ro = 0.7 state is also a stationary mode. 3.7 Bifurcations of tertiary states 93 Figure 3.42: Left: continuation in (αs, Ecf )-subspace for (βs, Re,Ro) = (1.5, 100, 0.7). Right: secondary stability of TVF (βs, Re,Ro) = (1.5, 100, 0.7). Figure 3.43: Left: continuation in (αs, Ecf )-subspace for (βs, Re,Ro) = (1.5, 100, 0.8). Right: secondary stability of TVF (βs, Re,Ro) = (1.5, 100, 0.8). Figure 3.44: wTWI solution with αs = 0.1 and (βs, Re,Ro) = (1.5, 100, 0.8). The twists become increasingly localized as αs → 0 and the solution approaches TVF with βs = 3. 94 3. Stability and bifurcations in rotating plane Couette flow 3.7.6 Bifurcation of wavy and oscillatory wavy vortex flows We trace the trajectories of WVF and oWVF for (αs, βs) = (0.6, 2.6) and (αs, βs) = (0.8, 3) in Figures 3.45 and 3.46. In Figure 3.45, for (αs, βs) = (0.6, 2.6) we find a trans- critical bifurcation from TVF towards WVF in the range 0.026 < Ro < 0.128 and a Hopf bifurcation towards oWVF over 0.384 < Ro < 0.7. The black curve which represents oWVF is the arithmetic mean of maximum and minium cross-flow energy that the solu- tion oscillates between during its temporal evolution, for each Ro. The oWVF solution branch is a stable attractor in this case, as it is always harmonically stable. The WVF solution branch is harmonically unstable to complex conjugate eigenmodes, prompting the existence of time-periodic quaternary states which we do not pursue here. For (αs, βs) = (0.8, 3) in Figure 3.46, the transcritical bifurcation towards WVF is similar to the previous case. The oWVF trajectory becomes more complex. Beginning with a Hopf bifurcation at Ro = 0.355, the solution reaches a turning point in (Ro,Ecf )-space at Ro = 0.438, where it loses stability and subsequently loops over parameter space, reach- ing a minimum rotation number Ro = 0.189. The state re-gains harmonic stability at Ro = 0.485 and re-joins TVF at Ro = 0.646. The periods of both oWVF states from our two bifurcation diagrams are plotted in Figure 3.47, alongside the period of the unstable linear mode which de-stabilizes TVF. Figure 3.47 suggests that in the (αs, βs) = (0.8, 3) geometry there should exist a set of states for 0.43 < Ro < 0.48 with period T ≈ 21.5, though we are unable to find any such periodic orbits to satisfactory accuracy. A situation arises at Ro = 0.45 where there do not exist any stable attractors in the do- main (αs, βs) = (0.8, 3). oWVF, TVF and the base flow are all unstable. Thus, a random perturbation settles into a quasi-chaotic motion, moving between the unstable TVF and oWVF, though not visiting the base flow. In Figure 3.48 we plot the cross-flow energy evolution for a random initial condition (a random distribution of Stokes modes) with initial energy Ecf (t = 0) = 10−6. Linear instability causes the perturbation to grow until nonlinear saturation when the fluid settles on TVF at t = 105. At t = 340 TVF breaks down resulting in a transient structure. The lower boundary of the Ecf oscillations of the transient state decreases steadily until t = 1331 when oWVF emerges. At t = 1615 oWVF breaks down and TVF re-appears, however, there is now enough background noise to cause TVF to de-stabilize on a much shorter timescale. By t = 1738 the flow has re- turned to the transient structure, and the cycle repeats again. The symmetries S and Ω are satisfied throughout the process. The dynamics are reminiscent of the regeneration cycle described by Kawahara & Kida (2001) for plane Couette flow. 3.7 Bifurcations of tertiary states 95 The stability properties of a harmonically stable and unstable oWVF solution, indicated by dots in Figures 3.45 and 3.46, are shown in Figure 3.49. We find that neither solution is globally stable. The harmonically stable solution is unstable to subharmonic modes with maximum at (α, β) = (0.45, 0.16). The frequency of the modes is tuned to the fundamental, suggesting they will oscillate with the same period as the oWVF solution. The instability of the harmonically unstable solution is strongest at the harmonic point. There is a region of spanwise subharmonic instability from 0.5 < β < 1, but it is a much weaker instability than the harmonic modes. 96 3. Stability and bifurcations in rotating plane Couette flow Figure 3.45: Tertiary bifurcations for (αs, βs, Re) = (0.6, 2.6, 100). The global stability of the solution indicated by the dot is calculated in Figure 3.49. Figure 3.46: Tertiary bifurcations for (αs, βs, Re) = (0.8, 3, 100). The global stability of the solution indicated by the dot is calculated in Figure 3.49. 3.7 Bifurcations of tertiary states 97 Figure 3.47: Period T of oWVF for (αs, βs, Re) = (0.6, 2.6, 100) and (αs, βs, Re) = (0.8, 3, 100). The solid line is the period of the orbit and the dashed line is the period of the secondary instability, 2pi|σi| . Figure 3.48: Cross-flow energy evolution of random initial condition for (αs, βs, Re,Ro) = (0.8, 3, 100, 0.45). Figure 3.49: Stability of harmonically stable (left) and unstable (right) oWVF solutions. σr ≥ 0 is coloured. Both states lose stability to subharmonic perturbations tuned to the fundamental frequency. 98 3. Stability and bifurcations in rotating plane Couette flow 3.7.7 Bifurcation of skewed vortex flow In Figure 3.50 we plot bifucations of streamwise-independent structures which bifurcate from TVF with βs = 4. We find that SVF exists for a short range of high and low Ro. At Ro = 0.085 SVF bifurcates from TVF in a transcritical bifurcation. At Ro = 0.235 SVF joins with the second Taylor vortex solution, TVF2. TVF2 is itself unstable to complex conjugate eigenmodes, prompting a Hopf bifurcation to time-periodic solutions. We do not include the unsteady solutions in Figure 3.50, as tracing the solutions became difficult due to a rapidly increasing period of oscillations. The second branch of SVF connects TVF and TVF2 at high Ro. Both the harmonically stable and unstable SVF solutions considered in Figure 3.51 lose stability to streamwise and spanwise subharmonic modes. The Ro = 0.1 solution is maximally unstable at (α, β) = (0.6, 0.2), losing stability to a complex conjugate pair of eigenmodes. An Eckhaus instability region with zero-phase modes arises close to β = 2, with local maximum (α, β) = (0, 0.2). The Ro = 0.85 solution is harmonically unstable to a complex-conjugate pair of modes. However, this instability is overtaken by zero-phase Eckhaus modes with maximum growth rate at half the spanwise wavenumber, (α, β) = (0, 2). Figure 3.50: Tertiary bifurcations for (αs, βs, Re) = (0, 4, 100) 3.7 Bifurcations of tertiary states 99 Figure 3.51: Stability of harmonically stable (top-left) and unstable (top-right) SVF solutions and TVF2 (bottom). σr ≥ 0 is coloured in (α, β). 100 3. Stability and bifurcations in rotating plane Couette flow 3.8 Transition modelling In this section we investigate the extent to which our stability analysis of secondary and tertiary states provide insight into the supercritical transitions observed in the experi- ments of Tsukahara et al. (2010) and Suryadi et al. (2013), who have mapped flow states in the (Ω ,Re)-plane. The experimenters guarded against hysteresis effects by arriving at parameter combinations in two ways; case-Re where Re is fixed while Ω is increased to its target value; and case-Ω where the rotation number is fixed and Re is slowly decreased from 1000 towards the desired value. Having settled the flow on a point in (Ω ,Re)-space the experimenters wait for some time to allow the fluid to settle before making observa- tions and characterizations of the flowfield. We analyze transition in two ways. Firstly, we perform a numerical analogue to the experiments: a long-time large periodic box simulation of a random initial flowfield. Sec- ondly, we introduce a simple model in §3.8.2 which makes a prediction of the final flow state without resorting to a large scale simulation of the Navier-Stokes equations, or in- corporating vortex splitting and merging, such as is discussed by Guo & Finlay (1991). Our first approach serves as an assessment of how well the Navier-Stokes equations with streamwise and spanwise periodic boundary conditions capture the dynamics of the ex- periments. As noted in Tsukahara et al. (2010), Ekman layers are likely to form at the lower boundary of their apparatus, causing a change in wavelengths between structures near the lower wall and those away from it. An Ekman layer is formed at a solid wall in a rotating fluid, and their effects have been discussed in Taylor-Couette experiments (Czarny et al. (2003), Hollerbach & Fournier (2004), Altmeyer et al. (2010), Heise et al. (2013)), however, we do not include the effects of Ekman layers in this analysis. The focus of this study is the transition between laminar flow states, prior to transitions to quasi or fully turbulent dynamics. As such, we fix Re = 100 and investigate changes in the flow structures with respect to the rotation parameter. Tsukahara et al. (2010) found four flow states in the region 0 < Ω < 30, which corresponds to 0 < Ro < 0.3 (see Figure 3.52 for a re-print of their state-space map). In increasing Ω , the flow states are labelled COU3Ds (spatio-temporally intermittent three-dimensional roll-cells) from 2 < Ω < 7, COU3D (three-dimensional roll-cells) from 7 < Ω < 12, COU2Dm (two- dimensional meandering roll-cells) from 12 < Ω < 15 and COU2Dh (two-dimensional roll-cells at high rotation number) from 15 < Ω < 30. Photographs of the flow states COU3Ds, COU3D and COU2Dh are re-printed from Tsukahara et al. (2010) in Figures, 3.53, 3.54 and 3.55. Suryadi et al. (2013) have extended the rotation parameters and found that for Re = 100: two-dimensional roll-cells persist over 30 < Ω < 70, and for 3.8 Transition modelling 101 Ω = 90 two distinct patterns are observed at different times. Photographs of the two flow states as observed in the latter experiments, which we will label (A) and (B), are printed in Figure 3.56. The photographs were kindly sent by private communication with the authors. Of the two flowfields, (A) has a much greater degree of streamwise structure in comparison with (B), which contains localized folds between streaky structures. Figure 3.52: Re-print of Figure 3 from Tsukahara et al. (2010), showing the flow states observed at different locations in the (Re,Ω) parameter space. We focus on the transition between states with fixed Reynolds number, Re = 100. 102 3. Stability and bifurcations in rotating plane Couette flow Figure 3.53: Re-print of Figure 10 from Tsukahara et al. (2010), showing the flow state COU3Ds for (Re,Ω) = (151,1.32). (a), (b), (c) and (d) correspond to images recorded at time intervals of 0, 12, 24 and 50 seconds. Figure 3.54: Re-print of Figure 9(a) from Tsukahara et al. (2010), showing the flow state COU3D for (Re,Ω) = (99.6,8.65). Figure 3.55: Re-print of Figure 14 from Tsukahara et al. (2010), showing the flow state COU2Dh for (Re,Ω) = (178,24.6). 3.8 Transition modelling 103 Figure 3.56: Photographs of the two distinct flow states observed experimentally by Suryadi et al. (2013) for (Re,Ω) = (100, 90). (Private communication). 3.8.1 Numerical procedure For each numerical simulation we have used a domain D, with fundamental wavenumbers αf = 0.2 and βf = 0.4 D = [0, 10pi]× [−1, 1]× [0, 5pi], (3.78) which is discretized as Nx × Ny × Nz = 256 × 49 × 128. This discretization is chosen such that the number of points per unit length is in keeping with other large domain simulations, for example Duguet et al. (2010) who used Nx×Ny×Nz = 2048×33×1024 in the very large domain [0, 800] × [−1, 1] × [0, 356]. The precise geometry of the ex- perimental apparatus is not reported for each case in Tsukahara et al. (2010), since the wall-normal gap is modified for each Re. However, the streamwise and spanwise lengths are reported as 1500mm and 360mm respectively, so a wall-normal width of 10mm, for example, would imply a domain of 300h × 2h × 72h. Thus D is much smaller than the 104 3. Stability and bifurcations in rotating plane Couette flow experimental domain, which can of course introduce a reason of discrepancy between our simulations and the experiments. A random initial perturbation of the base flow, u0, is generated using u0 = Nx/2∑ j=−Nx/2+1 Nz/2∑ k=−Nz/2+1 Ny−1∑ `=0 µj,k,sT`(y)e i(jαfx+kβf ), (3.79) where T`(y) is the `th Chebyshev polynomial of the first kind, and µj,k,s = λ(1− s) |j|+|k|+|`|, (3.80) where λ is a random number with λ ∈ [−1, 1]. s is a parameter controlling the smoothness of the field. We set s = 0.2 and initial energy E0 = 5× 10−2 for the initial perturbation. The same field is used as the initial condition for each simulation. 3.8.2 Model prediction We now introduce a simplified model of the transition scenario which we use to make a predicition of the final flow state at each Ω . We begin by calculating the spectral energy E0αβ of the initial disturbance, u0(x, y, z). For vαβ(y) = ∫ 2pi α 0 ∫ 2pi β 0 u0(x, y, z)ei(αx+βz) dzdx, (3.81) we have E0αβ = 1 2 ∫ 1 −1 vαβv ∗ αβ dy. (3.82) From linear theory, the growth rate ωαβ at which disturbances are linearly amplified about the linear base profiled is calculated. During the primary linear growth phase we expect the energy of each mode to grow as E(t) = E0αβe 2ωαβt. (3.83) We then define a saturation energy, Es, being the energy at which we expect nonlinear effects to dominate the dynamics. Since disturbances with α = 0 have the largest linear growth rates (cf §3.2), we assume that a spanwise wavenumber, βII, will emerge from the linear growth phase to dominate the secondary dynamics. To determine βII, we solve for the saturation time ts ts(β) = 1 2ωβ log ( Es E00β ) . (3.84) 3.8 Transition modelling 105 Then βII is defined to be the wavenumber which minimizes the saturation time ts(β II) = min β ts(β). (3.85) We then assume a secondary Taylor vortex solution with spanwise wavenumber βII to emerge, which we will denote, U II(α, β) = U II(0, βII). If the vortex is unstable to streamwise disturbances, a streamwise wavenumber αIII will appear which maximizes the growth rate σ of a linear perturbation to the Taylor vortex σ(αIII) = max α σ(α). (3.86) A third solution, U III(αIII, βII), can then be computed. From a further Floquet analysis on this third state we determine wavenumbers αIV and βIV which satisfy σ(αIV, βIV) = max α,β σ(α, β), (3.87) αIV ≤ αIII, βIV ≤ βII, (3.88) where σ is again the growth rate of a linear perturbation to the third state. A new state, U IV(αIV, βIV) is found by repeating U III over the wavenumber domain (αIV, βIV), adding a random perturbation to the flowfield, and then evolving towards a steady state or periodic orbit. The random perturbation selects the natural bifurcation away from U III towards U IV, without the need to perturb the system with specified eigenmodes from the Floquet analysis. This process is continued until a globally stable solution is found, or until finding a solution becomes too computationally expensive. We expect this solution to feature flowfields of the simulations and experiments. We have refrained from defining the progression of solutions as tertiary and quaternary in favour of III, IV, . . . , since we have defined tertiary solutions in §3.7 to be structures which originate in a bifurcation from secondary states, meaning IV and higher states could be classified as tertiary solutions in keeping with this convention. In the analysis which follows we choose saturation energy Es = 0.00125, on the assump- tion that velocities with amplitude As = √ 2Es = 0.05 give a reasonable estimate for the onset of nonlinear interactions. For the stability calculations we need only compute growth rates for wavenumbers that are harmonics of the fundamentals, since they are the only ones supported by the domain α = 0, 0.2, 0.4, 0.6, . . . (3.89) β = 0, 0.4, 0.8, 1.2, . . . . (3.90) 106 3. Stability and bifurcations in rotating plane Couette flow 3.8.3 Time-periodic flow: Ω = 2.5 For Ω = 2.5 our simulation gives an unsteady, quasi-periodic flowfield, with Ecf evolution reaching a local maximum approximately every 200−300 time units, shown in Figure 3.57 along with snapshots of the flowfield at certain times. The fluid is arranged into vortices with spanwise wavenumber β = 1.2 which oscillate in time, switching between wavy vortices and streamwise independent vortices during the oscillation. The streamwise variation of the vortices does not appear to be periodic on a smaller scale than the domain length, suggesting that the structure has streamwise wavenumber α = 0.2. The dynamics could be crudely interpreted as a small-domain version of the COU3Ds observed by Tsukahara et al. (2010) for 2 < Ω < 7, who describe an intermittent process of three-dimensional (wavy) vortices changing to two-dimensional (streamwise independent) vortices and so on. In the transition model we find that the wavenumber minimizing the saturation time is βII = 1.2, so we set the secondary solution, U II(0, 1.2), to be TVF. U II is unstable to streamwise periodic perturbations, and its least stable mode is stationary, with wavenumber αIII = 0.4. Thus U III(0.4, 1.2) is a steady WVF, and we find that U III(0.4, 1.2) loses stability in a harmonic Hopf bifurcation, which generates an oscillatory wavy vortex flow (oWVF)U IV(0.4, 1.2) with period Tp = 148.71. U IV(0.4, 1.2) has a subharmonic instability to complex conjugate modes with half the natural frequency of U IV, implying a structure will emerge which is periodic in time with the wavenumber combination (αV, βII) = (0.2, 1.2)., though the period of the solution is too high to analyzed here. However, the results of the transition model would appear to be in good agreement with the simulation which, as described above, contained oscillating vortices with wavenumber (α, β) = (0.2, 1.2). The model results are contained in Figure 3.58, where we plot the saturation time against spanwise wavenumber and the growth rates from the stability calculations about the solutions U II, U III and U IV, respectively. 3.8 Transition modelling 107 Figure 3.57: Simulation results for Ω = 2.5. The Ecf evolution, plotted at the top, clearly displays the unsteadiness of the flowfield, sustained over 2000 time units. From the snapshots of the flowfield in the four bottom panels, we can see that the wavy pattern of the vortices oscillates in time. cf with the oscillatory parts of COU3Ds in Figure 3.53. 108 3. Stability and bifurcations in rotating plane Couette flow Figure 3.58: Results of the transition model for Ω = 2.5. (a) Saturation time ts against spanwise wavenumber β. βII = 1.2 minimizes ts. In (b), (c) and (d) we plot the growth rates σr of the least stable perturbations of the solutions U II (TVF), U III (WVF) and U IV (oWVF) for all the subharmonic wavenumbers supported in D. The least stable mode throughout D is circled in each panel. 3.8 Transition modelling 109 3.8.4 Wavy vortex flow: Ω = 10 For Ω = 10, the simulation finishes in a stable WVF state with wavenumbers (α, β) = (0.4, 1.2). The cross-flow energy evolution and the final flowfield are shown in Figure 3.59, with the assembly of wavy vortex structures clearly visible in the velocity field. The stable WVF arrangement is in strong agreement with Tsukahara et al. (2010), who found stable wavy vortex flow, labelled COU3D in the range 7 < Ω < 12. The transition model additionally predicts a WVF final state. In the model, βII = 1.2 is selected by minimization of ts, and the TVF U II(0, 1.2) de-stabilizes to a streamwise periodic wave with αIII = 0.4. U II(0, 1.2) subsequently transitions to the WVF U III(0.4, 1.2) which is stable in D, as shown in Figure 3.60. We note that the streamwise wavenumber αIII = 0.4 produces a wavelength xL = 2piαIII ≈ 15.71, which is in good agreement with the wave- length of the wavy vortices of COU3D in Figure 3.54. Figure 3.59: Simulation results for Ω = 10. At the top is the Ecf evolution. On the bottom is the final flowfield in D, which consists of a stable configuration of wavy vortices. cf with the COU3D flowfield in Figure 3.54. 110 3. Stability and bifurcations in rotating plane Couette flow Figure 3.60: Transition model results for Ω = 10. In (a), (b), (c) and (d) we plot the saturation time ts against spanwise wavenumber β; the stability of U II(0, 1.2) in D; the stability of U III(0.4, 1.2); and the stable wavy vortex flowfield corresponding to U III(0.4, 1.2). The red circles in (a) and (b) indicate the wavenumbers selected in each step. 3.8 Transition modelling 111 3.8.5 Taylor vortex flow: Ω = 20 As the system rotation is increased to Ω = 20, the simulation in D finishes with a stable TVF with βs = 1.2. We interpret this as an example of the COU2Dh Tsukahara et al. (2010) observed for 15 < Ω < 30. The cross-flow energy evolution and the flowfield consisting of an array of Taylor vortices are plotted in Figure 3.61. The transition model once again selects the spanwise wavenumber βII, and we find that there are no secondary instabilities about the Taylor vortex at this spanwise wavelength and rotation parameter. The model results are shown in Figure 3.62. Thus the experiment, simulation and model all produce stable Taylor vortices in this case. Figure 3.61: Simulation results for Ω = 20. At the top is the Ecf evolution. On the bottom is the final flowfield: a stable arrangement of Taylor vortices. cf with the COU2Dh flowfield in Figure 3.55. 112 3. Stability and bifurcations in rotating plane Couette flow Figure 3.62: Transition model results for Ω = 20. (a) Saturation time ts against spanwise wavenumber β. (b) Secondary stability of the TVF U II(0, 1.2). 3.8.6 Taylor vortex flow: Ω = 50 From the energy evolution of the simulation in Figure 3.63 we see that for Ω = 50, the fluid has a lengthy unsteady phase over the first 3000 time units, after which there is an abrupt transition to a stable steady state. The unsteady flowfield contains a mixture of streaks and twist vortices, with a localized street of twist vortices settling in the region 0 < x < 10, 5 < z < 12 by t = 2000, and the steady flowfield is an arrangement of Taylor vortices with spanwise wavenumber β = 1.6, as is shown in Figure 3.63. The final flow structure agrees with the experimental observations of Suryadi et al. (2013), who report stable two-dimensional streamwise vortices persisting in the range 40 < Ω < 70, though structures similar to the unsteady flowfields, t = 1000, 2000, 3000 in Figure 3.63, have been observed at greater Reynolds numbers, Re = 150, 200. However, the transition model does not predict stable Taylor vortices in this case, instead predicting stable twist vortices. The model therefore offers an explanation as to why localized twist vortices appear in the simulation, though it is beyond the scope of the transition model to accurately describe the localized behaviour in D for t < 3000. Once again minimization of the saturation time gives the spanwise wavenumber βII = 1.2, and the Taylor vortices with this wavenumber lose stability to both wavy twist and twist modes. The wavy twist mode with streamwise wavenumber αIII = 1.6 is most unstable, generating the wTWI solution U III(1.6, 1.2), which in turn loses stability to a harmonic twist mode which initiates a bifurcation to the TWI solution U IV(1.6, 1.2). U IV(1.6, 1.2) is stable in D, as is shown in Figure 3.64, and the flowfields of both U III(1.6, 1.2) and U IV(1.6, 1.2) are displayed in Figure 3.65. 3.8 Transition modelling 113 Figure 3.63: Simulation results for Ω = 50. At the top is the Ecf evolution. On the bottom is the final flowfield: a stable arrangement of Taylor vortices. 114 3. Stability and bifurcations in rotating plane Couette flow Figure 3.64: Transition model results for Ω = 50. In (a), (b), (c) and (d) we plot the saturation time ts against spanwise wavenumber β; the stability of U II(0, 1.2) in D; the stability of U III(1.6, 1.2); and the stability of U IV(1.6, 1.2). The red circles in (a), (b) and (c) indicate the wavenumbers selected in each step. Figure 3.65: Flowfields of U III(1.6, 1.2) and U IV(1.6, 1.2) from the transition model for Ω = 50. 3.8 Transition modelling 115 3.8.7 Wavy twist vortex flow: Ω = 90 The large simulation and the transition model do not agree in the case Ω = 90, how- ever, the transition model predicts a streamwise periodic wavy twist vortex final state and the simulation finishes on a flowfield with localized folds, which are similar to the experimentally observed structures (A) and (B), respectively, in Figure 3.56. The results of the simulation are displayed in Figure 3.66. After saturation, the flowfield at t = 500 contains an upper and lower pair of folds along its centre, with streamwise indepenent streaks of wavenumber β = 1.6 completing the domain on either side of the folds. The folds resemble the long wavelength wavy twist vortex solution found in §3.7 in Figure 3.44. As the simulation evolves, the upper fold collapses, and a short wavelength wavy twist vortex can be seen at t = 1500 before the complete disappearance of the upper fold, leaving only the lower fold surrounded by Taylor vortices. From Floquet theory we fine that the final structure is stable. The transition model predicts the following transition of states: βII = 1.2 initiates TVF U II(0, 1.2) which loses stability to a wavy twist mode with αIII = 0.8; wTWIU III(0.8, 1.2) loses stability to subharmonic streamwise wavenumber αIV = 0.4 which bifurcates to- wards a second wTWI U IV(0.4, 1.2); a further subharmonic streamwise wavenumber αV = 0.2 de-stabilizesU IV(0.4, 1.2) prompting a bifurcation to a third wTWIUV(0.2, 1.2) which is stable in D. The stability results are given in Figure 3.67 and the flowfields of the three wTWI solutions are plotted in Figure 3.68. In Figure 3.69, the final flowfields from the model prediction and the simulation are plotted on a larger domain and placed next to their experimentally observed fields (A) and (B). We find that the model prediction provides a reasonable match for the behaviour seen in flowfield (A), while the simulation captures the folds observed in flowfield (B). 116 3. Stability and bifurcations in rotating plane Couette flow Figure 3.66: Simulation results for Ω = 90. At the top is the Ecf evolution. Beneath we plot flowfields for t = 500, 1000, 1500 and 2000, showing the decay of the upper fold and persistence of the lower fold. 3.8 Transition modelling 117 Figure 3.67: Transition model results for Ω = 90. In (a), (b), (c), (d) and (e) we plot the saturation time ts against spanwise wavenumber β; the stability of U II(0, 1.2) in D; the stability of U III(0.8, 1.2); the stability of U IV(0.4, 1.2); and the stability of UV(0.2, 1.2). The red circles in (a), (b) and (c) indicate the wavenumbers selected in each step. Figure 3.68: Flowfields of the wTWI states U III(0.8, 1.2), U IV(0.4, 1.2) and UV(0.2, 1.2). 118 3. Stability and bifurcations in rotating plane Couette flow Figure 3.69: Comparison between our model/simulation flowfields and experimentally observed flow structures. The model and simulation flowfields are doubled on to a grid which is roughly comparable to the experimental photograph size as determined by the 2h gap given in flowfield (B). The patterns in the model and simulation final flowfields are in reasonable agreement with flowfields (A) and (B), respectively. 3.8.8 Assessment Our attempts to understand the transition of flow states have had mixed success. The periodic-box simulations have captured the unsteady wavy vortex, steady wavy vortex and Taylor vortex flow regimes at parameter values in agreement with Tsukahara et al. (2010) and Suryadi et al. (2013), and there is some indication that the simulations re- produce the fold features observed by Suryadi et al. (2013) for Ω = 90. However, the unsteady wavy vortices (COU3Ds) are described by Tsukahara et al. (2010) as being spatio-temporally intermittent, which was not a feature produced by our simulations. In addition, we did not find any behaviour resembling COU2Dm, described in Tsukahara et al. (2010) as a “meandering” two-dimensional roll-cell structure, in our simulations. Discrepancies between simulation and experiment could have come about in three ways. Firstly, the periodic boundary conditions imposed in the simulation may not accurately represent the experimental domain, which would have solid walls at its boundaries. As mentioned previously, Ekman layers are known to form at rotating solid boundaries, and the effect of Ekman layers in the rotating plane Couette experiments has not, to our 3.8 Transition modelling 119 knowledge, been investigated. Secondly, the domain size in our simulations was much smaller than that of the experimental test section. A larger domain could introduce in- stabilities subharmonic to D, which affect the dynamics. Furthermore, a larger domain is likely necessary to observe spatio-temporal disruption such as COU3Ds. Thirdly, we chose a randomly generated flowfield as our initial disturbance. Clearly, in a deterministic dynamical system which supports chaotic dynamics, such as the Navier-Stokes equations, the final flowfield may be strongly dependent on the initial state. Even without chaos, the evolution of a system could be highly dependent on the initial energy spectrum, in such a way that wavenumbers with a high initial amplitude may dominate the dynamics. Our initial condition was chosen with a degree of smoothness, which gives the velocity field a physically acceptable appearance but prevents larger wavenumbers from having a high energy excitation. An alternative approach may be to first simulate a turbulent field at a higher Reynolds number and subsequently use this as an inital condition at the lower Reynolds number, since this is exactly what is done in experiments. Simula- tion a fully turbulent flowfield, however is an expensive task. Hence, a study which can incorporate non-periodic boundary conditions, a larger domain size and a greater range of initial conditions might better replicate some of the flowfields observed experimentally. The transition model introduced in §3.8.2 has a varied performance when compared both to simulations and experiment. In the cases Ω = 2.5, 10, 20 there is encouraging agree- ment between all three of experiment, simulation and model. For Ω = 50, 90 the model fails to predict the final flowfield of the simulations, though a structure similar to the predicted state is observed transiently for Ω = 50, and a structure observed in the ex- periments of Suryadi et al. (2013) resembles the model prediction for Ω = 90. Given the simplicity of the model, this could be construed as reasonable success. Furthermore, given the purpose of the model is to provide insight into the complicated dynamics that a fluid experiences on evolution under the Navier-Stokes equations, which solving the equa- tions in a direct numerical simulation does not always offer, we believe it to be adequate. Multiple stable structures can exist as solutions to the equations at given parameters, thus an improvement to the model procedure could be made by choosing not one but a range of wavenumbers which are selected at each stage of transition, and forming a probabilistic model of the emergent states. We note that the states oWVF and SVF from §3.7 do not appear in any of the simulations, model bifurcation sequences or experimental results. Our model results indicate that, for the chosen initial flowfield, the dominant secondary wavenumber is always βII = 1.2. However, from §§3.6 and 3.7 it is clear that βII > 1.2 would be required for oWVF or SVF to appear. Hence, to assess the effect of these solutions on large-domain transition, 120 3. Stability and bifurcations in rotating plane Couette flow an initial condition which selected a larger secondary wavenumber in a simulation would be necessary. 3.9 Conclusion In this chapter we have characterized a range of the secondary and tertiary nonlinear states that emerge in RPCF. In §3.5 we introduced three secondary states, TVF, TVF2 and OVF. We produced bifurcation diagrams of these states to investigate their behaviour as the rotation parameter Ro is varied, finding that while TVF and TVF2 undergo rela- tively straightforward bifurcations in Ro, OVF bifurcates with an upper and lower branch solution which merge as the streamwise wavenumber αs is increased and leave behind a disconnected low-Ro and high-Ro set of solutions. The least stable perturbation of the primary, linear Couette profile in RPCF initiates a bifurcation to TVF, meaning TVF is the most likely secondary structure to be observed. Therefore, in §3.6 we performed an extensive Floquet stability analysis on TVF, at three Reynolds numbers identified to be in the transitional flow regime by the experimental study of Tsukahara et al. (2010). The stability properties of TVF are shown to be depen- dent on Re, Ro and the wavenumber of the vortex, and instabilities increase in range and magnitude as Re is increased, giving indication of a wide range of tertiary instabilities possible over the experimentally observed transition region. Tertiary states which bifurcate from TVF are studied in terms of their bifurcation and stability behaviour in §3.7. With focus on bifurcation in Ro, we see WVF bifurcate at low Ro, oWVF bifurcate for mid-Ro and TWI and wTWI bifurcate at high Ro. We find a streamwise independent structure SVF, which joins TVF and TVF2 in its bifurcation in rotation number. wTWI is found to emerge in a saddle-node bifurcation at low α, and after a continuation of the solution in αs the solution branch moves towards the TVF with double the spanwise wavenumber β, suggesting a connection between TVF and wTWI that has not previously been discussed. A regeneration cycle is found to occur for (αs, βs, Re,Ro) = (0.8, 3, 100, 0.45). After TVF with β = 3 breaks down, no stable tertiary state exists and the fluid moves between the unstable TVF, an unstable oWVF and a transient state with a slowly varying amplitude. A similar process has been found to occur in Kawahara & Kida (2001) in non-rotating Couette flow at Re = 400. Our finding could mean that the process “originates” in RPCF and can be continued over to the non-rotating flow in a manner similar to Nagata (1990, 2013). Finally in §3.8 we explored the extent to which local stability analysis of the secondary 3.9 Conclusion 121 and tertiary structures from §3.6 and 3.7 can be used to describe the flow regimes observed in a large periodic-box simulation, and the experiments of Tsukahara et al. (2010) and Suryadi et al. (2013). We performed direct numerical simulation of the Navier-Stokes equations at specific parameter values, finding good agreement between the simulations and experiments. A simple transition model was proposed to predict the form of the final flowfield by assuming a transition between least stable states, with reasonable consistency between experiment and simulation. 122 3. Stability and bifurcations in rotating plane Couette flow Chapter 4 Non-modal growth in rotating plane Couette flow In this chapter we continue our analysis of stability in rotating plane Couette flow (RPCF) by extending existing non-modal growth analyses across the the entire subcritical param- eter space. Our aim is to identify the maximally amplified, optimal linear disturbances in subcritical (Re,Ro)-space. There are three known linear energy amplification mecha- nisms present in RPCF; the lift-up mechanism, the anti lift-up mechanism and the Orr mechanism. In determining the maximally amplified linear disturbances, we can indicate where, in parameter space, each growth mechanism dominates the transient dynamics. We find that for strong rotation, whether cyclonic or anti-cyclonic, the linear optimals two-dimensionalize and grow via the Orr mechanism. The chapter outline is as follows: in §4.1 we introduce the mathematical background for non-modal growth and the compu- tational procedures we used to determine it, in §4.2 we present our findings of the optimal disturbances under both cyclonic and anti-cyclonic rotation and in §4.3 we discuss the relevance of our analysis to the problem of subcritical hydrodynamic transition. 4.1 Non-modal growth The equations governing RPCF are stated in §3.1, and we do not repeat them here. Consider the linear stability equations (3.8). Fixing α, β,Re,Ro ∈ R we can solve for ω` ∈ C satisfying Aqˆ` = ω`Bqˆ`, (4.1) with the operators A and B defined as in §3.2 and the vector qˆ`(y) defined as qˆ` = (uˆ`, p`). We can then construct an arbitrary velocity perturbation as an expansion of the first L, say, modes u(x, y, z, t) = ei(αx+βz) L∑ `=1 a`uˆ`(y)eω`t. (4.2) 123 124 4. Non-modal growth in rotating plane Couette flow The physical perturbation kinetic energy E˜(t) through the streamwise and spanwise periodic domain [0, 2piα ]× [−1, 1]× [0, 2pi β ], of a velocity u, is given by E˜(t) = 2ρpi2 αβ ∫ 1 −1 uHu dy, (4.3) where ρ is the constant density of the fluid and H denotes the Hermitian transpose. We seek to maximize the gain in kinetic energy at arbitrary time t over all inital velocity profiles. To this end we introduce the kinetic energy-density E(t) = αβE˜ 2ρpi2 = ∫ 1 −1 uHu dy = ‖u(t)‖2. (4.4) For aˆ = (a1, a2, . . . , aL) and Λ = diag(ω1, ω2, . . . , ωL), E(t) can be written E(t) = aˆHeΛ H tMeΛtaˆ, (4.5) where eΛt is the matrix exponential of Λt, and M is the L× L matrix Mjk = ∫ 1 −1 uˆHk uˆj dy. (4.6) M is both Hermitian (M = MH) and positive definite. Since M is Hermitian, its singular value decomposition can be written M = UMΣMU H M , (4.7) where ΣM = diag(σ1, σ2, . . . , σL) is a diagonal matrix of the singular values of M in descending order σ1 > σ2 > . . . > σL, and UM is a unitary matrix. For √ ΣM = diag( √ σ1, √ σ2, . . . , √ σL) we have F = UM √ ΣM , (4.8) which allows us to factor M into M = FHF . Hence we have ‖u(t)‖2 = aˆHeΛ H tFHF eΛtaˆ = ‖F eΛtaˆ‖22, (4.9) 4.1 Non-modal growth 125 where ‖.‖2 is the Euclidean 2-norm. The maximum energy amplification G is simply G(t) = max u0 6=0 E(t) E(0) = max u0 6=0 ‖u(t)‖2 ‖u0‖2 , (4.10) = max aˆ6=0 ‖F eΛtaˆ‖22 ‖F aˆ‖22 , (4.11) = max aˆ6=0 ‖F eΛtF−1F aˆ‖22 ‖F aˆ‖22 , (4.12) = ‖F eΛtF−1‖22. (4.13) G(0) = 1 by definition. Linear non-modal energy growth is supported in a flow if G(t) > 1 for any t. G represents the maximum possible energy amplification at time t of an initial perturbation, though different initial perturbations may be responsible for maximum growth at distinct times. Suppose, having found the maximum energy amplification achieveable at time tq, we wish to find the form of the optimal initial perturbation and also its form at tq. Let q0 denote the optimal perturbation, with q0 = e i(αx+βz) L∑ `=1 a`qˆ`, and ‖q0‖ = 1. (4.14) Then since q0 will have energy ‖F e ΛtqF−1‖22 at tq, we have Baˆ = ‖F eΛtqF−1‖2aˆe Λtq = µaˆeΛtq , (4.15) where B = F eΛtqF−1. The constant µ = ‖F eΛtqF−1‖2 is simply the principal singular value of B. The singular value decomposition of B lets us write BVB = UBΣB, (4.16) for unitary matrices VB and UB, and the matrix of singular values of B, ΣB = diag(σ1, σ2, . . . , σL). If we denote the principal eigenvectors of VB and UB as v1 and u1 respectively, then since µ = σ1, we have Bv1 = σ1u1, (4.17) from which we can identify aˆ = v1. The linear optimal q0 can then be found by expanding aˆ over the basis vectors qˆ1 . . . qˆL. Our optimal growth formulation follows that of Schmid & Henningson (2001) which uti- lizes the singular value decomposition. An alternative method exists in the variational Lagrangian formulation, and such an approach allows for a nonlinear or time-dependent 126 4. Non-modal growth in rotating plane Couette flow base flow about which optimal growth can be calculated. The variational method has been used to calculate optimal disturbances in complex flows by, for example, Corbett & Bottaro (2001) and Gue´gan et al. (2006) and to calculate nonlinear optimals by Pringle et al. (2012), Rabin et al. (2012) and Duguet et al. (2013). However, since the base flow under consideration in this work is parallel and time-independent, we felt the singular value decomposition approach was appropriate for its ease of implementation. To compute G, we discretize the linear system (3.8) in the wall-normal direction using an expansion of N Chebyshev polynomials. The resultant finite eigenproblem is solved for eigenvalues ω` and eigenvectors qˆ`, for ` ∈ [0, L]. The eigenvectors are used to compute the matrix M and from the singular value decomposition of M we compute F . Thus we calculate G using the matrix norm in equation (4.13). In Figure 4.1 we plot G as a function of mode truncation L and wall-normal discretization N , with t, α, β, Re and Ro held fixed. The curves for N = 49 and N = 65 are almost indistinguishable, and there is little change in G for L > 60. This is an indication that G is well-converged at these parameters. Accordingly, in our calculations we use discretization/truncation (N,L) = (49, 60). Figure 4.1: Convergence of G at various wall-normal discretizations N . (t, α, β, Re,Ro) = (21.6, 0.1644, 1.6914, 1000, 1.01). 4.1.1 Maximum optimal growth Note that G is a five parameter function, G = G(t, α, β, Re,Ro). In recognition of this, G can be maximized over wavenumbers α, β and time t. We define the temporal maxima 4.1 Non-modal growth 127 Gmax and tmax Gmax(α, β,Re,Ro) = max t G(t, α, β, Re,Ro) = G(tmax, α, β, Re,Ro) (4.18) and the wavenumber maxima G, αmax and βmax G(Re,Ro) = max α,β Gmax(α, β,Re,Ro) = G(Tmax, αmax, βmax, Re,Ro). (4.19) G indicates the greatest linear amplification possible at a given Re, Ro combination and we shall refer to the initial condition which reaches this energy amplification as the linear optimal. Tmax is used in (4.19) to denote the time taken to achieve G, and is to be differentiated from tmax, which is the time taken to reach the maximum G for an arbitrary wavenumber combination. Optimization routines in MATLAB are used to find both Gmax and G. Finding Gmax requires optimization over one variable, so the golden mean search fminbnd is used over the time horizon 0 ≤ t ≤ Re. In some cases G can have local maxima. By varying the time horizon of the optimization from one iteration to the next, we increase the chances of finding global maxima. Maximizing G requires optimization over two variables; for this the interior-point algorithm within the function fmincon is used, subject to the constraints α, β ≥ 0. In order to lessen the likelihood of finding a local maximum of Gmax, rather than the global maximum, five random wavenumber pairs are chosen in the domain 0 ≤ α ≤ 2, 0 ≤ β ≤ 2 as initial conditions for the optimization. The algorithm is then applied to each inital condition, with G presumed to be the maximum of the outcomes of the five iterations. In cases where αmax < 0.1, we perform golden mean searches over β with α = 0 to check for a streamwise independent maximum. A golden mean search over α is similarly performed when βmax < 0.1. In all optimizations, iterations k are continued until k = K satisfying |fK − fK−1| < 10−10, for f = Gmax or f = G. 128 4. Non-modal growth in rotating plane Couette flow 4.2 Optimal growth in RPCF Figure 4.2: Contours of G in (Ro,Re). The black curves are contours of G = 100.125, 100.25, 100.375, . . . , 104.5 The solid red curve represents the energy stability bound- ary and the dashed red curve is the neutral stability boundary The central result of this chapter is presented in Figure 4.2, where we plot level curves of G in the (Ro,Re) parameter space. Contours of G, ranging from 100.125 to 104.5 in steps of 100.125, are computed from a grid of 72 × 23 points spanning the subcritical (Ro,Re)-space. The neutral curve, which denotes the boundary between linearly stable and unstable modes is plotted as the red dashed line. An unstable perturbation has unbounded maximum growth, G = ∞, since unstable perturbations grow exponentially in time. For this reason we do not include contours in the unstable region. The thick red line in Figure 4.2 corresponds to the energy stability curve, which we define for fixed Ro as the maximum Re for which G = 1. Below the energy stability curve non-modal growth is not supported, thus the energy of any linear disturbance is strictly decaying. The curve is traced via a bisection method, where we fix Ro and calculate the corre- sponding Re to a tolerance of  = 10−6. The energy stability Reynolds number ReE, defined here as the minimum Reynolds number on the energy stability curve, is found to be ReE ≈ 20.67. By comparison, the critical Reynolds number Recr, i.e. the minimum Re for which the flow is neutrally stable, is Recr ≈ 20.69. Both numbers are plotted in Figure 4.3, which is a blow-up of the 0 < Ro < 1 region at low Re. Both Recr and ReE 4.2 Optimal growth in RPCF 129 are found at Ro = 0.5, perhaps because both the neutrally stable and the energy stable perturbations along this section of the curve are streamwise oriented vortices, with α = 0 and β ≈ 1.55. From Figure 4.2 it can be seen that for strong rotation, with |Ro| far from the neu- tral curve, linear energy amplification is severely weakened and the value of G becomes independent of rotation. This serves as the main result of this chapter. Figure 4.3: The energy stability curve (solid) and the neutral curve (dashed) are plotted with the energy stability Reynolds number ReE and the critical Reynolds number Recr indicated. 130 4. Non-modal growth in rotating plane Couette flow 4.2.1 Cyclonic rotation and the non-rotating case Figure 4.4: Optimal growth in the cyclonic rotation regime. Contours are equivalent to those in Figure 4.2. Under cyclonic rotation G decreases sharply for Re held fixed and Ro decreasing from zero. In the absence of rotation we find that G ∝ Re2, and we plot this relationship in Figure 4.5 along with the wavenumbers αmax and βmax. This scaling is indictative of energy amplification via the lift-up mechanism (see, for example, Landahl (1980), Gustavsson (1991), Schmid & Henningson (2001)), which refers to extraction of energy from the mean shear profile as low-energy streamwise oriented vortices are ‘lifted-up’ into high-energy streamwise streaks. The mechanism is necessarily three-dimensional. As a representative case, the flowfield of the optimal with G = 959.5 at (Ro,Re) = (0, 900) is depicted in Figure 4.6. Note that αmax = 0.0381, and from Figure 4.5 that for Re ≥ 100, αmax 6= 0, which indictates that the optimal is not streamwise independent but weakly oblique. This trend is in agreement with the observations of Butler & Farrell (1992), who noted that the disturbance achieving the most growth in plane Couette flow used both the lift-up mechanism in addition to the Orr mechanism, Orr (1907). The streamwise wavenumber follows αmax ∝ Re−1, a scaling derived by Chapman (2002). The Orr mechanism is two-dimensional: streamwise periodic rolls (u, v) lean against the mean shear and are amplified as they ‘swing’ to align with the shear. However, the combined optimal decays faster than the pure streamwise vortex. Using the streamwise, wall-normal 4.2 Optimal growth in RPCF 131 and spanwise energy fractions, defined as Eu = ∫ 1 −1 u ∗u dy ∫ 1 −1 u Hu dy , Ev = ∫ 1 −1 v ∗v dy ∫ 1 −1 u Hu dy , Ew = ∫ 1 −1w ∗w dy ∫ 1 −1 u Hu dy , (4.20) with ∗ denoting the complex conjugate, the optimal at t = 0 has Ev = 0.3282 and Ew = 0.6807 respectively, with a low amount of energy in the streamwise velocity Eu = 0.0011. At t = Tmax however, the energy budget has changed to Eu = 0.9991, Ev = 0.0001 and Ew = 0.0008, which highlights the transition from a streamwise vortex to a high-energy streak. Figure 4.5: Left: G in Re at Ro = 0. The dashed line depicts the scaling G ∝ Re2. Right: wavenumbers αmax (blue) and βmax (red) for Ro = 0. Figure 4.6: (z, y)-profiles of the optimal disturbance at (Ro,Re) = (0, 900) with (αmax, βmax) = (0.0381, 1.6). The velocity at t = 0 (top) and t = Tmax = 105.7 (bottom) is shown. Each velocity is represented as a colour-plot of the streamwise velocity u, with red and blue indicating positive and negative, and a vector plot of the wall-normal and spanwise velocities (v, w). 132 4. Non-modal growth in rotating plane Couette flow For Ro ≈ −20 we find G ∝ Re2/3, in agreement with the scaling found by Yecko (2004) for two-dimensional disturbances under Keplerian anti-cyclonic rotation Ro = 4/3. This scaling is indicative of a disturbance which grows purely under the Orr mechanism, and has been shown to hold for oblique disturbances in planar flows by Chapman (2002) and Heaton & Peake (2007). The Re dependence is shown in Figure 4.7, though larger Re than shown in Figure 4.2 are needed to resolve the scaling. An example of the veloc- ity field is displayed in Figure 4.8 for G = 11.96 at (Ro,Re) = (−20, 900). Since the optimal is two-dimensional, with w = 0, a vector plot representation is used in Figure 4.8. Figure 4.7: Left: G in Re at Ro = −20. The dashed line depicts the scaling G ∝ Re2/3. Right: wavenumbers αmax (blue) and βmax (red) for Ro = −20. Figure 4.8: (x, y)-profiles of the optimal disturbance at (Ro,Re) = (−20, 900) with (αmax, βmax) = (1.23, 0). The velocity at t = 0 (left) and t = Tmax = 8.34 (right) is shown. Each velocity is represented as a vector plot of the streamwise and wall-normal velocities (u, v). We conclude that for Ro far enough from zero, optimal perturbations can no longer be amplified in an Orr and lift-up mechanism combination and instead rely solely upon the Orr mechanism. As Ro decreases, we find that G monotonically decreases towards the 4.2 Optimal growth in RPCF 133 two-dimensional optimals of the Orr mechanism. In Figure 4.9 we plot G and its cor- responding optimal wavenumbers over −20 < Ro < 0 at Re = 500. G drops sharply from G = 296.02 to G = 9.547 over −2 < Ro < 0, beyond which it decreases much more gradually until it reaches a minimum value of G = 7.551 when Ro ≈ −14.46. To gain in- sight into the structure of the optimal perturbations as Ro decreases, we plot the energy fractions against Ro of the optimal disturbances at t = 0 and t = Tmax in Figure 4.10. In the Figure it can be seen that non-rotating type optimals, with energy transfer from v and w to u at t = Tmax, survive for only small Ro. The energy fraction Eu for the optimal at Tmax falls rapidly until Ro ≈ 2 when it reaches a turning point and begins to steadily increase. Over the same region Ew swiftly grows before beginning a steady decrease at the same turning point. From −13.82 < Ro < −2, the optimals transfer energy from an energetically streamwise dominated initial disturbance to an energetically spanwise dominated disturbance. Over −14.46 < Ro < −13.82, the greatest energy fraction of the perturbation at Tmax is in the streamwise velocity, with Eu becoming greater than Ew. We examine the form of the optimal disturbance at (Ro,Re) = (−0.2, 500) for which there is little streamwise energy transfer as Eu(0) = 0.416 and Eu(Tmax) = 0.463. The optimal wavenumbers in this case are (αmax, βmax) = (0.89, 0.81) and their closeness al- lows for the true aspect ratio of the velocity fields to be accurately depicted and easily viewed. In Figure 4.12 we plot the evolution of G(t) in blue, and in red we plot the evolution of the specific initial condition which maximizes G in time. In Figures 4.13 and 4.14 we plot cross-sections of the perturbation velocities in the (x, y) and (z, y) planes respectively, with each plot at a point in time corresponding to the black dots in Figure 4.12. The counter-rotating vortices in Figures 4.13 and 4.14 are seen to change their orientation in the streamwise and spanwise directions as time evolves. This is indication of a strong dependence on the Orr mechanism for growth at the relatively minor devi- ation from the non-rotating case, Ro = −0.2. For Ro < −14.46 our algorithm fails to find a three-dimensional disturbance with greater amplification that the two-dimensional Orr-type optimal, resulting in a discontinuity in the curves of αmax and βmax. In Figure 4.11, a refined surface plot of Gmax in the wavenumber domain reveals two peaks, one for β = 0 and one for β small but non-zero, each at slightly separate α. As rotation is further decreased, the peak with β 6= 0 decays whilst the β = 0 is constant. In this way there is a discontinuous transition between three-dimensional and two-dimensional optimals. That large rotation forces the optimal to be two-dimensional can be viewed as a transient version of the Proudman theorem (Proudman (1916)) on the two-dimensionalization of steady, inviscid fluids dominated by Coriolis rotation. It should be emphasized that since the Orr and lift-up mechanisms are not unique to ro- 134 4. Non-modal growth in rotating plane Couette flow tating Couette flow, we find that no new growth mechanisms strong enough to maximize G emerge under cyclonic rotation. The effect of the introduction of cyclonic rotation on the flow can therefore be viewed as being the reversal of the dominance of the lift-up mechanism. Figure 4.9: G for Re = 500. Left: G against Ro. Right: Optimal wavenumbers αmax and βmax against Ro. The transition from lift-up mechanism dominated optimals to Orr mechanism optimals is evident as Ro decreases. Figure 4.10: Energy fractions of optimal disturbances at t = 0 (left) and t = Tmax (right) for Re = 500. Eu, Ev and Ew represent the fraction of total perturbation energy in u, v and w respectively. 4.2 Optimal growth in RPCF 135 Figure 4.11: Colour-plot of Gmax in the (α, β)-plane for Ro = −14 (left) and Ro = −14.8 (right). The switch between G with β 6= 0 and β = 0 can be seen from one panel to the next. The discontinuity in Figures 4.9 and 4.10 is a result when this switch, occuring between −15 < Ro < −14. Figure 4.12: Gain curve for Re = 500, Ro = −0.2, α = αmax = 0.89, β = βmax = 0.81. The blue curve is the gain G(t), which is the maximum growth at time t for an arbitrary initial condition. The red curve is the evolution of the initial condition which maximizes the growth at t = Tmax = 11.66. The black dots indicate the times of the velocity fields shown in Figures 4.13 and 4.14 136 4. Non-modal growth in rotating plane Couette flow Figure 4.13: Perturbation velocity fields for Re = 500, Ro = −0.2, α = αmax = 0.89, β = βmax = 0.81. Each image is a vector plot of the streamwise and wall-normal velocities u and v, and a colour-plot of the spanwise velocity w, with red and blue indicating positive and negative, in the (x, y)-plane at z = 0. The velocity vectors clearly ‘swing’ in the x direction as time evolves. Figure 4.14: Perturbation velocity fields for Re = 500, Ro = −0.2, α = αmax = 0.89, β = βmax = 0.81. The fields are represented by vector plots of the spanwise and wall- normal velocities w and v, and colour-plots of the streamwise velocity u, with red and blue indicating positive and negative, in the (z, y)-plane at x = 0. The velocity vectors ‘swing’ in the z direction as time evolves. 4.2 Optimal growth in RPCF 137 4.2.2 Subcritical anti-cyclonic rotation Figure 4.15: Optimal growth in the subcritical anti-cyclonic rotation regime. The black curves are contours of G = 100.125, 100.25, 100.375, . . . , 104.5 The solid red curve represents the energy stability boundary. Our calculations of optimal perturbations in the anti-cyclonic regime repeat the trend of §4.2.1 insofar as G drops rapidly as Ro is increased from Ro = 1. However, unlike cyclonic rotation, the addition of anti-cyclonic rotation introduces a linear growth mechanism not present in the non-rotating flow. Optimals for Ro = 1 are amplified by the anti lift-up mechanism identified by Antkowiak & Brancher (2007). The effect is named “anti lift-up” because high energy streamwise vortices are generated from an initial streamwise oriented streak; a process opposite to the lift-up mechanism. Rincon et al. (2007) explained the analogy between the two growth mechanisms by analysis of the linearized equations governing the streamwise velocity u and the streamwise vorticity ξ = ∂w∂y − ∂v ∂z , under the assumption of streamwise independence, ∂u ∂t = (Ro− 1)v + 1 Re ∇2u, (4.21a) ∂ξ ∂t = −Ro ∂u ∂z + 1 Re ∇2ξ, (4.21b) where we have used that U ′ = 1 for the linear shear profile. For Ro = 0 the wall-normal velocity v of the streamwise vorticity field ξ acts as a forcing on the streamwise velocity u, driving a streak. When Ro = 1, the ∂u∂z term acts as a forcing on ξ, driving a streamwise 138 4. Non-modal growth in rotating plane Couette flow vortex field. The scaling G ∝ Re2 holds for anti lift-up optimals, as it did for lift-up optimals, and is shown in Figure 4.16. In optimizing for G, we find that the optimal dis- turbances for Ro = 1 take αmax = 0, which is consistent with the calculations of Yecko (2004). The Ro = 1 case is therefore in contrast to the non-rotating optimals for which αmax is small but non-zero. As Ro increases, the anti lift-up effect loses influence on the optimals and beyond a critical Ro, the optimal perturbations have two-dimensionalized and grow via the Orr mechanism. Following our analysis for cyclonic rotation, we fix Re = 500 and in Figure 4.17 the evolution of G and the optimal wavenumbers are plot- ted as Ro increases. G decays rapidly from G = 291.7 to G = 16.52 over the range 1 < Ro < 1.05, whereupon it slowly decays until Ro = 1.276 when G = 7.55 and the Ro- independent Orr mechanism dominates. The optimal wavenumbers in the right panel of Figure 4.17 and the energy fractions in Figure 4.18 give a sense of the form of the optimal disturbance at each Ro. The spanwise wavenumber βmax increases from βmax = 1.658 to βmax = 1.777 over 1 ≤ Ro ≤ 1.018, while αmax remains zero, and from Figure 4.18 we see that energy is transfered from u to v and w indicating that the anti lift-up mechanism is preferred over this Ro range. Between Ro = 1.018 and Ro = 1.275 the optimals are oblique, with αmax increasing and βmax decreasing for greater rotation. The perturba- tions behave in the same way as the oblique cyclonic optimal depicted in Figures 4.12, 4.13 and 4.14, though now they are amplified in a mixture of the anti lift-up and the Orr mechanisms. The spanwise energy fraction Ew of the initial disturbance reaches a maximum at Ro ≈ 1.15, after which it decreases, reaching zero at the onset of two- dimensionalization. The transition from three-dimensional to two-dimensional optimals is noticeably smoother than the transition for cyclonic rotation in Figures 4.9 and 4.10, and occurs over a much shorter range of Ro. Thus the effect of anti-cyclonic rotation on optimal linear energy growth mechanisms is the introduction and dominance of the anti lift-up mechanism at Ro = 1, and the sub- sequent erosion in strength of three-dimensional growth mechanisms towards dominance of the Orr mechanism as Ro is increased beyond unity. 4.2 Optimal growth in RPCF 139 Figure 4.16: Left: Variation of G with Re for Ro = 1. The dashed line depicts the scaling G ∝ Re2. Right: wavenumbers αmax (blue) and βmax (red) for Ro = 0. Figure 4.17: Variation of G for Re = 500. Left: G against Ro. Right: Optimal wavenum- bers αmax and βmax against Ro. The transition from lift-up mechanism dominated opti- mals to Orr mechanism optimals is evident as Ro decreases. 140 4. Non-modal growth in rotating plane Couette flow Figure 4.18: Energy fractions of optimal disturbances at t = 0 (left) and t = Tmax (right) for Re = 500. Eu, Ev and Ew represent the fraction of total perturbation energy in u, v and w respectively. 4.2.3 Rotational independence of the Orr mechanism The two-dimensional optimals observed for strong cyclonic and anti-cyclonic rotation in §4.2.1 and §4.2.2 appear to be independent of rotation. We now show that because they are spanwise independent they are therefore necessarily independent of Ro. The spanwise independent perturbation equations reduce to ∂u˜ ∂t + U ∂u˜ ∂x + v˜U ′ = − ∂p˜ ∂x + 1 Re ∇2u˜+Ro v˜, (4.22a) ∂v˜ ∂t + U ∂v˜ ∂x = − ∂p˜ ∂y + 1 Re ∇2v˜ −Ro u˜, (4.22b) ∂w˜ ∂t + U ∂w˜ ∂x = 1 Re ∇2w˜, (4.22c) ∂u˜ ∂x + ∂v˜ ∂y = 0, (4.22d) where ′ = ddy and U = y. Clearly the spanwise velocity w˜ has decoupled from u˜ and v˜; it can therefore play no role in an energy growth mechanism and it is henceforth ignored. Assuming wave-like dependence [u˜, v˜, p˜](x, y) = [u, v, p](y)eωt+iαx, (4.23) 4.2 Optimal growth in RPCF 141 and then using equation (4.22d) to eliminate u, equations (4.22a) and (4.22b) can be combined to eliminate p˜ and to derive a non-constant coefficient fourth order ODE for v v′′′′ + A(y)v′′ +B(y)v = 0, (4.24) where A(y) = Re(iαy − ω)− 2α2, (4.25a) B(y) = α4 − α2Re(iαy − ω). (4.25b) From the no-slip condition applied to the channel walls, v must satisfy v = v′ = 0 on y = ±1. (4.26) We note in passing that a simpler equation governing the spanwise vorticity could be derived but it is not clear what boundary conditions should be applied to it, since infor- mation would be needed concerning ∂u∂y on y = ±1. Equation (4.24) is independent of Ro, hence the evolution of any two-dimensional perturbation (u, v) will not be affected by the introduction of system rotation. To highlight the effect of this on optimal growth, we re-print Figure 4.2 with added contours of Gmax for the two-dimensional Orr mechanism growth in Figure 4.19. These contours are calculated by constraining the optimization over wavenumbers, described in §4.1.1, such that β = 0. In the Figure, it can clearly be seen that G tends towards the Orr mechanism growth for |Ro| large. Figure 4.19: Contours of G in (Ro,Re), as in Figure 4.2 with Ro-independent Orr mech- anism growth marked in dashed blue. 142 4. Non-modal growth in rotating plane Couette flow 4.3 Discussion Summarizing the results of §4.2, we find that the effect of strong cyclonic or anti-cyclonic rotation is to switch the dominance of the lift-up and anti lift-up energy growth mecha- nisms to the two-dimensional Orr mechanism. We now discuss what relevance this may have on transition to turbulence in flows with strong rotation. However, we must be careful in drawing conclusions from the behaviour of transients alone. Two recent ap- proaches have been successful in explaining aspects of the transitional process in shear flows: nonlinear optimization over initial conditions, and analysis of finite-amplitude nu- merical solutions of the Navier-Stokes equations. It has been shown, using a variational technique over the fully nonlinear Navier-Stokes equations, that transition is not efficiently triggered by finite-amplitude counterparts of the linear transients, but rather by inherently nonlinear states which undergo some lo- calization during their energy growth (Rabin et al. (2012), Pringle et al. (2012)). These structures have distinct growth phases, using combinations of the lift-up mechanism and the Orr mechanism to drive their energy amplification. Thus we can, tentatively, con- jecture that in the absence of one of the linear growth mechanisms the structure of the minimal seed would be much changed, if such a seed were found to exist at all. A new route to transition would be required, one that relies solely on the Orr mechanism. Our results however do not preclude the existence of strong linear growth mechanisms about a modified base flow, though how the base flow would be initially modified would remain in question. Exact coherent structures, or finite-amplitude solutions to the Navier-Stokes equations, are conjectured to be key components to an understanding of turbulent dynamics in sub- critical shear flows. Exact coherent structures are thought to form an invariant set about which chaotic dynamics are supported, and chaotic dynamics are thought to be related to turbulence. There has been recent success in finding exact coherent structures in plane Couette and plane Poiseuille flows (see Nagata (1990, 2013), Waleffe (1998, 2003), Gibson et al. (2009), Itano & Generalis (2009)). Furthermore, exact coherent structures are thought to play a role in what is termed a self-sustaining process (Waleffe (1995, 1997, 1998)), which drives and maintains turbulent dynamics in shear flows. The self- sustaining process consists of three phases: weak streamwise vortices re-distribute the streamwise momentum to create spanwise fluctuations in the mean streamwise velocity; spanwise inflections in the streamwise velocity cause a three-dimensional streamwise pe- riodic instability to develop; and nonlinear interaction driven by the three-dimensional instability re-energizes the streamwise vortices such that the process can be repeated. 4.3 Discussion 143 The lift-up mechanism is thought to be crucial in the first step of the cycle, as energy is transfered from streamwise vortices and amplified into streamwise streaks. Exact coher- ent structures are thought to result from this action, and then lose stability to streamwise periodic disturbances, as per the second phase of the self-sustaining process. Hence, our results show that in the presence of strong rotation such a cycle would not be possible unless streamwise momentum could be re-distributed, and an exact coherent structure formed, in the absence of the lift-up mechanism. Rincon et al. (2007) attempted to con- tinue exact coherent structures from plane Couette flow into the cyclonic and subcritical anti-cyclonic rotation regimes. They found that structures could be continued into the cyclonic regime for very small rotation rates, however they could not find any structures which could be continued into the anti-cyclonic rotation regime for Re < 500. If it is the case that the lift-up mechanism is crucial to the generation of exact coherent states, then our results suggest that no such structures would be found for subcritical anti-cyclonic ro- tation and strong cyclonic rotation. Whether the anti lift-up mechanism could be used to generate an exact coherent structure of another form remains unanswered. We conclude by noting that Lesur & Longaretti (2005) observed transition in numerical simulations for (Re,Ro) = (3000, 1), offering some hope that coherent structures might indeed exist for high Reynolds numbers. 144 4. Non-modal growth in rotating plane Couette flow Chapter 5 Nonnormality and optimal growth of the Papkovitch-Fadle operator In this chapter we address the maximum spatial energy amplification that can be achieved by an arbitrary edge condition expressed as a finite series of Papkovitch-Fadle eigenfunc- tions. The layout of the chapter is as follows. In §5.1 we state the governing equations and give their formal solution. Our optimal growth framework is defined in §5.2 and we present computations of the optimal disturbances, for both an unconstrained disturbance and a disturbance with constrained corner shear-stress. In §5.3 we conduct an asymp- totic analysis of how the energy inner product behaves with respect to the interaction of long and short wavelength modes. In §5.4 we extend our optimal growth analysis to spatially developing, two-dimensional plane Poiseuille flow. Finally in §5.5 we contrast our findings for spatially developing optimal disturbances to temporal optimals through an asymptotic analysis of the energy inner product in temporal plane Couette flow. 5.1 Governing equations and modal solution We consider the solution of the biharmonic operator problem on a semi-infinite Cartesian domain Ω ∆2Ψ = 0, Ω = {x, y : −1 ≤ y ≤ 1, 0 ≤ x ≤ ∞} (5.1a) ∂Ψ ∂x (x,±1) = ∂Ψ ∂y (x,±1) = 0, (5.1b) ∂Ψ ∂x , ∂Ψ ∂y → 0, as x→∞ (5.1c) where ∆ is the Laplacian operator ∆ = ∂ 2 ∂x2 + ∂2 ∂y2 . This is also known as the Papkovitch- Fadle operator. The boundary conditions (5.1b) and (5.1c) are motivated by conditions of no-slip on the walls in the Stokes edge problem (see Figure 5.1), for which Ψ can be 145 146 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator considered the streamfunction. A formal solution can be written as Ψ(x, y) = ∑ n cnψn(y)e iknx, (5.2) for coefficients cn which specify the conditions at the edge Γ = {x = 0,−1 ≤ y ≤ 1}. The eigenfunctions ψn(y) are the Papkovitch-Fadle eigenfunctions, and they can be even or odd in y. The even eigenfunctions are written as ψn(y) = ( y − 1 ) sinh [ kn ( y + 1 )] + ( y + 1 ) sinh [ kn ( y − 1 )] , (5.3a) kn + sinh kn cosh kn = 0, kn ∈ C, (5.3b) and for odd eigenfunctions we have ψn(y) = ( y − 1 ) sinh [ kn ( y + 1 )] − ( y + 1 ) sinh [ kn ( y − 1 )] , (5.4a) kn − sinh kn cosh kn = 0, kn ∈ C. (5.4b) The dispersion relations (5.3b) and (5.4b) are transcendental equations which we solve by Newton iteration. The solutions occur as reflections in each quadrant of the complex plane k = ±kr + iki, k = ±kr − iki, (5.5) for kr, ki ∈ R. We only include modes for which Im{k} > 0 in the solution expansion, to ensure that (5.1c) is satisfied. A selection of the modes are plotted in Figure 5.2. To describe each mode, we use the notation kn = kr + iki, k−n = −kr + iki so that we can list the eigenvalues such that Im{k1} < Im{k2} < . . . . (5.6) In addition, neglecting modes with Im{k} < 0 renders the system parabolized in x. Therefore we can treat the Papkovitch-Fadle operator (5.1) as an evolution operator. 5.1 Governing equations and modal solution 147 Figure 5.1: Stokes edge problem: The Papkovitch-Fadle operator (5.1) can be thought of as a mathematical description of a spatially developing, two-dimensional, viscous fluid. The fluid is driven by an arbitrary upstream (x < 0) stirring, which creates an edge profile at the channel entrance. The Papkovitch-Fadle operator governs the downstream (x > 0) development of the fluid along the channel. Figure 5.2: Eigenvalue spectrum in the complex k-plane. Blue and red dots denote the even and odd modes respectively. 148 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.2 Nonnormality and optimal energy growth The Papkovitch-Fadle eigenfunctions are known to exhibit nonnormal characteristics. We investigate this by first introducing a norm based on the kinetic energy density of a fluid in Stokes flow. Using the expressions for the fluid velocities u = Ψy = ∑ n anun(y)e iknx = ∑ n cnψ ′ n(y)e iknx, (5.7a) v = −Ψx = ∑ n anvn(y)e iknx = − ∑ n ikncnψn(y)e iknx, (5.7b) and un = (un, vn), we define the energy inner product Emn = 〈un,um〉 = 1 2 ∫ 1 −1 u∗mun + v ∗ mvn dy, (5.8) equivalently Emn = 〈un,um〉 = 1 2 ∫ 1 −1 (ψm) ∗ψn + 1 k∗mkn (ψ′m) ∗ψ′n dy. (5.9) The Papkovitch-Fadle eigenfunctions are non-orthogonal with respect to this inner prod- uct. More precisely, odd modes are non-orthogonal to other odd modes, and similarly for even modes. Odd and even modes are orthogonal with respect to each other. If we let pn and qn denote odd and even eigenfunctions from (5.4a) and (5.3a) respectively, then we have the following orthogonality relations 〈pn,pm〉 = 4 sinhα α3 ( coshα− cosh β ) − 4 α2 ( 1− coshα cosh β ) + cosh 2kn k∗m ( 1 α − 1 β ) + cosh 2k∗m kn ( 1 α + 1 β ) , m 6= −n (5.10a) 〈qn, qm〉 = 4 sinhα α3 ( coshα + cosh β ) − 4 α2 ( 1 + coshα cosh β ) + cosh 2k∗m k∗m ( 1 α − 1 β ) + cosh 2kn kn ( 1 α + 1 β ) , m 6= −n (5.10b) 〈pn, qm〉 = 0, ∀m,n (5.10c) where α = kn + k∗m and β = kn − k ∗ m. If km = k−n we have α = 0, in which case 〈pn,p−n〉 = 4 3 ( 2 + cosh 2kn ) − 2 sinh 2kn kn , (5.11a) 〈qn, q−n〉 = 4 3 ( 2− cosh 2kn ) − 2 sinh 2kn kn . (5.11b) 5.2 Nonnormality and optimal energy growth 149 The non-orthogonality of the eigenfunctions is a consequence of the nonnormality of the Papkovitch-Fadle operator. Transient dynamics are an important feature of a nonnormal system and here we investigate the maximum energy amplification that the fluid can experience and the form of the optimal edge condition. A velocity u spanned by the first 2N modes is expressed as u(x, y) = N∑ n 6=0, n=−N anu˜n(y)eiknx, (5.12) where 〈u˜n, u˜n〉 = 1. In this notation u˜n can be either p˜n or q˜n depending on whether n corresponds to an odd or even mode. Using the expressions a = (a−N , . . . , aN), a ∈ C2N (5.13) Mmn = 〈u˜n, u˜m〉, M ∈ C2N×2N (5.14) Λ =     eik−Nx . . . eikNx     , Λ ∈ C 2N×2N (5.15) the energy norm is defined as E(x) = ‖u(x)‖2E = 1 2 ∫ 1 −1 uHu dy = aHΛHMΛa, (5.16) where H denotes the Hermitian transpose. The physical, kinetic energy-density E˜ over the semi-infinite strip is related to the energy norm via the relation E˜ = 1 2 ∫ ∞ 0 ∫ 1 −1 uHu dy dx = ∫ ∞ 0 ‖u(x)‖2E dx. (5.17) Thus E˜ is always constant, as the flow is steady. For a given streamwise location x we compute the maximum energy gain G(x) = max u E(x) E(0) = max u0 6=0 ‖u(x)‖2E ‖u0‖2E . (5.18) To determine G(x), we formulate a Lagrangian optimization problem for a. The La- grangian is L = aHPa− λ(aHMa− 1), (5.19) 150 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator where P = ΛHMΛ. Taking variations with respect to the Lagrange multipliers, we have δL δλ = aHMa− 1, (5.20) δL δaH = Pa− λMa, (5.21) δL δa = aHP − λaHM. (5.22) Setting (5.20) to zero fixes E(0) = 1. Equating (5.21) or (5.22) to zero gives an eigenvalue problem for right eigenvector a or left eigenvector aH respectively. We solve Pa = λMa, (5.23) for the largest eigenvalue λ and its eigenvector a. λ corresponds to the maximum, or optimal, gain G achieved at a point x, and a is the coefficient vector which describes the optimal edge velocity. 5.2.1 Optimal growth and convergence We begin by considering the optimal growth characteristics for a modal truncation with no restrictions on the coefficient vector a. Given that even and odd modes are orthogonal, we restrict our attention to even modes since calculations show them to be responsible for higher gains. For clarity we list the modes from m = 1 . . .M rather than from n = −N . . .N using m(n) =    2|n| − 1, n > 0 2|n|, n < 0. The absence of conditions on a means that there is no guarantee that any quantities at the edge Γ = {x = 0,−1 ≤ y ≤ 1}, which are defined by a series in a, will be well-behaved. However, since the computations in this work use a finite expansion of eigenfunctions, we are not concerned with the convergence of an infinite eigenfunction series at the edge or elsewhere. The quantities of interest are the edge velocities and streamfunction u0 = M∑ m=1 iam km ψ˜ ′ m(y), v0 = M∑ m=1 amψ˜m(y), Ψ0 = M∑ m=1 iam km ψ˜m(y), (5.24) where ′ = ddy and ψ˜m(y) is normalized such that Emm, as defined in equation (5.9), is unity. We find that the gain G(x) converges slowly as a function of modal truncation M . Results are displayed in Figures 5.3 and 5.4. As M increases, the magnitude of the optimal coefficients, |am|, become large for m ≈ M . This suggests that the growth mechanism 5.2 Nonnormality and optimal energy growth 151 is reliant on interactions between high wavenumber modes. We address what drives this effect in §5.3. To quantitatively investigate the rate of convergence of G, we fix x and compute a finite-difference derivative approximation δGδM , where δG δM = G(M + 1)−G(M) M + 1−M = G(M + 1)−G(M), M ∈ N. (5.25) Results are plotted in Figure 5.4. δGδM is dependent on x, however for x > 0.1 we have the approximate scaling δGδM ∝ M −2. At the trunction M = 104, G is reasonably well converged for x ≥ 0.1 and less well converged as x approaches zero. Based on our calculations it is unclear whether G will converge in the limit x → 0. Convergence may require a greater truncation as x approaches zero, or a point may be reached where the gain does not converge at all. Trefethen (1997) stated that the numerical abscissa ω, which can be thought of as the initial gradient of G(x), is unbounded: ω = lim x→0 dG dx =∞. (5.26) In Figure 5.6 we show that ω scales with the modal expansion number M , implying an in- finite eigenfunction expansion would be required to capture the limit ω → 0. Nonetheless, the steep inital gradient of G(x) is visible in the right-panel of Figure 5.3 for M = 103. Trefethen (1997) and Trefethen & Embree (2005) also noted that the repeated oscilla- tions of the gain curve in Figure 5.3 correspond to a sequence of counter-rotating Moffatt vortices in Stokes flow. Figure 5.3: On the left we plot G(x) (solid line) for M = 103. The dashed line is the energy decay rate of the leading mode. On the right is a blow-up for low x with M = 10, 102 and 103 corresponding to the blue, red and black lines respectively. 152 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.4: G and δGδM against M for x = 10 −3, 0.01, 0.05. G is not well converged for M ≤ 104, for the values of x considered here. Figure 5.5: G and δGδM against M for x = 0.1, 0.15, 0.2. The dashed line in the right panel is the slope M−2. G reaches satisfactory convergence for M = 104. 5.2 Nonnormality and optimal energy growth 153 Figure 5.6: On the left we plot |am| at four different truncations M for x = 0.1. On the right is the numerical abscissa ω against mode truncation M . Figure 5.7: Edge profiles in y optimising G at x = 0.1. Along each row from the top are u0, v0 and Ψ0 and the columns are the truncations M = 10, 102, 103 and 104. 154 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.8: Blow-up of the lower wall region of the edge profiles in Figure 5.7. The top and bottom rows are u0 and v0, and the columns are the truncations M = 10, 102, 103 and 104. Figure 5.9: Optimal profiles in y at x = 0.1. Along each row from the top are u(0.1, y), v(0.1, y) and Ψ(0.1, y) and the columns are the truncations M = 10, 102, 103 and 104. 5.2 Nonnormality and optimal energy growth 155 In Figures 5.7 and 5.9 we plot the wall-normal profiles of u, v and Ψ which optimize the energy gain at x = 0.1. We plot the profiles both at the edge Γ and their form as they have evolved to x = 0.1. While Ψ is well resolved throughout the domain, u0 and v0 display sharp oscillatory behaviour near the walls (y = ±1). This is a consequence of the large amplitude, large wavenumber modes selected by the energy growth mechanism. Increasing the truncation appears to dampen the near-wall oscillations in v0 but leads to increasingly sharp behaviour in u0 at the walls, though a blow-up of the lower wall region in Figure 5.8 confirms that the eigenfuctions are always continuous. The sharp near-wall gradient in u0 suggests that the optimal edge profile tends towards a slip con- dition in streamwise velocity as the modal resolution is increased. The no-slip condition at the edge, Ψy(0,±1) = 0, can be relaxed to include a slip profile at the edge, though this has no effect on our findings. Further downstream, the large wavenumber modes become less prominent in the profiles at x = 0.1 as they are very highly damped, thus the eigenfunctions have no oscillatory behaviour near the walls. The optimal velocities and streamfunction are plotted in the (x, y)-plane in Figure 5.10. From the Figure, steady Moffatt vortices are seen to emanate from the corners of the domain. We also plot the streamfunction on a log scale in Figure 5.11, to allow the succession of Moffatt vortices which cause the oscillations in the gain curve of Figure 5.3 to be seen. 156 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.10: Shading plots of Ψ (top), u (bottom-left) and v (bottom-right) for x = 0.1. The plots do not have a true aspect ratio, so that the Moffatt vortices can be clearly seen. Figure 5.11: Shading plot of | log(Ψ)|. Plotting on the log-scale allows the repeated oscillations of Moffatt eddies to be visualized. 5.2 Nonnormality and optimal energy growth 157 5.2.2 Optimal growth with bounded wall shear It is found that as the modal resolution of the optimal edge condition found in §5.2.1 is improved, the fluid exerts increasingly high stress on the channel corners. We define τ to be the corner shear stress for normalized viscosity µ = 1 τ = µ ∂u ∂y ∣ ∣ ∣ ∣ y=±1,x=0 = ∂2Ψ ∂2y ∣ ∣ ∣ ∣ y=±1,x=0 . (5.27) Figure 5.12 suggests that |τ | scales linearly with the eigenfunction expansion coefficient M , which would imply unbounded stress if an infinite eigenfunction expansion were to be considered. This may be physically questionable, therefore we are motivated to seek a solution Ψτ , for which τ is bounded at the corner. The Lagrangian optimization problem in equation (5.19) is re-formulated such that |τ | is normalized to unity. Writing aτ for the coefficient vector of the constrained solution, we impose the constraint |τ |2 = aHτ Baτ = a H τ b Hbaτ = 1, (5.28) where bn = iψ˜′′n(±1) kn = 4i cosh2 kn √ 〈qn, qn〉 . (5.29) The Lagrangian now becomes Lτ = a H τ Paτ − λ(a H τ Maτ − 1)− µ(a H τ Baτ − 1). (5.30) Setting variations of Lτ to zero yields the double-eigenvalue problem Paτ = λMaτ + µBaτ . (5.31) Techniques are available for the direct solution of (5.31), such as Blum & Chang (1978) and Ji (1996), however, for its ease of implementation we solve for aτ and the energy gain G(x) using the active-set constrained nonlinear optimization routine fmincon in MATLAB. The routine makes an approximation for the Hessian of the Lagrangian, which is used to define a search direction along which a line search can be conducted. The addition of the new constraint renders the calculations more computationally expensive than in §5.2.1 and owing to this we can calculate a solution as far as truncation number M = 103. The solutions share qualitiative features with the unconstrained optimals (cf Figures 5.7 and 5.9 with 5.13), suggesting that the growth mechanism is not reliant on large shear stress in the corners of the domain. In Figure 5.12 we note that Ψτ experiences less energy amplification than the unconstrained optimal, however the difference narrows as truncation number M is increased. 158 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.12: On the left: Corner stress |τ | against expansion number M for optimals at x = 0.1. The red dashed curve indicates the slope with gradient one. On the right: G against expansion number M for x = 0.1. The blue and red lines are the constrained and unconstrained optimals respectively. Figure 5.13: Wall-normal profiles of the constrained optimal which maximizes G(x) at x = 0.1. The top row has the edge profiles u0, v0 and Ψ0. On the second row we plot u(0.1, y), v(0.1, y) and Ψ(0.1, y). 5.3 Short wavelength asymptotics 159 5.3 Short wavelength asymptotics In this section we derive asymptotic expressions for short wavelength modes, or equiv- alently, modes with large wavenumber. Our attention is restriced to even modes since they are found to be responsible for the largest transient growth in §5.2.1. The aim of this analysis is to investigate the energy inner product interactions of modes in the short wavelength limit, in order to provide an insight into why the optimal edge coefficients utilize the non-orthogonality of large wavenumber modes. We start by splitting the dispersion relation into real and imaginary parts. For a mode k we have k + sinh k cosh k = 0. (5.32) Let λ = 2k. Then λ+ sinhλ = 0, (5.33) =⇒ (λr + sinhλr cosλi) + i(λi + coshλr sinλi) = 0. (5.34) Hence the following two relations must be satisfied sinhλr = − ( 1 cosλi ) λr, (5.35) sinλi = − ( 1 coshλr ) λi. (5.36) In order for equation (5.35) to be satisfied for large |λ|, cosλi must be small and negative. Therefore we write λi = µ+ ε1 + ε2 + . . . , for µ = ( 2n− 1 2 ) pi, n ∈ N, (5.37) for small parameters ε1 and ε2 with ε2 < ε1. Equation (5.36) becomes − sin ( µ+ ε1 + ε2 + . . . ) = − 1 coshλr ( µ+ ε1 + ε2 + . . . ) . (5.38) Since sinµ = −1 and cosµ = 0, we have coshλr ( 1− 1 2 ( ε1 + ε2) 2 + . . . ) = µ+ ε1 + ε2 + . . . , (5.39) 160 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator =⇒ coshλr = ( µ+ ε1 + ε2 + . . . )( 1− 1 2 ( ε1 + ε2) 2 + . . . )−1 , (5.40) = ( µ+ ε1 + ε2 + . . . )( 1 + 1 2 ( ε1 + ε2) 2 + . . . ) , (5.41) = µ+ ε1 + ε2 + µε21 2 + 2µε1ε2 + µε22 2 + . . . , (5.42) = µ+ δr + . . . . (5.43) Let ε1 + ε2 = δi. From equation (5.35) we get ( (µ+ δr + . . .) 2 − 1 ) 1 2 cos ( µ+ δi + . . . ) = − cosh−1 ( µ+ δr + . . . ) , (5.44) µ ( 1 + 2δr µ − 1− δ2r µ2 + . . . ) 1 2 sin ( δi + . . . ) = − log ( µ+ δr + ( (µ+ δr + . . .) 2− 1 ) 1 2 + . . . ) , (5.45) δi + . . . = − 1 µ ( log ( 2µ+ 2δr + . . . ) ) , (5.46) = − log 2µ µ − δr µ2 + . . . , (5.47) =⇒ ε1 + ε2 + . . . = − log 2µ µ − 1 µ2 ( ε1 + ε2 + µε21 2 + 2µε1ε2 + µε22 2 + . . . ) . (5.48) From this we conclude ε1 = − log 2µ µ , (5.49) ε2 = − ε21 2µ = − (log 2µ)2 2µ3 . (5.50) Thus we have that λi = µ− log 2µ µ +O ( (log µ)2 µ3 ) , (5.51) and from equation (5.42) λr = log 2µ+ (log 2µ)2 2µ2 +O ( log µ µ2 ) . (5.52) 5.3 Short wavelength asymptotics 161 Now we can determine how the energy inner product for even modes behaves for |k| → ∞. Recalling equation (5.10b) from §5.1, we have Emn = 4 sinhα α3 ( coshα + cosh β ) − 4 α2 ( 1 + coshα cosh β ) + cosh 2k∗m k∗m ( 1 α − 1 β ) + cosh 2kn kn ( 1 α + 1 β ) , m 6= −n (5.53) where α = kn + k∗m and β = kn − k ∗ m. We investigate modal interactions with respect to the normalized energy inner product. 5.3.1 Normalization factor First we find an expression for the normalization factor to leading order E2n = Enn, n→∞. (5.54) Suppose kn = 1 2 ( log 2µ+ (log 2µ)2 2µ2 +O ( log µ µ2 ) + i ( µ− log 2µ µ +O ( (log µ)2 µ3 ))) . (5.55) Then α = 2(kn)r = λr and β = i2(kn)i = iλi. Therefore coshα = sinhα = µ+ (log 2µ)2 2µ +O ( log µ µ ) , (5.56) cosh β = − log 2µ µ +O ( (log µ)2 µ3 ) . (5.57) So for the energy norm we have E2n = ( 4µ2 (log 2µ)3 +O ( 1 log µ )) + ( 4 log 2µ +O ( 1 (log µ)2 )) + ( − 2 log 2µ +O ( log µ µ2 )) + ( − 2 log 2µ +O ( log µ µ2 )) . (5.58) Hence E2n = 4µ2 (log 2µ)3 +O ( 1 log µ ) , n→∞ (5.59) where µ = (2n− 12)pi. 162 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.3.2 Interaction of short wavelength modes Now consider the leading order behaviour of Emn for two short wavelength modes, i.e. with |kn| and |km| large, but not equal. We investigate the interaction of modes with m = n− a, with a ∈ N, a = O(1), n→∞. (5.60) Expressions for wavenumbers kn and km are kn = log 2µ 2 + api µ + (log 2µ)2 4µ2 +O ( log µ µ2 ) +i ( µ 2 +api− log 2µ 2µ +O ( (log µ)2 µ3 )) , (5.61) km = log 2µ 2 + (log 2µ)2 4µ2 +O ( log µ µ2 ) + i ( µ 2 − log 2µ 2µ +O ( (log µ)2 µ3 )) . (5.62) Therefore we have α = kn + k ∗ m, (5.63) = log 2µ+ api µ + (log 2µ)2 2µ2 + iapi +O ( log µ µ2 ) . (5.64) β = kn − k ∗ m, (5.65) = api µ + i ( µ+ api − log 2µ µ ) +O ( (log µ)2 µ3 ) . (5.66) coshα = (−1)a coshλr, (5.67) = (−1)a ( µ+ (log 2µ)2 2µ +O ( log µ µ )) . (5.68) sinhα = (−1)a sinhλr, (5.69) = (−1)a ( µ+ (log 2µ)2 2µ +O ( log µ µ )) . (5.70) cosh β = cosh ( api µ + i ( µ+ api − log 2µ µ +O ( (log µ)2 µ3 )) , (5.71) = (−1)a+1 ( log 2µ µ + iapi µ ) +O ( (log µ)2 µ3 ) . (5.72) cosh 2k∗m = cosh ( λr − iλi ) , (5.73) = − log 2µ+ iµ+O ( (log µ)2 µ ) . (5.74) cosh 2kn = cosh ( λr + 2api µ + i ( λi + 2api ) ) , (5.75) = − log 2µ+ iµ+O ( (log µ)2 µ ) . (5.76) 5.3 Short wavelength asymptotics 163 Using these expressions, the energy inner product Emn becomes Emn = ( 4µ2 (log 2µ)3 − i 12apiµ2 (log 2µ)4 +O ( 1 log µ )) + ( 4 log 2µ +O ( 1 (log µ)2 )) + ( − 2 log 2µ +O ( 1 (log µ)2 )) + ( − 2 log 2µ +O ( 1 (log µ)2 )) . (5.77) Thus to leading order we have Emn = 4µ2 (log 2µ)3 − i 12apiµ2 (log 2µ)4 +O ( 1 log µ ) . (5.78) Re-introducting the normalized energy inner product 〈q˜n+a, q˜n〉 = Emn EmEn , (5.79) we see that 〈q˜n+a, q˜n〉 = 1− i 3api log 2µ +O ( (log µ)2 µ2 ) , a = O(1), n→∞ (5.80) where µ = (2n− 12)pi. Hence in the limit |k| → ∞, modes are non-orthogonal. 5.3.3 Interaction of long and short wavelength modes We now determine the effect on the energy inner product of the interaction between long and short wavelength modes. We keep one mode fixed, and O(1), and take the limit of the second wavenumber becoming large. In terms of kn and km this means kn = log 2µ 2 + (log 2µ)2 4µ2 +O ( log µ µ2 ) + i ( µ 2 − log 2µ 2µ +O ( (log µ)2 µ3 )) , km = O(1). (5.81) We start by treating cosh kn = cosh { log 2µ 2 + (log 2µ)2 4µ2 +O ( log µ µ2 ) + i ( µ 2 − log 2µ 2µ +O ( (log µ)2 µ3 ))} , (5.82) = √ 2µ 2 ( 1 √ 2 ) + i √ 2µ 2 ( − 1 √ 2 ) +O(µ− 1 2 ), (5.83) = √ µ 2 (1− i) +O(µ− 1 2 ). (5.84) 164 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Therefore we have coshα = cosh(kn + k ∗ m), (5.85) = cosh kn cosh k ∗ m + sinh kn sinh k ∗ m, (5.86) = 1 2 √ µ(1− i) cosh k∗m + 1 2 √ µ(1− i) sinh k∗m +O(µ − 12 ), (5.87) =⇒ coshα = ek ∗ m 2 (1− i) √ µ+O(µ− 1 2 ). (5.88) Similarly sinhα = ek ∗ m 2 (1− i) √ µ+O(µ− 1 2 ), (5.89) cosh β = e−k ∗ m 2 (1− i) √ µ+O(µ− 1 2 ). (5.90) Using these expressions on the energy inner product (5.10b) we have Emn = ( 16(e2k ∗ m + 1) µ2 +O ( 1 µ3 )) − ( 8i µ +O ( 1 µ2 )) + ( 8i µ +O ( log µ µ2 )) + ( 8 cosh 2k∗m µ2 +O ( 1 µ3 )) . (5.91) The O(µ−1) term cancels, meaning Emn = O ( log µ µ2 ) , n→∞. (5.92) Hence we have 〈q˜n, q˜m〉 = O ( (log µ) 5 2 µ3 ) , (5.93) where µ = (2n− 12)pi. 5.3.4 Short wavelength growth mechanism The results of §§5.3.2 and 5.3.3 can be summarized by the following two leading order expressions for the energy norm (i) 〈q˜n, q˜m〉 ∼ 0, m = O(1), n→∞ (5.94a) (ii) 〈q˜n+a, q˜n〉 ∼ 1− i 3api log 2µ , a = O(1), n→∞. µ = ( 2n− 1 2 ) pi (5.94b) From (i) we have that, in the limit of large wavenumber, short wavelength modes are orthogonal to large wavelength modes. (ii) implies that the energy inner product of two short wavelength modes is non-zero, and modes become increasingly parallel on the slow (log 2µ)−1 scale. We interpret this to mean that short wavelength modes are able to 5.3 Short wavelength asymptotics 165 transfer energy between one another, but they transfer little energy to large wavelength modes. Therefore we believe that the optimal growth calculations in §5.1 exploit this energy transfer by including increasingly short wavelength modes in the expansion for the optimal edge profile. We expect that as the number of modes in the expansion is increased, the (n + 1)th mode, say, becomes almost parallel to the nth mode, and its influence on the optimal edge profile is diminished since its effects are be spanned by the nth mode. Furthermore, since the short wavelength modes are highly damped, as wavelengths are shortened modes are less able to affect the domain far from the edge. This offers an explanation as to why the location of the maximum energy growth, xmax, decreases with mode truncation M , as seen in the right panel of Figure 5.3. Figure 5.14: On the left panel we plot E2n against n. Red dots represent numerically calculated values and the blue line is our leading order estimate. On the right panel we have the norm 〈q˜n+1, q˜n〉 against n. The blue and red lines are the real and imaginary parts of the norm, and the blue and red dashed lines are the real and imaginary leading order estimates respectively. 166 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.4 Spatial Poiseuille flow To investigate the role that the growth mechanism observed in the previous chapters plays in flows with convection, we address the problem of optimal growth in two-dimensional spatially developing Poiseuille flow of a viscous fluid. For the parabolic Poiseuille base flowU = U(y)ex = (1−y2)ex and velocity and pressure perturbations (u, p) = [u˜, v˜, p˜](x, y), the linearized Navier-Stokes equations for steady fluid motion in a semi-infinite channel are U ∂u˜ ∂x + v˜U ′ = − ∂p˜ ∂x + 1 Re ∇2u˜, (5.95a) U ∂v˜ ∂x = − ∂p˜ ∂y + 1 Re ∇2v˜, (5.95b) ∂u˜ ∂x + ∂v˜ ∂y = 0, (5.95c) where ′ denotes differentiation with respect to y and ∇2 is the Laplacian operator. The Reynolds number Re = Uchν is based on the centre velocity Uc, the channel half-height h and the kinematic viscosity of the fluid ν. Equations (5.95) are non-dimensionalized in terms of these quantities. We solve on the domain Ω as defined in §5.1 and impose the no-slip boundary conditions on the walls y = ±1. From this formulation, Stokes flow and the Papkovitch-Fadle operator can be recovered in the limit Re → 0. The linear spatial stability of the system is found by Fourier transforming the perturbations in the streamwise co-ordinate x such that [u˜, v˜, p˜](x, y) = [u, v, p](y)eiαx, (5.96) for α ∈ C. Equations (5.95) now become U iαu+ vU ′ = −iαp+ 1 Re (D2 − α2)u, (5.97a) U iαv = −p′ + 1 Re (D2 − α2)v, (5.97b) iαu+ v′ = 0, (5.97c) where D = ddy . For a given Re, we can determine eigenvalues αn ∈ C which govern the stability of the perturbations. The equations (5.95) are elliptic in space, thus solutions of (5.97) yield upstream and downstream travelling modes. We know from temporal stability theory that Poiseuille flow is linearly stable for Re < 5772 (Orszag (1971), Schmid & Henningson (2001)), so we can ignore upstream modes by neglecting eigenval- ues with Im{αn} < 0 provided Re is sufficiently small. In this way we can parabolize the system and restrict our attention to disturbances evolving downstream from the edge, Γ = {x = 0,−1 ≤ y ≤ 1}, of the domain. 5.4 Spatial Poiseuille flow 167 A matrix method is used to find the eigenmodes of (5.97). Introducing u¯ = αu and v¯ = αv, (5.98) we can recast the system into the linear eigenproblem Lq = αMq, (5.99) where the vector q = [u v p u¯ v¯]T and the operators L and M are L =         1 ReD 2 −U ′ 0 0 0 0 1ReD 2 −D 0 0 0 −D 0 0 0 0 0 0 1 0 0 0 0 0 1         , (5.100) M =         iU 0 i 1Re 0 0 iU 0 0 1Re i 0 0 0 0 1 0 0 0 0 0 1 0 0 0         . (5.101) We then discretize the system in y using Chebyshev polynomials, following the procedure detailed in Schmid & Henningson (2001), and solve the resulting generalized eigenvalue problem for α and q. 5.4.1 Spatial spectra In Figure 5.15 we show that the leading eigenvalue from our calculations at low Re shows good convergence to the leading Papkovitch-Fadle eigenvalue, which we interpret as some validation of the code for arbitrary Re. In Figure 5.16 we plot spatial spectra at various Re. We see that as Re→ 0 the spectra are qualitatively similar to the Papkovitch-Fadle operator eigenvalues. The emergence of stationary downstream propagating modes is investigated in Figure 5.17. We find that for 7 < Re < 8, a set of upstream propagating stationary modes become increasingly damped as Re is increased. The modes then ‘turn around’ in space, changing from highly damped upstream to high damped downstream modes. In so doing, the imaginary part of the wavenumber switches from negative to positive infinity, meaning that the stationary upstream modes emerge from infinity. 168 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.15: Comparison between the leading eigenvalue of our finite Re computations and the Papkovitch-Fadle spectrum. The dots correspond to calculations at the given Re, and the dashed line is the Papkovitch-Fadle mode. Figure 5.16: Spatial spectra at various Re. Downstream/upstream propagating modes are coloured blue/red. Top row: Re = 1, 5, 7. Bottom row: Re = 10, 100, 1000. 5.4 Spatial Poiseuille flow 169 Figure 5.17: Emergence of stationary downstream modes from infinity. Down- stream/upstream modes are coloured blue/red. A selection of stationary upstream modes are circled for Re = 7. As Re is increased, the upstream modes become more damped before they ‘turn around’ in space to become highly damped downstream modes. The stationary downstream modes appear in our Figure at Re = 7.5. They become less damped as Re is further increased. 170 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.4.2 Optimal growth and convergence We now turn our attention to optimal growth. The formulation is the same as in §5.2, but since we no longer have a general analytic solution the energy norm (5.16) is now computed via numerical integration. Having solved (5.97) for un and vn where Lqn = αnMqn and qn = [un vn pn u¯n v¯n] T, (5.102) the energy inner product Emn = 〈un,um〉 = 1 2 ∫ 1 −1 u∗mun + v ∗ mvn dy, (5.103) is numerically approximated using Clenshaw-Curtis quadrature (Clenshaw & Curtis (1960); Boyd (2001)), which exploits the choice of Chebyshev polynomials for discretization in y. We then compute the spatial energy gain G(x) as before. Each mode in our calculations has a wall-normal discretization of N = 601 Chebyshev polynomials, which ensures that we can include well-resolved highly damped modes. Figure 5.18 compares curves of the energy amplification G(x) at different Re. Each curve is computed for a perturbation expansion with M = 500 modes. We see that for Re = 1, 10, G(x) is qualitatively similar to the corresponding curves for the Stokes prob- lem shown in Figure 5.3. At larger Re however, the characteristic oscillations in G(x) are no longer present, since the leading eigenvalue has zero real part. Figure 5.18 also compares the energy growth regions for each flow. A notable feature in common with Stokes flow is that there is a sharp growth gradient near the edge x = 0, for every Re considered, suggesting that the unbounded numerical abscissa of the Papkovitch-Fadle operator may extend to flows with Re large. With respect to our trunction of M = 500, we define the quantities Gmax and xmax Gmax = max x G(x) = G(xmax). (5.104) However, these maxima may change subject to an increase in the modal truncation. 5.4 Spatial Poiseuille flow 171 Figure 5.18: Gain curves in x. Top row: Re = 1, 10. Bottom row: Re = 100, 1000. To assess the convergence of G(x) four Reynolds numbers, Re = 1, 10, 100, 1000, are chosen and in each case x is fixed at a value close to xmax. In Figure 5.19 we plot the results. Note that the increases in G are stepped in M . This is because the spatial spectra consist of even and odd modes which are orthogonal, as was the case for the Papkovitch- Fadle operator in §5.2, in equations (5.10). This means if the optimal perturbation is comprised of even modes, say, then odd modes cannot contribute to the energy growth. We compare the magnitude of the optimal perturbation coefficients |am| and filter out modes with markedly lower magnitudes (see Figure 5.20), which make no contribution to the optimal growth. On this reduced basis of modes we can consider G at fixed x to be a function of M , and we re-introduce δGδM from §5.2.1. δG δM can be used as a diagnostic to compare convergence at each Re. δGδM is qualitatively similar to the Stokes case for Re = 1, 10. The following scalings are computed and plotted in Figure 5.21 Re = 0 : δG δM ∝M−2, x = 0.1, (5.105a) Re = 1 : δG δM ∝M−1.9, x = 0.055, (5.105b) Re = 10 : δG δM ∝M−1.7, x = 0.02. (5.105c) For Re = 100, 1000, G grows rapidly in the range 1 . M . 100. For M & 100 G 172 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator increases on a much slower scale; this is shown in on the lower row of Figure 5.19. For these Reynolds numbers we plot δGδM in Figure 5.22. We calculate the convergence rates Re = 100 : δG δM ∝M−2, x = 1.6, Re = 1000 : δG δM ∝M−2.5, x = 3.5. (5.106) Figure 5.19: Convergence of G in modal truncation M . The streamwise coordinate is fixed at x = 0.055, 0.02, 1.6 and 3.5 for Re = 1, 10, 103 and 103 respectively. Figure 5.20: Absolute value of the optimal coefficients |am| for Re = 1, 10, 100, 1000. Modes on the lower order of magnitude are filtered out of the optimal eigenmode expan- sion to form a reduced basis. 5.4 Spatial Poiseuille flow 173 Figure 5.21: Convergence rates δGδM for Re = 0 (top), Re = 1 (bottom-left) and Re = 10 (bottom-right). The dashed curves are the numerically calculated exponents in equations (5.105) Figure 5.22: Convergence rates δGδM for Re = 100 (left) and Re = 1000 (right). The dashed curves are the numerically calculated exponents in equations (5.106) 174 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.4.3 Optimal disturbances In keeping with the trends observed in the previous section, the optimal disturbances at low Re share similar qualitative characteristics with the optimals in the Stokes problem. We plot the modal form of the optimals in Figure 5.23, and their form in the (x, y)- plane in Figure 5.24. For larger Re the optimals still possess oscillatory features at the corners of the domain, however it is less pronounced than for the low Re optimals. At xmax, the optimals at large Re have maximum streamwise velocity u away from the wall, suggesting they are extracting energy from the mean shear in a manner similar to the two-dimensional Orr mechanism in temporal flows. Figure 5.23: Wall normal and optimal profiles with velocities u and v in blue and red, respectively. The top row has the edge optimals, and the bottom row has the optimals near xmax for Re = 1, 10, 100, 1000. 5.4 Spatial Poiseuille flow 175 Figure 5.24: Streamwise velocity u for Re = 1 (top) and Re = 1000 (bottom) in the (x, y)-domain. 176 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator 5.5 Temporal plane Couette flow We contrast our findings in the Papkovitch-Fadle problem with an analysis of temporally evolving and streamwise independent plane Couette flow. For this flow we find that numerical computations of the transient perturbation energy density growth converge quickly with respect to the modal truncation of the calculation. An analytic solution is available which allows us to see what is happening to the modes. We consider the motion of an incompressible, viscous fluid constrained between two infi- nite plates which move in opposite directions at equal speeds. The no-slip boundary con- dition is applied on each of the walls. The governing equations are non-dimensionalized with the wall speed U˜ , the channel half-height h˜ and the kinematic viscosity of the fluid ν. We seek a solution to the linearized perturbation equations with base velocity U = yex. Furthermore we restrict the solution to streamwise independent and spanwise periodic perturbations such that [u, v, w, p](y, z, t) = [u, v, w, p](y, t)eiβz. (5.107) If we further assume that the velocity components have exponential dependence in time, ∝ e−σt, the linearized Navier-Stokes equations are −σu+ v = 1 Re (u′′ − β2u), (5.108a) −σv = −p′ + 1 Re (v′′ − β2v), (5.108b) −σw = −iβp+ 1 Re (w′′ − β2w), (5.108c) v′ + iβw = 0, (5.108d) where D denotes differentiation with respect to y and the Reynolds number is defined as Re = U˜ h˜ν . The system can be further reduced to v′′′′ + (σRe− 2β2)v′′ − β2(σRe+ β2)v = 0, (5.109a) u′′ + (σRe− β2)u = Re v, (5.109b) subject to v(±1) = v′(±1) = u(±1) = 0. Using a Green’s function approach we derive a general solution for this system with three classes of modes. Writing µ2 = σRe− β2 we have Class (a) : µ tanµ+ β tanh β = 0, (even modes), 5.5 Temporal plane Couette flow 177    u = Reµ2+β2 ( cosh β cosµy − cosµ cosh βy ) + Re coshβ2µ cosµ ( y cosµ sinµy − sinµ cosµy ) , v = cosh β cosµy − cosµ cosh βy, w = − iβ ( µ cosh β sinµy + β cosµ sinh βy ) . Class (b) : µ cotµ− β coth β = 0, (odd modes)    u = Reµ2+β2 ( sinh β sinµy − sinµ sinh βy ) + Re sinhβ2µ sinµ ( cosµ sinµy − y sinµ cosµy ) , v = sinh β sinµy − sinµ sinh βy, w = iβ ( µ sinh β cosµy − β sinµ cosh βy ) . Class (c) : µ = npi 2    u = sinµ(y + 1), v = 0, w = 0. (u-modes) Given β ∈ R, the transcendental dispersion relations of the even and odd modes can be solved via a Newton iteration algorithm for µ ∈ R. The flow is found to be stable (σ > 0) ∀β. We will list the modes in terms of increasing µ. For example, µan denotes the n th a-mode with σan = (µan) 2 + β2 Re . (5.110) Similarly, the corresponding eigenfunctions or eigenmodes are denoted uan. Note that without loss of generality we can assume µ > 0, since the value of the decay rate depends on the square of µ. Despite the asymptotic stability of the flow significant transient energy density growth has been shown to exist at moderate Reynolds numbers (Butler & Farrell (1992), Reddy & Henningson (1993)). This growth is a consequence of the nonnormality of the linearized Navier-Stokes operator and the non-orthogonal nature of the eigenmodes of the system with respect to the energy inner product. Writing the kinetic energy through the fluid domain as E = 1 2 ∫ 1 −1 u∗u+ v∗v + w∗w dy, (5.111) we introduce the modal energy inner product 〈uλ,uµ〉 = 1 2 ∫ 1 −1 u∗λuµ + v ∗ λvµ + w ∗ λwµ dy. (5.112) Hence if we have an arbitrary velocity disturbance expanded in terms of the eigenmodes 178 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator of the system u(y, z, t) = ∞∑ n=1 κnuneiβz−σnt = ∞∑ j=1 ajuaj e iβz−σjt + ∞∑ k=1 bkubke iβz−σkt + ∞∑ `=1 c`uc`e iβz−σ`t, (5.113) then the appropriate energy norm is ‖u(t)‖2 = 〈u,u〉. (5.114) In terms of the modal energy we find the following orthogonality relations 〈uaλ,u a µ〉 6= 0, (5.115a) 〈uaλ,u b µ〉 = 0, ∀µ, λ (5.115b) 〈uaλ,u c µ〉 6= 0, (5.115c) 〈ubλ,u b µ〉 6= 0, (5.115d) 〈ubλ,u c µ〉 6= 0, (5.115e) 〈ucλ,u c µ〉 = 0 for µ 6= λ. (5.115f) The energy gain G(t) = max u E(t) E0 = max u0 ‖u(t)‖2 ‖u0‖ , (5.116) can be calculated in the same manner as in §5.2, with exact expressions for the energy in- ner product. The gain, and corresponding optimal disturbance, can also be computed via collocation in y as in Schmid & Henningson (2001). Computations plotted in Figure 5.25, show that the optimal initial condition which reaches Gmax, where Gmax = maxtG(t), is comprised of class (a) and (c) modes. The value of Gmax is found to be well converged for a modal truncation N  104, as can be seen in Figure 5.26. Since orthogonal modes can make no contribution to the optimal growth, we assume that in addition to being highly damped, modes with large σn orthogonalize so that they make a negligible contribution to the maximum optimal growth. Accordingly, we assess how the relevant modal energy inner products behave in the asymptotic limit of large decay rate. We begin with exact expressions for the inner products of interest 〈uaλ,u c µ〉 = Reµ cosλ cosh β 2 (cos 2µ− 1) ( 1 (λ2 − µ2)2 + 1 λ2 + β2 ( 1 λ2 − µ2 + 1 µ2 + β2 )) . (5.117) 5.5 Temporal plane Couette flow 179 〈uaλ,u a µ〉 = Re2 2(µ2 + β2)(λ2 + β2) { cosh2 β ( sin(µ− λ) µ− λ + sin(µ+ λ) µ+ λ ) + cosµ cosλ ( 1 + sinh 2β 2β )} + Re2 cosh2 β 8µλ cosµ cosλ { sin2(µ− λ) (µ− λ)2 − sin2(µ+ λ) (µ+ λ)2 + 2 cosµ cosλ ( cos(µ− λ) (µ− λ)2 − cos(µ+ λ) (µ+ λ)2 + sin(µ+ λ) (µ+ λ)3 − sin(µ− λ) (µ− λ)3 )} + Re2 cosh β 4µ cosµ(λ2 + β2) { cosµ cosh β ( sin(µ+ λ) (µ+ λ)2 + sin(µ− λ) (µ− λ)2 − cos(µ+ λ) µ+ λ − cos(µ− λ) µ− λ ) − cosµ cosλ ( 2(β sinh β sinµ− µ cosh β cosµ− cosh β sinµ) µ2 + β2 ) − sinµ cosh β ( sin(µ+ λ) µ+ λ + sin(µ− λ) µ− λ )} + Re2 cosh β 4λ cosλ(µ2 + β2) { cosλ cosh β ( sin(µ+ λ) (µ+ λ)2 − sin(µ− λ) (µ− λ)2 − cos(µ+ λ) µ+ λ + cos(µ− λ) µ− λ ) − cosµ cosλ ( 2(β sinh β sinλ− λ cosh β cosλ− cosh β sinλ) λ2 + β2 ) − sinλ cosh β ( sin(µ+ λ) µ+ λ + sin(µ− λ) µ− λ )} , for µ 6= λ. (5.118) 〈uaµ,u a µ〉 = Re2 2(µ2 + β2)2 { cosh2 β ( 1 + sin 2µ 2µ ) + cosµ2 ( 1 + sinh 2β 2β )} + Re2 cosh2 β 8µ2 cos2 µ { cos2 µ ( 1 3 − cos 2µ 2µ2 + sin 2µ 4µ3 ) − sin2 2µ 4µ2 + sin2 µ } + Re2 cosh β 2µ cosµ(µ2 + β2) { cosµ cosh β ( sin 2µ 4µ2 − cos 2µ 2µ ) − sinµ cosh β ( 1 + sin 2µ 2µ ) − cos2 µ ( 2(µβ sinh β sinµ− µ2 cosh β cosµ+ β sinh β cosµ) µ(µ2 + β2) )} + sinh 2β cos2 µ 2β + cosh2 β 2 ( 1 + sin 2µ 2µ + µ2 β2 ( 1− sin 2µ 2µ )) + 2µ cosµ cosh β β(µ2 + β2) ( β cosh β sinµ− µ cosµ sinh β ) . (5.119) 〈ucµ,u c µ〉 = 1 2 . (5.120) Note that the terms appearing in 〈uaλ,u a µ〉 are all proportional to Re 2, whereas 〈uaµ,u a µ〉 contains terms which are not. This indicates that as the Stokes limit, Re → 0, is ap- proached, the inner products orthogonalize, in contrast with the case for spatially evolving 180 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator disturbances. Returning to the problem at hand, we pass to the high decay rate limit and keep Re and β fixed, which gives µan ∼ npi − ε, n→∞. ε = β tanh β npi . (5.121) Thus 〈uaµ,u a µ〉 = pi2 cosh2 β 2β2 n2 +O(1), n→∞. (5.122) This means that the normalization factor for the a-modes has the asymptotic form Naµ = pi cosh β √ 2β n+O(1), n→∞ (5.123) whilst for the c-modes the normalization factor is simply a constant, N c = √ 2. Hence we can introduce asymptotically normalized modes u˜aµ = uaµ Naµ and u˜cµ = ucµ N c , (5.124) such that 〈u˜aµ, u˜ a µ〉 ∼ 1 and 〈u˜ c µ, u˜ c µ〉 = 1. (5.125) Now we can look at how the modal interactions influence the energy inner product. For the 〈u˜aµ, u˜ c λ〉 inner product we choose three relevant large n limits to investigate (i) µan = npi +O(n −1), µc = λ = O(1), (5.126a) (ii) µa = λ = O(1), µcn = ( n+ 1 2 ) pi, (5.126b) (iii) µan = npi +O(n −1), µcn = ( n+ 1 2 ) pi. (5.126c) For limit (i) we have 〈u˜aµan , u˜ c λ〉 = (−1) n+1 ( λβRe pi3(λ2 + β2) ) n−3 +O(n−5), n→∞. (5.127) Taking limit (ii) we get 〈u˜aλ, u˜ c µcn 〉 = ( − √ 2(λ2 + β2)Re cosλ cosh β Naλpi 5 ) n−5 +O(n−6), n→∞ (5.128) 5.5 Temporal plane Couette flow 181 where Naλ = √ 〈uλ,uλ〉. Limit (iii) gives 〈u˜aµan , u˜ c µcn 〉 = (−1)n+1 ( βRe pi4 ) n−2 +O(n−3), n→∞. (5.129) Now considering the 〈u˜aµ, u˜ a λ〉 inner product, we investigate the following limits (iv) µn = npi +O(n −1), λ = O(1), n→∞, (5.130a) (v) µn = npi +O(n −1), µn−p = (n− p)pi +O((n− p) −1), (5.130b) for p ∈ Z and p = O(1), n→∞. We have dropped the superscripts for µa in the above expressions, since both modes are class a in this case. For limit (iv) we have 〈u˜µn , u˜λ〉 = An −3 +O(n−5), n→∞ (5.131) where A = βRe2 √ 2pi3(λ2 + β2)Naλ ( β sinh 2β sinλ 2λ −cosh2 β ( cosλ+ sinλ λ ) −cosλ ( 1+ sinh 2β 2β )) . (5.132) In limit (v) we have, to leading order 〈u˜µn , u˜µn−p〉 = (−1) pβ 2Re2 2p2pi6 n−4 +O(n−5), n→∞. (5.133) Hence it can be seen that in all five limits considered, the asymptotic energy inner products (i), (ii), (iii), (iv) and (v) decay for highly damped modes. Though this is not an exhaustive list of large n limits, we believe that the behaviour shown for the limits considered will extend to others. We posit that this is a key component in the good convergence of optimal disturbances in plane Couette flow. 182 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator Figure 5.25: On the left is the σ spectrum in the complex plane, with β = 2 and Re = 500. On the right we plot the absolute value of the optimal coefficients |κn| for the optimal with β = 2, Re = 500 and t = 100. •, ◦ and ∗ are the a, b and c modes respectively. Figure 5.26: We plot G (left) and δGδN (right), for β = 2, Re = 500 and t = 100, against a modal truncation N comprised of a and c modes. 5.6 Conclusion 183 5.6 Conclusion Our focus in this chapter has been the spatial optimal growth of an arbitrary, finite expansion of the Papkovitch-Fadle eigenfunctions under boundary conditions motivated by Stokes flow. By analysis of the energy gain G(x) as a function of the number of Papkovitch-Fadle eigenmodes included in the velocity expansion, we have shown that G(x) is not well converged for x < 0.1, and that a large number of modes, M = 104, are required to guarantee convergence for x > 0.1. We analyze the form of the optimal edge condition, and find that increasingly large wavenumber modes are selected by the optimization procedure which finds the optimal edge condition for a given downstream location. The large wavenumber modes make a prominent contribution to the optimal edge condition at its corners, inducing highly oscillatory velocity fluctuations, and ex- erting a large shear stress on the corner walls. In light of this, we solve a constrained optimization problem in which the corner stress is bounded at unity, and we find that though less energy amplification is possible, the optimal edge condition retains the use of large wavenumber modes in its spatial energy growth. To gain an insight into why the optimal edge velocity is reliant on large wavenumber modes, we undertake an asymptotic analysis of the energy inner product in the large wavenumber limit. We show that the inner product of two large wavenumber modes is O(1), but that the energy inner product interactions of large and low wavenumber modes are asymptotically small. This suggests that there exists a mechanism by which large wavenumber modes instigate energy growth by transferring energy between one another. That such a mechanism exists does not imply that large wavenumber modes definitively will contribute to the maximum energy amplification, rather, that they could make a contribution. Therefore, our asymptotic results are to be interpreted as some justifica- tion for the inclusion of large wavenumber modes in the optimal edge condition, but not an indication that they must be present. Moving away from Stokes flow, we assess the impact of the large wavenumber growth mechanism on convected flows, choosing two-dimensional spatially developing Poiseuille flow to examine. We find that qualitiative features of the Stokes optimal edge condition are carried over into the finite Re flow; highly oscillatory velocity fluctuations are found near the channel corners, and the energy gain G(x) converges slowly as a function of the modal resolution. Finally, we analyzed temporal optimal growth in plane streamwise-independent plane Couette flow, for which large optimal energy amplification is well-known and no con- 184 5. Nonnormality and optimal growth of the Papkovitch-Fadle operator vergence issues have been reported. We find an analytic solution to the system, which allows us to perform an analogue of the asymptotic analysis of the spatial energy inner product, on the interactions of high and low damped modes. In contrast to our spatial results, we find that the temporal energy inner product invariably decays in the all high damping limits considered, meaning that highly damped modes cannot contribute to non-modal energy growth. This difference highlights the dichotomy between the spatial and temporal frameworks in fluid dynamics. Chapter 6 Conclusions & further work By way of concluding this thesis we will re-cap and contextualize the conclusions from the four results chapters §§2, 3, 4 and 5. In addition, with an outlook on potential future work we will discuss the ways in which our results can be developed. Each of the flows considered throughout the thesis will be discussed in turn, in §§6.1, 6.2 and 6.3. 6.1 Annular Poiseuille-Couette flow In §2 we contribute to existing stability analyses in annular Poiseuille-Couette flow (APCF) with a study of the non-modal linear energy growth. Beginning with a lin- ear stability analysis, we demonstrate that the Rayleigh stability criterion can be used to show that APCF is stable to inviscid perturbations for inner cylinder velocity V < V ∗. However, the Rayleigh stability criterion is only a necessary condition for instability, and we find no indication that the flow loses stability to inviscid waves for V ≥ V ∗. Vis- cous waves are known to de-stabilize APCF (Sadeghi & Higgins (1991), Gittler (1993), Walton (2003, 2004, 2005)), and by maximizing over axial and azimuthal wavenumbers we produce neutral curves in (Re, V )-space, for suitably defined Reynolds number Re. Though the neutral curves are not monotonic in (Re, V ), the trend is for non-zero V to stabilize the laminar base flow. There appear to be cut-off velocities for positive and negative V , beyond which all perturbations are linearly stable, and our calculations are in agreement with the cut-off velocity for radius ratio η = 0.5 and V > 0 derived by Sadeghi & Higgins (1991). Our non-modal growth analysis focuses on subcritical (Re, V ) space, and we find that our conclusion differs from linear stability results: non-zero V has the tendency to de-stabilize the laminar profile, in the sense that the maximum energy gain G increases with |V |, provided |V | is sufficiently large. This is demonstrated in a series of Figures in which contours of G are plotted in (Re, V ) space for a range of η. The level of maximum growth is found to be of a similar order of magnitude regardless of the radius ratio η. Non-axisymmetric modes are found to be amplified by the well-known 185 186 6. Conclusions & further work lift-up mechanism, while axisymmetric modes are amplified via the Orr mechanism. The presence of both energy growth mechanisms is a strong indication that a subcritical route to transition may be possible in a similar manner to that which is currently understood in the canonical subcritical shear flows: plane Poiseuille flow, plane Couette flow and Hagen-Poiseuille flow. Recent studies such as Pringle et al. (2012), Rabin et al. (2012) and Duguet et al. (2013) indicate, via optimization over the fully nonlinear Navier-Stokes equations, that the minimal amplitude perturbations which instigate transition amplify themselves via combinations of the Orr mechanism and the lift-up mechanism. Addi- tionally, we find that the experimentally observed (Re, V ) transition values from Shands et al. (1980) lie at similar values of G, which we find to be suggestive of a link between the lift-up mechanism and transition to turbulence. Further study in APCF might address the existence of exact coherent structures in sub- critical (Re, V ) space. Wong & Walton (2012) have already found axisymmetric solutions which emerge in a bifurcation from linearly unstable axisymmetric perturbations, and can be continued into the linearly stable regime. This procedure could be extended to non-axisymmetric solutions, whose existence would provide further evidence for the link between exact coherent structures and transition to turbulence in subcritical shear flows. The relationship between the structures and the asymptotic nonlinear solutions of Wal- ton (2003, 2005) could be investigated. An asymptotic analysis in search of linear cut-off velocities for V < 0, following the approach of Sadeghi & Higgins (1991) who treated V > 0, might prove interesting. Such a study could provide asymptotic confirmation that large magnitude inner cylinder velocities stabilize linear modes for all Re, regardless of whether the inner cylinder motion opposes or aligns with the axial pressure gradient. 6.2 Rotating plane Couette flow In §3 we study the stability and bifurcation sequences of nonlinear solutions in supercrit- ical rotating plane Couette flow (RPCF). We present bifurcations of steady secondary solutions in rotation number and cross-flow energy space, (Ro,Ecf ). In addition to well- known Taylor vortex solutions, we discuss steady oblique vortex flows (OVF), finding upper and lower solutions which bifurcate from unstable modes of the primary flow. The steady OVF are related to the travelling wave spiral vortex flows discussed for Taylor- Couette flow (TCF) by, for example, Hoffmann et al. (2009), Altmeyer et al. (2010) and Deguchi & Altmeyer (2013). Our analysis contributes to the understanding of these oblique flows, and elucidates the differences in bifurcation behaviour between the oblique and streamwise independent flow states. 6.2 Rotating plane Couette flow 187 An extensive stability analysis of Taylor vortex flow (TVF) in §3.6 yields an insight into the bifurcation sequences which develop as Taylor vortices lose stability. Stability calculations across a wide range of TVF allows the characterization of the instabilities affecting TVF and the tertiary flows into which they develop. Alongside the tertiary flows which have already been discovered in RPCF, such as wavy vortex flow (WVF) and twisted vortex flows (TWI and wTWI), in §3.7 we find an oscillating vortex struc- ture (oWVF) and a streamwise indepedent short spanwise wavelength structure which we have named skewed vortex flow (SVF). Bifurcation analyses in Ro and streamwise wavenumber αs reveal that for small αs wTWI emerges in a saddle-node bifurcation from TVF. A solution branch of wTWI is found to exist which connects TVF states with βs = 1.5 and βs = 3, suggesting a more intricate relationship between the two flows than has previously been appreciated. A bifurcation scenario is presented for oWVF where no stable solution exists and flow moves between unstable states in a quasi-chaotic orbit. There is perhaps some novelty that such an orbit can exist in a transitional flow regime, and a related process could form the basis for the turbulent dynamics found at larger Re. In a bifurcation analysis of SVF, we find that it connects the first and second TVF which bifurcate from the primary flow. The stability properties of representative cases for all of these tertiary flows have been analyzed, suggesting the existence of quaternary structures to be investigated in future research. The existence of such a range of tertiary and higher order nonlinear structures could, at larger Reynolds numbers, be important in the understanding of fully turbulent flow. The relevence of tertiary states to transitional RPCF is assessed in §3.8. Through nu- merical simulations on a domain which supports a range of wavenumbers and through a model which simplifies the transition process based on bifurcation sequences, we analyze the states observed in the experiments of Tsukahara et al. (2010) and Suryadi et al. (2013). We find good agreement between model, simulation and experiment in most cases. However, a notable divergence between model and simulation occurs for Ro = 0.9, where the model predicts a final flowfield consisting of wTWI and the simulation finishes in a final flowfield with a localized twist vortex amongst streamwise independent TVF. Despite this deviation, the model and simulation flowfields are good approximations for two experimentally observed flowfields at Ro = 0.9. Further studies could include both a Floquet stability analysis of secondary oblique vor- tices and a study of the bifurcation scenarios which follow when it loses stability. We have some preliminary results in this direction, having found a tertiary oblique wave structure reminiscent of wavy spiral vortex flow of TCF discussed in Hoffmann et al. (2009) and Deguchi & Altmeyer (2013), shown in Figure 6.1. We name this structure wavy oblique 188 6. Conclusions & further work vortex flow wOVF. Figure 6.1: Two wOVF flowfields for (αs, βs, Re,Ro) = (0.2, 2, 100, 0.1) which bifurcate from steady OVF at the same parameters. It would be interesting to explore the connection between wOVF and the turbulent stripe phenomenon in non-rotating plane Couette flow (PCF). Turbulent stripes are a curious feature known to the fluid dynamics community since they were reported in experiments by Coles (1965), whose images caught the attention of Richard Feynman, and prompted him to coin the phrase “barber-pole turbulence” in description of the stripes (Feynman (1964)). Turbulent stripes are found in experiments and simulations of PCF (see, for example, the recent studies of Prigent et al. (2002), Duguet et al. (2010), Duguet & Schlatter (2013)), however, the origin of the oblique structure remains an open question. If wOVF or a similar state can be continued from RPCF to PCF, it could serve as a co- herent structure from which obliqueness of the stripes could originate, and about which the turbulent dynamics could be organized. A non-modal growth analysis across subcritical RPCF has been conducted in §4. By optimization over streamwise and spanwise wavenumbers α and β, we compute the maximum linear energy amplification, G, available to perturbations across subcritical (Ro,Re) space. We find that the effect of strong cyclonic (Ro < 0) or anti-cyclonic (Ro > 0) rotation is to hinder the lift-up and anti lift-up mechanisms, such that optimal perturbations rely on the Orr mechanism for energy growth. The restriction of the anti lift-up mechanism in anti-cyclonic RPCF has been noted before by Yecko (2004) and Rincon et al. (2007), whose results were interpreted in the context of transition in astro- physical accretion disks. Our contribution is therefore to show that strong cyclonic RPCF behaves comparably to strong anti-cyclonic RPCF, with respect to optimal disturbances. We posit that the absence of the lift-up mechanism in strong cyclonic RPCF will alter the transition scenario as it is currently understood in PCF. The lift-up mechanism is 6.3 The Papkovitch-Fadle operator and spatial non-modal growth 189 thought to be relevant in generating the streak component of the self-sustaining process (Waleffe (1995, 1997, 1998)) and vortex-wave interaction theory (Hall & Smith (1991); Hall & Sherwin (2010), Deguchi et al. (2013), Blackman et al. (2013)), which have been proposed to explain the existence of exact coherent structures in the canonical subcritical shear flows. Furthermore a lift-up mechanism growth phase has been highlighted in the nonlinear optimals found by Pringle et al. (2012), Rabin et al. (2012) and Duguet et al. (2013), for Hagen-Poiseuille flow and PCF. Hence, the absence of the lift-up mechanism would, in theory, prevent these transition scenarios from developing and a new under- standing would be required if indeed turbulence was shown to develop in RPCF with strong rotation. Further study could be made into transition in subcritical RPCF: in the three regimes of large cyclonic rotation, large anti-cyclonic rotation and Ro = 1. In the latter case, the anti lift-up mechanism and the Orr mechanism co-exist and Lesur & Longaretti (2005) found that transition occured in numerical simulations. A novel route to transition may then be found for this subcritical shear flow. A search for exact coherent structures, renewing the efforts of Rincon et al. (2007), may prove illuminating. Bifurcations from the secondary and tertiary structures of §3 in supercritical RPCF would be the ideal starting point for such a search. Additionally, a study similar to Reddy et al. (1998), with an analysis of the secondary stability of transiently amplified structures, could be carried out on the optimal disturbances of subcritical RPCF. 6.3 The Papkovitch-Fadle operator and spatial non- modal growth In §5 we have studied the non-modal stability properties of the Papkovitch-Fadle oper- ator with application to Stokes flow. We show that convergence of the maximum linear spatial energy amplification in modal truncation is slow, due to reliance of the optimal edge condition on highly damped modes. We demonstrate that the corner shear-stress of the optimal edge condition scales linearly with the number of modes used to expand the optimal disturbance. However, we solve a constrained optimization problem in which the corner stress is fixed at unity, and we show that the constrained optimal edge condi- tion has similar qualitative features to the unconstrained case. An asymptotic analysis of the energy inner product interactions of the Papkovitch-Fadle eigenfunctions reveals that short wavelength modes can transfer energy between one another, but not to a long wavelength mode. This result offers an explanation as to why optimization for the maxmimum energy amplification selects increasingly short wavelength modes. The impli- 190 6. Conclusions & further work cations of our results on convected flows are investigated in an analysis of optimal growth in two-dimensional spatial Poiseuille flow, for which the Papkovitch-Fadle operator can be considered the Re → 0 limit. We find that the optimals in spatial Poiseuille flow share qualitative characteristics with the optimal edge conditions for the Papkovitch- Fadle operator, suggesting that the short wavelength energy growth mechanism is not unique to Stokes flow. Finally, our results are contrasted with a treatment of temporal plane Couette flow, where it is shown that no analogus highly-damped modal coupling does not occurs. A continuation of our ideas on short wavelength energy transfer in spatially develop- ing flows may serve to complete our understanding of the mechanisms by which fluid structures are energetically amplified. Henningson & Schmid (1994) pointed out that the kinetic energy norm need not behave in the same way with respect to spatially de- veloping disturbances as it does for temporal disturbances. A new approach could be to construct a novel disturbance measure which is not affected by a short wavelength energy transfer, based on the bi-orthogonality properties of the Papkovitch-Fadle eigen- functions, though the physical interpretation of such a norm would need to be considered. Temporal growth mechanisms are well understood, and have been discussed at length in this thesis. A richer understanding of spatial growth mechanisms, in any norm, would serve to complement our existing knowledge on the important issue of spatio-temporal devlopment in fluid flow. Bibliography Alfredsson, P. H. & Tillmark, N. 2005 Instability, transition and turbulence in plane Couette flow with system rotation. In IUTAM Symposium on Laminar Turbulent Transition and Finite Amplitude Solutions (ed. T. Mullin & R. R. Kerswell), Fluid Mechanics and Its Applications , vol. 77. Springer. Altmeyer, S., Hoffmann, Ch., Heise, M., Abshagen, J., Pinter, A., Lu¨cke, M. & Pfister, G. 2010 End wall effects on teh transition between Taylor vortices and spiral vortices. Phys. Rev. E 81. Andereck, C. D., Liu, S. S. & Swinney, H. L. 1986 Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech. 164, 155–183. Andersson, P., Berggern, M. & Henningson, D. S. 1999 Optimal disturbances and bypass transition in boundary layers. Phys. Fluids 11 (1). Antkowiak, A. & Brancher, P. 2007 On vortex rings around vortices: an optimal mechanism. J. Fluid Mech. 578, 295–304. Antonijoan, J. & Sa´nchez, J. 2000 Transitions from Taylor vortex flow in a co- rotating Taylor-Couette system. Phys. Fluids 12. Bergstro¨m, L. 1993 Optimal growth of small disturbances in pipe Poiseuille flow. Phys. Fluids A 5. Blackman, H. M., Hall, P. & Sherwin, S. J. 2013 Lower branch equilibria in Couette flow: the emergence of canonical states for arbitrary shear flows. J. Fluid Mech. 726, 35–97. Blum, E. K. & Chang, A. F. 1978 A numerical method for the solution of the double eigenvalue problem. IMA J. Appl. Math. 22, 29–42. Boyd, J. P. 2001 Chebyshev and Fourier Spectral Methods . Dover, New York. Butler, K. M. & Farrell, B. F. 1992 Three-dimensional optimal perturbations in viscous shear flow. Phys. Fluids 4. 191 Chapman, S. J. 2002 Subcritical transition in shear flows. J. Fluid Mech. 451, 35–97. Chevalier, M., Schlatter, P., Lundbladh, A. & Henningson, D. S. 2007 A pseudo-spectral solver for incompressible boundary layer flows. Tech. Rep. TRITA- MEK 2007:07. KTH Mechanics, Stockholm, Sweden. Clenshaw, C. W. & Curtis, A. R. 1960 A method for numerical integration on an automatic computer. Numer. Math. 2 (197). Coles, D. 1965 Transition in circular Couette flow. J. Fluid Mech. 21, 385–425. Corbett, P. & Bottaro, A. 2001 Optimal linear growth in swept boundary layers. J. Fluid Mech. 435, 1–23. Cowley, S. J. & Smith, F. T. 1985 On the stability of Poiseuille-Couette flow: a bifurcation from infinity. J. Fluid Mech. 156, 83–100. Czarny, O., Serre, E., Bontoux, P. & Lueptow, R. M. 2003 Interaction between Ekman pumping and the centrifugal instability in Taylor-Couette flow. Phys. Fluids 15. Davey, A., DiPrima, R. C. & Stuart, J. T. 1968 On the instability of Taylor vortices. J. Fluid Mech. 31, 17–52. Deguchi, K. & Altmeyer, S. 2013 Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow. Phs. Rev. E 87. Deguchi, K., Hall, P. & Walton, A. 2013 The emergence of localized vortex-wave interaction states in plane Couette flow. J. Fluid Mech. 721, 58–85. Duguet, Y., Monokrousos, A., Brandt, L. & Henningson, D. S. 2013 Minimal transition thresholds in plane Couette flow. Phys. Fluids 25. Duguet, Y. & Schlatter, P. 2013 Oblique laminar-turbulent interfaces in plane shear flows. Phys. Rev. Lett. 110. Duguet, Y., Schlatter, P. & Henningson, D. S. 2010 Formation of turbulent patterns near the onset of transition in plane Couette flow. J. Fluid Mech. 650, 119– 129. Eagles, P. M. 1971 On stability of Taylor vortices by fifth-order amplitude expansions. J. Fluid Mech. 49, 529–550. Eckhaus, W. 1965 Studies in Nonlinear Stability Theory . Springer. Ellingsen, T. & Palm, E. 1974 Stability of linear flow. Phys. Fluids 18. 192 Fadle, J. 1941 Die selbstpannungs-eigenwertfunktionen der quadratischen scheibe. In- gen. Arch. 11, 125–149. Feynman, R. P. 1964 Lecture Notes in Physics, Vol. 2 . Addison-Wesley. Frei, Ch., Lu¨scher, P. & Wintermantel, E. 2000 Thread-annular flow in vertical pipes. J. Fluid Mech. 410, 185–210. Gibson, J. F. 2008 Channelflow: a spectral Navier-Stokes simulator in C++. Tech. Rep.. Georgia Institute of Technology. Gibson, J. F., Halcrow, J. & Cvitanovic´, P. 2009 Equilibrium and travelling-wave solutions of plane Couette flow. J. Fluid Mech. 683, 243–266. Gittler, Ph. 1993 Stability of axial Poiseuille-Couette flow between concentric cylin- ders. Acta Mech. 101, 1–13. Gregory, R. D. 1979 Green’s functions, bi-linear forms and the completeness of the eigenfunctions for the elastostatic strip and wedge. J. Elasticity 9, 283–309. Gregory, R. D. 1980a The semi-infinite strip x ≥ 0,−1 ≤ y ≤ 1; completeness of the eigenfunctions when φxx(0, y), φyy(0, y) are prescribed. J. Elasticity 10, 57–80. Gregory, R. D. 1980b The traction boundary value problem for the elastostatic semi- infinite strip; the existence of solutions and completeness of the Papkovitch-Fadle eigenfunctions. J. Elasticity 10, 295–327. Gue´gan, A., Schmid, P. J. & Huerre, P. 2006 Optimal energy growth and optimal control in swept Hiemenz flow. J. Fluid Mech. 566, 11–45. Guo, Y. & Finlay, W. H. 1991 Splitting, merging and wavelength selection of vortices in curved and/or rotating channel flow due to Eckhaus instability. J. Fluid Mech. 228, 661–691. Gustavsson, L. H. 1991 Energy growth of three-dimensional disturbances in plane Poiseuille flow. J. Fluid Mech. 224, 241–260. Hall, P. & Sherwin, S. J. 2010 Streamwise vortices in shear flows: harbingers of transition and the skeleton of coherent structures. J. Fluid Mech. 661, 178–205. Hall, P. & Smith, F. 1991 On strongly nonlinear vortex/wave interactions in boundary-layer transition. J. Fluid Mech. 227, 641–666. Heaton, C. J. & Peake, N. 2007 Transient growth in vortices with axial flow. J. Fluid Mech. 587, 271–301. 193 Hegseth, J. J., Baxter, G. W. & Andereck, C. D. 1996 Bifurcations from Taylor vortices between corotating concentric cylinders. Phys. Rev. E 53, 507–521. Heise, M., Hoffmann, Ch., Will, Ch., Altmeyer, S., Abshagen, J. & Pfister, G. 2013 Co-rotating Taylor-Couette flow enclosed by stationary disks. J. Fluid Mech. 716, 507–521. Henningson, D. S. & Schmid, P. J. 1994 A note on measures of disturbance size for spatially evolving flows. Phys. Fluids 6, 2862–2864. Hiwatashi, K., Alfredsson, P. H., Tillmark, N. & Nagata, M. 2007 Exper- imental observations of instabilities in rotating plane couette flow. Phys. Fluids 19, 48–103. Hoffmann, Ch., Altmeyer, S., Pinter, A. & Lu¨cke, M. 2009 Transitions between Taylor vortices and Spirals via wavy Taylor vortices and wavy Spirals. New J. Phys. 11. Hollerbach, R. & Fournier, A. 2004 End-effects in a rapidly rotating cylindrical Taylor-Couette flow. MHD Couette flows: Experiments and Models 776, 114–121. Itano, T. & Generalis, S. C. 2009 Hairpin vortex solution in planar Couette flow: a tapestry of knotted vortices. Phys. Rev. Lett. 102, 114–501. Ji, X. 1996 On 2D bisection method for double eigenvalue problems. J. Comput. Phys. 126, 92–96. Joseph, D. D. 1977 The convergence of biorthogonal series for biharmonic and Stokes flow edge problems. part I. SIAM J. Appl. Math. 33, 337–347. Joseph, D. D. & Sturges, L. 1978 The convergence of biorthogonal series for bihar- monic and Stokes flow edge problems. part II. SIAM J. Appl. Math. 34, 7–26. Kawahara, G. & Kida, S. 2001 Periodic motion embedded in plane Couette turbu- lence: regeneration cycle and burst. J. Fluid Mech. 449, 291–300. Koschmieder, E. L. 1993 Be´nard cells and Taylor vortices . Cambridge University Press. Kreilos, T. & Eckhardt, B. 2012 Periodic orbits near the onset of chaos in plane Couette flow. Chaos 22. Landahl, M. T. 1980 A note on the algebraic instability of inviscid shear flows. J. Fluid Mech. 98, 243–251. 194 Lesur, G. & Longaretti, P.-Y. 2005 On the relevance of hydrodynamic turbulence to accretion disk transport. Astron. Astrophys. 444, 25–44. Lezius, D. K. & Johnston, J. P. 1976 Roll-cell instabilities in rotating laminar and turbulent channel flows. J. Fluid Mech. 77, 153–175. Liu, R. & Liu, Q. S. 2012 Non-modal stability in sliding Couette flow. J. Fluid Mech. 484, 505–544. Longaretti, P.-Y. 2002 On the phenomenology of hydrodynamic shear turbulence. Astrophys. J. 576, 587–598. Luchini, P. 2000 Reynolds-number-independent instability of a boundary layer over a flat surface: optimal perturbations. J. Fluid Mech. 484, 289–309. Moffatt, H. K. 1963 Viscous and resisitive eddies near a sharp corner. J. Fluid Mech. 18. Nagata, M. 1986 Bifurcations in Couette flow between almost corotating cylinders. J. Fluid Mech. 169, 229–250. Nagata, M. 1988 On wavy instabilities of the Taylor-vortex flow between corotating cylinders. J. Fluid Mech. 188, 585–598. Nagata, M. 1990 Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity. J. Fluid Mech. 217, 519–527. Nagata, M. 1998 Tertiary solutions and their stability in rotating plane Couette flow. J. Fluid Mech. 358, 357–378. Nagata, M. 2013 A note on the mirror-symmetric coherent structure in plane Couette flow. J. Fluid Mech. 727, R1. Nagata, M. & Kawahara, G. 2004 Three-dimensional periodic solutions in rotating/non-rotating plane Couette flow. In Advances in Turbulence X, proceedings of the Tenth European Turbulence Conference (ed. H. I. Andersson & P. A˚. Krogstad), CIMNE , vol. 77. Narayan, R., Goldreich, P. & Goodman, J. 1987 Physics of modes in a differen- tially rotating system - analysis of the shearing sheet. R. Astron. Soc. 228, 1–41. Navier, C. L. M. H. 1822 Memoire sur les lois du mouvement des fluides. Mem. Acad. Sci. Inst. France 6, 389–440. 195 Orr, W. M. F. 1907 The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. Part I: A perfect liquid. Part II: A viscous liquid. Proc. R. Irish Acad. A 27, 9–138. Orszag, S. A. 1971 Accurate solutions of the Orr-Sommerfeld stability equation. J. Fluid Mech. 50, 689–. Papkovitch, P. F. 1940 U¨ber eine form der lo¨sung des byharmonischen problems fu¨r das rechteck. Russian Acad. Sci. Dokl. Math. 27, 334–338. Prigent, A., Gre´goire, G., Chate´, H., Dauchot, O. & van Saarloos, W. 2002 Large-scale finite-wavelength modulation within turbulence shear flows. Phys. Rev. Lett. 89. Pringle, C. T., Willis, A. P. & Kerswell, R. R. 2012 Minimal seeds for shear flow turbulence: using nonlinear transient growth to touch the edge of chaos. J. Fluid Mech. 702, 415–443. Proudman, J. 1916 On the motion of solids in a liquid possessing vorticity. Proc. R. Soc. Lond. 92, 408–424. Rabin, S. M. E., Caulfield, C. P. & Kerswell, R. R. 2012 Triggering turbulence efficiently in plane Couette flow. J. Fluid Mech. 712, 244–272. Rayleigh, L. 1880 On the stability, or instability, of certain fluid motions. Proc. Lond. Math. Soc. 174, 57–70. Reddy, S. C. & Henningson, D. S. 1993 Energy growth in viscous channel flows. J. Fluid Mech. 252, 209–283. Reddy, S. C., Schmid, P. J. & Baggett, J. S. Henningson, D. S. 1998 Stability of streamwise streaks and transition thresholds in plane channel flows. J. Fluid Mech. 365, 269–303. Reshotko, E. & Tumin, A. 2001 Spatial theory of optimal disturbances in a circular pipe flow. Phys. Fluids 13, 991–996. Reynolds, O. 1883 An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels. Phil. Trans. R. Soc. Lond. 174, 935–982. Rincon, F., Ogilvie, G. I. & Cossu, C. 2007 On self-sustaining processes in Rayleigh- stable rotating plane Couette flows and subcritical transition to turbulence in accretion disks. Astron. Astrophys. 463, 817–832. 196 Romanov, V. A. 1973 Stability of plane-parallel Couette flow. Funct. Anal. Appl. 7, 137–146. Saad, Y. & Schultz, M. H. 1986 GMRES - A generalized minimal residual algorithm for solving nonsymmetric linear-systems. SIAM J. Sci. Stat. Comp. 7, 856–869. Sadeghi, V. M. & Higgins, B. C. 1991 Stability of sliding Couette-Poiseuille flow in an annulus subject to axisymmetric and asymmetric disturbances. Phys. Fluids A 3, 2092–2104. Sa´nchez, J., Net, M., Garc´ıa-Archilla, B. & Simo´, C. 2004 Newton-Krylov continuation of periodic orbits for Navier-Stokes flows. J. Comput. Phys 201, 13–33. Schmid, P. J. & Henningson, D. S. 1994 Optimal energy density growth in Hagen- Poiseuille flow. J. Fluid Mech. 277, 197–225. Schmid, P. J. & Henningson, D. S. 2001 Stability and Transition in Shear Flows . Springer. Shands, J., Alfredsson, H. & Lindgren, E. R. 1980 Annular pipe flow subject to axial motion of the inner boundary. Phys. Fluids 23, 2144–2145. Shartman, E., Hantao, J., Burin, M. J. & Goodman, J. 2012 Stability of quasi- Keplerian shear flow in a laboratory experiment. Astron. Astrophys. 543, 1–13. Smith, R. C. T. 1952 The bending of a semi-infinite strip. Aust. J. Scient. Res. 5, 227. Spence, D. A. 1982 A class of biharmonic end-strip problems arising in elasticity and Stokes flow. IMA J. Appl. Math. 30, 107–139. Stokes, G. G. 1845 On the theories of internal friction of fluids in motion, and of the equilibrium and motion of elastic solids. Trans. Cam. Phil. Soc. 8, 279–319. Suryadi, A., Tillmark, N. & Alfredsson, P. H. 2013 Rotating plane Couette flow at high rotation number. In 65th Annual Meeting of the APS Division of Fluid Dynamics , , vol. 57. Taylor, G. I. 1923 Stability of a viscous liquid contained between two rotating cylin- ders. Phil. Trans. R. Soc. Lond. 223, 289–343. Tempelmann, D., Hanifi, A. & Henningson, D. S. 2010 Spatial optimal growth in three dimensional boundary layers. J. Fluid Mech. 646, 5–37. Tevzadze, A. G., Chagelishvili, G. D., Zahn, J.-P., Chanishvili, R. G. & Lominadze, J. G. 2003 Hydrodynamic shear turbulence in stratified keplerian disks: 197 transient growth of small-scale 3D vortex mode perturbations. Astron. Astrophys. 407, 779–786. Trefethen, L. N. 1997 Pseudospectra of linear operators. SIAM Review 39, 383–406. Trefethen, L. N. & Embree, M. 2005 Spectra and Pseudospectra. Princeton Uni- versity Press. Trefethen, L. N., Trefethen, A. E., Reddy, S. C. & Driscoll, T. A. 1993 Hydrodynamic stability without eigenvalues. Science 261, 578–584. Tsukahara, T., Tillmark, N. & Alfredsson, P. H. 2010 Flow regimes in plane Couette flow with system rotation. J. Fluid Mech. 648, 5–33. Viswanath, D. 2007 Recurrent motions within plane Couette turbulence. J. Fluid Mech. 580, 339–358. Waleffe, F. 1995 Transition in shear flows. Nonlinear normality versus non-normal linearity. Phys. Fluids 7, 3060–3066. Waleffe, F. 1997 On a self-sustaining process in shear flows. Phys. Fluids 9, 883–900. Waleffe, F. 1998 Three-dimensional coherent states in plane shear flows. Phys. Rev. Lett. 81, 4140–4143. Waleffe, F. 2003 Homotopy of exact coherent structures in plane shear flows. Phys. Fluids 15, 1517–1534. Walton, A. G. 2003 The nonlinear stability of thread-annular flow at high reynolds number. J. Fluid Mech. 477, 227–257. Walton, A. G. 2004 Stability of circular Poiseuille-Couette flow to axisymmetric dis- turbances. J. Fluid Mech. 500, 169–210. Walton, A. G. 2005 The linear and nonlinear stability of thread-annular flow. Phil. Trans. R. Soc. A 363, 1223–1233. Weisshaar, E., Busse, F. H. & Nagata, M. 1991 Twist vortices and their instabil- ities in the Taylor-Couette system. J. Fluid Mech. 226, 549–564. Wong, A. H. & Walton, A. H. 2012 Axisymmetric travelling waves in annular Couette-Poiseuille flow. Q. J. Appl. Math. 65, 293–311. Yecko, P. A. 2004 Accretion disk instability revisited: transient dynamics of rotating shear flow. Astron. Astrophys. 425, 383–395. 198