Two dimensional hydrodynamic instabilities in shear flows by Anirban Guha M.A.Sc., Mechanical Engineering, University of Windsor, ON, Canada, 2008 B.E., Mechanical Engineering, Jadavpur University, India, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June 2013 c© Anirban Guha 2013 Abstract Hydrodynamic instabilities occurring in two dimensional shear flows have been investigated. First, the process of resonant interaction between two progressive interfacial waves is studied. Such interaction produces exponentially grow- ing instabilities in idealized, homogeneous or density stratified, inviscid shear layers. It is shown that two oppositely propagating interfacial waves, hav- ing arbitrary initial amplitudes and phases, eventually phase-lock, provided they satisfy a particular condition. Three types of shear instabilities - Kelvin Helmholtz, Holmboe and Taylor have been studied. The above-mentioned condition provides a range of unstable wavenumbers for each instability type, and this range matches the predictions of the canonical normal-mode based linear stability theory. The non-linear evolution of Kelvin-Helmholtz (KH) instability has been studied. The commonly known manifestation of KH is in the form of spiral billows. However, KH evolving from a piecewise linear shear layer is remark- ably different; it is characterized by elliptical vortices of constant vorticity connected via thin braids. Using direct numerical simulation and contour dynamics, it is shown that the interaction between two counter-propagating vorticity waves is solely responsible for this KH formation. The oscillation of the vorticity wave amplitude, the rotation and nutation of the elliptical vortex, and straining of the braids have been investigated. Finally, the linear stability of plane Couette-Poiseuille flow in the presence of a cross-flow is studied. The base flow is characterized by the cross flow Reynolds number, Reinj and the dimensionless wall velocity, k. Corresponding to each dimensionless wall velocity, k ∈ [0, 1], two ranges of Reinj exist where ii Abstract unconditional stability is observed. In the lower range of Reinj , for modest k we have a stabilization of long wavelengths leading to a cut-off Reinj. As Reinj is increased, we see first destabilization and then stabilization at very large Reinj. Analysis of the eigenspectrum suggests the cause of instability is due to resonant interactions of Tollmien-Schlichting waves. iii Preface The second chapter of this thesis has been submitted to a peer-reviewed jour- nal, while the third and fourth chapters are already published journal papers. The contributions of the co-authors are outlined as follows: Chapter 2 has been submitted for publication with myself as the first au- thor and G. A. Lawrence as the second. I was responsible for developing the idea and formulating the theory, validating it using Matlab and writing the manuscript. G. A. Lawrence edited and revised the manuscript. Chapter 3 has been published in Physical Review E. I am the first author of this paper, the second and third authors respectively being M. Rahmani and G. A. Lawrence. I was responsible for developing the idea, writing the codes and preparing the manuscript. M. Rahmani was responsible for analyzing the DNS data, while G. A. Lawrence helped with writing the manuscript as well as providing an overall guidance. Chapter 4 has been published in the Journal of Fluid Mechanics with myself as the first author and I. A. Frigaard as the second. I was responsible for developing the idea, writing the codes and preparing the manuscript. I. A. Frigaard helped me with writing the manuscript and providing mathematical insights. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Hydrodynamic stability theory . . . . . . . . . . . . . . . . . 4 1.2.1 Linearization and normal mode analysis . . . . . . . . 6 1.2.2 Non-modal stability . . . . . . . . . . . . . . . . . . . 7 1.3 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Understanding the mechanism of shear instabilities from wave interaction perspective . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Linear wave(s) at an interface . . . . . . . . . . . . . . . . . . 13 2.2.1 Vorticity waves . . . . . . . . . . . . . . . . . . . . . . 16 v Table of Contents 2.2.2 Interfacial internal gravity waves . . . . . . . . . . . . 18 2.3 Interaction between two linear interfacial waves . . . . . . . . 21 2.4 Homogeneous and stratified shear instabilities . . . . . . . . . 30 2.4.1 The Kelvin-Helmholtz instability . . . . . . . . . . . . 30 2.4.2 The Taylor instability . . . . . . . . . . . . . . . . . . 36 2.4.3 The Holmboe instability . . . . . . . . . . . . . . . . . 40 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3 Non-linear Kelvin-Helmholtz instability in a piecewise linear shear layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2 Non-linear simulation . . . . . . . . . . . . . . . . . . . . . . 51 3.2.1 Contour dynamics . . . . . . . . . . . . . . . . . . . . 53 3.2.2 Direct numerical simulation . . . . . . . . . . . . . . . 58 3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 60 3.3.1 Pre-saturation and saturation phases . . . . . . . . . . 60 3.3.2 Early post-saturation phase . . . . . . . . . . . . . . . 62 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4 Linear stability of a Plane Couette-Poiseuille flow in presence of cross-flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Stability of plane Couette-Poiseuille flow with cross-flow . . . 71 4.2.1 Streamwise Reynolds number . . . . . . . . . . . . . . 77 4.2.2 The stability problem . . . . . . . . . . . . . . . . . . 78 4.2.3 Characteristic effects of varying k and Reinj . . . . . . 80 4.3 PCP flows and the effects of small Reinj . . . . . . . . . . . . 86 4.3.1 Long wavelength approximation . . . . . . . . . . . . 89 4.3.2 Effects of asymmetry of the velocity profile . . . . . . 90 4.3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.4 Intermediate Reinj and short wavelength instabilities . . . . . 97 vi Table of Contents 4.4.1 Behaviour of preferred modes for intermediate Reinj. . 101 4.5 Stability and instability at large Reinj. . . . . . . . . . . . . . 105 4.5.1 Linear energy balance at Reinj,2 . . . . . . . . . . . . 107 4.5.2 Eventual stabilization at Reinj,3 . . . . . . . . . . . . . 112 4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . 121 5.2.1 Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.2.2 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.2.3 Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.3 Future research . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.3.1 Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.3.2 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.3.3 Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . 124 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 vii List of Tables 4.1 Critical values for k = 0.5 and increasing Reinj. . . . . . . . . 88 4.2 Cut-off values, k1, and wavespeed cr,crit, for increasing Reinj. . 90 4.3 Cut-off values evaluated for shorter wavelength instabilities for Re = 106. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 viii List of Figures 1.1 (a) Wall bounded turbulence in channel flow (Green, 2006) and (b) turbulence produced in the ocean mixing layer (Smyth & Moum, 2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Schematic of the interfacial wave interaction mechanism. The deformation and associated vertical velocity of each wave is shown by the same color. Interaction imposes an additional ver- tical velocity (shown by different color). The horizontal arrow associated with a wave indicates the intrinsic wave propaga- tion direction. Both the waves are counter-propagating (move against the background velocity at that location). . . . . . . . 24 2.2 The figure shows that any initial condition (R0,Φ0) finally yields the resonant configuration (RNM ,ΦNM), provided Eq. (2.41) is satisfied. The case depicted here is Kelvin-Helmholtz insta- bility (interaction between two vorticity waves) corresponding to α = 0.4. Any other shear instability will show qualita- tively similar characteristics. (a) Φ versus t corresponding to Φ0 = −pi,−pi/2, 0,ΦNM , pi/2 and pi. The value of R0 is held constant, and is equal to 2. (b) R versus t corresponding to R0 = 1/2, 1(RNM), 2 and 5. The value of Φ0 is held constant, and is equal to −pi/2. . . . . . . . . . . . . . . . . . . . . . . . 27 ix List of Figures 2.3 (a) The setting leading to the Kelvin-Helmholtz instability. The velocity profile in Eq. (2.44) is shown on the left, while the vorticity waves (marked by “V”) are shown on the right. (b) Linear stability diagram of the Kelvin-Helmholtz instability (γ denotes the modal growth rate). . . . . . . . . . . . . . . . . . 31 2.4 Phase portrait of Kelvin-Helmholtz instability corresponding to α = 0.4. The system has two equilibrium points - one unstable (◦) and the other stable (•). Φ is the phase difference between the lower and upper waves, while R represents the ratio of the upper wave amplitude to the lower wave amplitude. . . . . . . 35 2.5 Phase portrait when the two interacting vorticity waves do not produce KH (α = 2). . . . . . . . . . . . . . . . . . . . . . . . 36 2.6 (a) The setting leading to the Taylor instability. The velocity and density profiles in Eq. (2.55) are shown on the left, while the gravity waves “G” are shown on the right. (b) Linear stability diagram of the Taylor instability. The contours represent the growth rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.7 Phase portrait of Taylor instability corresponding to an unsta- ble combination of α and J . Here α = 0.2 and J = 0.7264. . . 41 2.8 (a) The setting leading to the Holmboe instability. The velocity and density profiles in Eq. (2.67) are shown on the left, while the vorticity wave “V” and the gravity wave “G” are shown on the right. (b) Linear stability diagram of the Holmboe instability. The contours represent the growth rate. . . . . . . . . . . . . . 42 2.9 Phase portrait of Holmboe instability corresponding to α = 1 and J = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 x List of Figures 3.1 (a) 3D DNS performed to capture the complete turbulent dissi- pation of a KH billow ensuing from a hyperbolic tangent velocity profile. False color is added to aid visualization. (b) Spanwise averaged mean velocity profile corresponding to each instant. (c) The magenta line is the continuous velocity profile obtained from Eq. (3.7), while the thick black line below it is the piece- wise linear profile from Eq. (2.44). (d) CVWs (exaggerated) existing at the vorticity discontinuities. . . . . . . . . . . . . 52 3.2 Contour dynamics simulation of a piecewise linear shear layer by Pozrikidis & Higdon (1985) for (a) α = 0.639, (b) α = 0.5 and (c) α = 0.0625. . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Time evolution of Kelvin-Helmholtz instability - comparison be- tween DNS and CD. . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Temporal variation of the wave amplitude a. The straight blue line is the prediction from linear theory, the black line corre- sponds to CD while the magenta line represents DNS results. The green and red lines in the inset respectively show the vari- ation of the ellipse aspect ratio r and the angular rotation rate ω with time. These variations are obtained by solving Eqs. (3.9)-(3.10). The black markers indicate the corresponding data points measured from DNS. . . . . . . . . . . . . . . . . . . . 60 3.5 Formation of winding filaments around the elliptical vortex dur- ing the late post-saturation phase. . . . . . . . . . . . . . . . . 65 4.1 Mean Velocity Distribution for k = 0.5 and Reinj = 0, 1, 4, 8, 15 and 30 (Reinj = 30 marked with a ) . . . . . . . . . . . . 76 4.2 (a) Maximal velocity gradient, |Du|max, plotted against Reinj for k = 0.35, 0.5, 0.65, (k = 0.65 marked with a ). The thick line indicates where the maximum is attained at y = −1; otherwise at y = 1. (b) Variation of D2u with y for k = 0.5 for : Reinj = 0, (); Reinj = 4, (◦); Reinj = 8, (×); Reinj = 12, (). 77 xi List of Figures 4.3 Eigenspectrum of (α,Re) = (1, 6000) by varying k andReinj. 40 least stable modes are considered. (a) Effect of increasing k from 0 to 1 in steps of 0.01, keeping Reinj = 0. (b) Effect of in- creasing Reinj from 0 to 100 in steps of 0.05, keeping k = 0 (PP flow). The symbols in (a) and (b) are similar and are denoted as follows: k = 0 or Reinj = 0 by (), k = 1 or Reinj = 100 by (◦) and intermediate k or Reinj by (.). Note that the PP flow spectrum is represented by the  in both figures, and shows the vertical family of S-modes, the branch of A-modes (diagonally upwards from centre to left) and branch of P-modes (diagonally upwards from centre to right). . . . . . . . . . . . . . . . . . 82 4.4 (a) Effect of increasing Reinj on the stability of PCP flow, for (α,Re) = (1, 6000) and different values of k = 0 (), 0.5 (◦), 1 (×). (b) Maximal growth rate for increasing Reinj at different Re, (k = 0.5 and the step in values of Re between curves is 104). 83 4.5 Maximal growth rate versus Reinj at Re = 40000: (a) Reinj,2 & Reinj,3 for k = 0 (◦) to 1 (); (b) Reinj,1 for k = 0 (◦) to 0.6 (). Step size is 0.2 in both figures. . . . . . . . . . . . . . . . . . . 85 4.6 Critical values for k = 0.5: (a) neutral stability curves for Reinj = 0 (×), 0.3 (◦) and 0.53 (); (b) variation in cr,crit with Reinj. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.7 (a) Long wave NSC’s showing the dependence of λ on k for Reinj = 0 (dashed line), 0.3 (dash-dot), 0.5 (solid), 0.7 (dash- dot-dot) and 1 (long dash); (b) ci,crit versus k for λ = 2.5×10−5, and Reinj = 0 (dashed line), 0.3 (dash-dot), 1 (solid), 1.2 (dash- dot-dot) and 1.3 (long dash). . . . . . . . . . . . . . . . . . . 89 xii List of Figures 4.8 Eigenspectrum for (k, α, Re) = (0.5, 0.2, 31656) (a) Reinj = 0.5 (Critical Conditions) and (b) Reinj = 23.5. Symbol ◦ indi- cates the eigenspectrum from the O-S equation while  indi- cates the spectrum obtained by neglecting the additional cross- flow inertial term. . . . . . . . . . . . . . . . . . . . . . . . . 91 4.9 k1 as a function of Reinj,1 (shown by ) and the radius ratio, η (shown by ) in ACP flow (Sadeghi & Higgins (1991)) . . . . 93 4.10 Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to criticality at Reinj= (a) 0, (b) 0.2, (c) 0.4 and (d) 0.6. In all the cases, k = 0.5. Dash- dot-dot line with symbol  represents T2, dashed line with filled 4 represents 1 Re T3 and solid vertical line represents the location of the critical layer. . . . . . . . . . . . . . . . . . . . . . . . . 96 4.11 Neutral Stability Curves (NSCs) for (a) k = 0 and (b) 0.18 at different Reinj. The symbols indicate different values of Reinj and are as follows: × → Reinj = 0, ◦ → Reinj = 6 in (a) and 4 in (b), → Reinj = 12 in (a) and 8 in (b) . . . . . . . . . . . 98 4.12 Shorter wavelength cut-off showing k1 as a function ofReinj,1. The flow is linearly stable for Re ≤ 106 above the curve. The values in Table 4.3 are marked by . . . . . . . . . . . . . . . . . . . 100 xiii List of Figures 4.13 Behavior of preferred modes (belonging to different wavelengths and denoted by alphabets ‘A’-‘D’) under the influence of cross- flow with Re = 106. Symbols  and ◦ respectively imply the starting and the ending position of the preferred mode in the ci, cr plane, whereas the dots (‘.’) trace the locus. The difference in Reinj between consecutive dots is 0.1. (a) k = 0. Mode ‘A’ has α = 0.001 and is traced for Reinj=[0,25]. Mode ‘B’ has α = 3.5227 and is traced for Reinj=[15,21] (shown in the inset), the position at Reinj = 15 is marked by ‘∗’. (b) k = 0.5. Mode ‘C’ has α = 0.01 and is traced for Reinj=[0,30]. Mode ‘D’ has α = 2.5 and is traced for Reinj=[7,30]. . . . . . . . . . . . . . 103 4.14 Isovalues of the normalised perturbation stream functions (ψ′) for the preferred modes at Re = 106 under different Reinj. The streamwise extent of the domain is one wavelength. Corre- sponding to k = 0.5, the long wavelength mode ‘C’ is shown for (a) Reinj = 0.1 (unstable), (b) Reinj = 1 (stable) and (c) Reinj = 30 (unstable). Corresponding to k = 0, ψ ′ for two dif- ferent preferred modes, viz. ‘A’ and ‘B’ are shown. The shorter wavelength mode ‘B’ (α = 3.5227) is shown for (d) Reinj = 15 (unstable) and (e) Reinj = 21 (stable). The longer wavelength mode ‘A’ (α = 0.001) is shown for (f) Reinj = 25 (unstable). . 104 4.15 (a) NSC of PP flow (k = 0) when Reinj → Re−inj,2. The different values of Reinj are 22.5 (dashed line with ×), 23 (dash-dot line with ) and 24 (dash-dot-dot line with ◦). Near cut-off, αRe is constant along the upper and lower branches. (b) Long wave NSCs showing the dependence of log10 λ on k. The different values of Reinj are 22.4 (), 23.5 (◦), 24 () and 25 (×). Cut- off is achieved over the entire range of k, i.e. [0, 1]. . . . . . . . 106 xiv List of Figures 4.16 (a) Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to criticality at Reinj = 25. Dash-dot-dot line with symbol  represents T2, dashed line with filled 4 represents 1 Re T3 and solid vertical line represents the location of the critical layer. (b) Reynolds Stress τ dis- tribution at criticality for Reinj = 0 (denoted by  symbol), Reinj = 0.6 (denoted by ×) and Reinj = 25 (denoted by ◦). The location of the critical layers are shown by solid lines with corrosponding symbols. . . . . . . . . . . . . . . . . . . . . . . 107 4.17 Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to mode ‘C’ at Reinj= (a) 0, (b) 0.6, (c) 1, (d) 3, (e) 10 and (f) 30. Dash-dot-dot line represents T2, solid line represents 1 Re T3. . . . . . . . . . . . . 109 4.18 Non-dimensional mean perturbation kinetic energy profiles for mode ‘C’ at different Reinj. Solid lines with symbols denote the unstable modes. For each Reinj, q has been scaled by its maximum value. . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.19 Variation of k with Reinj. The filled  symbols show the long wavelength cut-off achieved for 0.7 ≥ k ≥ 0.19. The filled ◦ symbols show the shorter wavelength cut-off for 0.19 > k ≥ 0 evaluated numerically for Re = 106. The filled  symbols imply the second long wavelength cut-off. The shaded region depicts the entire zone of unconditional linear stability. . . . . . . . . 115 5.1 Schematic of a vorticity interface - left half shows unperturbed velocity field, while the right half depicts infinitesimal interfacial displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 xv Acknowledgments I humbly acknowledge the two major sources of financial support that I re- ceived during my Ph.D. studies, the Four Year Fellowship of UBC and the Faculty of Applied Science Graduate Award. I would like to thank my Ph.D. supervisor, Prof. Greg Lawrence for his continuous encouragement, support and invaluable guidance. I am indebted to him for giving me the research freedom, and always having confidence in my abilities. This thesis would not be anywhere near the product it currently is if I was not given the opportunity to grow as an independent researcher. I would also like to thank my supervisory committee members - Bernard Laval, Neil Balmforth and Ian Frigaard for their help and valuable comments. The general atmosphere of UBC has provided a great learning experience. I took many courses from other Departments, e.g. Mechanical, Mathematics, Physics and Earth and Ocean Sciences (again, thanks to Greg for giving me this freedom), and during this time I had the opportunity to interact with many great minds scattered across the campus. I took a course on Hydrodynamic stability with Prof. Ian Frigaard; the course project actually became the fourth chapter of my Ph.D. thesis. I am indebted to Prof. Frigaard for his help and guidance, and invaluable advise on academia in general. I am also indebted to xvi Acknowledgments Prof. Neil Balmforth for being kind enough to lend his valuable time whenever needed. The long discussions we had on waves and instabilities have helped me a lot with my research. The Environmental Fluid Mechanics weekly lunch meetings also played a key role in my academic development - it has been a great peer learning experience. I would like to thank our EFM members Mona Rahmani, Ted Tedford, Jeff Carpenter and Andrew Hamilton, as well as Sarah Hormozi from the Institute of Applied Mathematics for their encouragement, help and support. Last but not the least, I am indebted to Prof. Kraig Winters of Scripps Institute of Oceanography. I was lucky to have him as my external examiner, and the advise I had from him is simply invaluable. Special thanks for my wife, Tanaya, who is also a Ph.D. student at UBC. Long ago, when I just completed my Undergraduate degree and had no in- tention of pursuing higher studies, she made me believe that I am a “Ph.D. material”. She has always been there by my side, and has kept on motivating and inspiring me everyday. I am also thankful to my Indian friends in Van- couver - Santanu Mitra and Anuradha Mitra, whose love, encouragement and hospitality I greatly treasure. Finally, I would like to thank my mother from the bottom of my heart. She undertook a lot of difficulties for my upbringing. I am what I am because of her. My journey as a Ph.D. student has been an eventful and satisfying ex- perience. I believe I have grown as an academic, a philosopher, and most importantly, as a human being. xvii Dedication To the three lovely ladies - my mother, my wife, and my motherland... xviii Chapter 1 Introduction It was easier to know it than to explain why I know it. — Sherlock Holmes in A Study in Scarlet. Velocity gradient, also known as velocity shear (or simply, shear), is ubiq- uitous in natural and industrial flows. As an example, when a fluid flows past a stationary obstacle, its velocity varies from zero at the obstacle wall to the free stream velocity. Therefore a shear layer is produced adjacent to the wall, which is technically known as the boundary layer. Flows over topography or over aircraft wings are examples of boundary layer flows. Shear flow is not necessarily wall bounded. In atmospheric and oceanic flows there are situations where the upper fluid is left (right) moving while the lower fluid is moving to the right (left). This is known as free shear flow and is typically observed in exchange dominated flows, e.g. oceanic gravity currents, estuarine flows and sea breezes. Natural exchange flows are often density strat- ified, hence understanding such shear flow dynamics requires consideration of the effects of density variations. Shear flows can become hydrodynamically unstable, resulting in a process characterized by the growth of wavelike perturbations. These perturbations 1 1.1. Governing equations can grow at an exponential rate, transforming the base flow from a laminar to turbulent state. Hence hydrodynamic instability is often regarded as the precursor of turbulence. During the instability mechanism energy is extracted from the large-scale motions and transferred to the smaller scales. The smallest scale processes are responsible for dissipating the mechanical energy into heat. If the flow is density stratified, the resulting mechanism leads to the mixing of the density field. Turbulence and mixing are important aspects of many problems in me- teorology, oceanography and several branches of engineering (e.g. mechanical, aerospace, environmental, chemical). In the aerospace industry, understanding boundary layer instabilities is crucial for reducing turbulent drag over an air- craft. In the fields of oceanography and atmospheric sciences, understanding free shear instabilities is important for developing accurate models for predict- ing weather and climate. Fig. 1.1 shows examples of wall bounded (boundary layer) and free shear turbulence. 1.1 Governing equations In this thesis, we will only be considering two-dimensional flows. Let x be the streamwise direction and z be the cross-streamwise/vertical direction1. We as- sume the fluid to be incompressible and non-diffusive. Hence the conservation 1We replace z by y in Chapters 3 and 4, since the flows under consideration are most likely to occur in the horizontal plane. 2 1.1. Governing equations Figure 1.1: (a) Wall bounded turbulence in channel flow (Green, 2006) and (b) turbulence produced in the ocean mixing layer (Smyth & Moum, 2012). of volume (continuity equation) can be written as ∂u ∂x + ∂w ∂z = 0 (1.1) where (u,w) is the velocity field. Let us denote the density, pressure and kine- matic viscosity by ρ, P , and ν respectively. Furthermore, we will only consider flows where density variations are either small or negligible. The effect of small density variations is included through Boussinesq approximations (Kundu & Cohen, 2004), in which it is assumed that the variation in density is negligible except when ρ is multiplied by gravity, g. Following Kundu & Cohen (2004), 3 1.2. Hydrodynamic stability theory the two-dimensional, Boussinesq momentum equations are written as follows: ∂u ∂t + u ∂u ∂x + w ∂u ∂z = − 1 ρ0 ∂P ∂x + ν4u (1.2) ∂w ∂t + u ∂w ∂x + w ∂w ∂z = − 1 ρ0 ∂P ∂z − g ρ ρ0 + ν4w (1.3) where ρ0 is a reference density, and 4 ≡ ∂2/∂x2 +∂2/∂z2. The fluid, assumed non-diffusive, also satisfies the conservation of mass: ∂ρ ∂t + u ∂ρ ∂x + w ∂ρ ∂z = 0 (1.4) 1.2 Hydrodynamic stability theory To understand the hydrodynamic stability of a given fluid flow, we follow the procedure outlined in Drazin & Reid (2004). We will be considering a background flow U(z) (i.e. parallel to x and varying only in z) and background density ρ̄(z). The background flow is perturbed with a small two-dimensional perturbation velocity of (u ′ , w ′ ), density of ρ ′ , and pressure of P ′ , so that the total two-dimensional velocity field, (u,w), density field, ρ, and pressure field, P , are u = U + u ′ , w = w ′ , ρ = ρ̄+ ρ ′ , P = P̄ + P ′ , (1.5) where P̄ is the hydrostatic background pressure2, implying ∂P̄ /∂z = −ρ̄g. 2The quantity ∂P̄ /∂x = 0 in Chapters 2 and 3, whereas in Chapter 4 ∂P̄ /∂x = negative constant. 4 1.2. Hydrodynamic stability theory Using these relations along with Eq. (1.5) in Eqs. (1.2)-(1.3), we obtain ∂u ∂t + u ∂u ∂x + w ∂u ∂z = − 1 ρ0 ∂P ′ ∂x + ν4u (1.6) ∂w ∂t + u ∂w ∂x + w ∂w ∂z = − 1 ρ0 ∂P ′ ∂z − g ρ ′ ρ0 + ν4w (1.7) We can use the continuity equation (Eq. (1.1)) to define a streamfunction such that u = ∂ψ ∂z ; w = −∂ψ ∂x (1.8) Taking ∂(1.6)/∂z - ∂(1.7)/∂x and using Eq. (1.8), we formulate the vorticity equation: ∂4ψ ∂t + u ∂4ψ ∂x + w ∂4ψ ∂z = g ρ0 ∂ρ′ ∂x + ν44ψ (1.9) We define a characteristic length scale l, velocity scale δU and density scale δρ. This allows us to represent Eq. (1.9) in complete non-dimensional form: ∂4ψ ∂t + u ∂4ψ ∂x + w ∂4ψ ∂z = J ∂ρ′ ∂x + 1 Re 44ψ (1.10) where J is the bulk Richardson number and Re is the Reynolds number, and are defined as follows: J = δρ gl ρ0 (δU) 2 ; Re = δU l ν (1.11) 5 1.2. Hydrodynamic stability theory 1.2.1 Linearization and normal mode analysis We decompose ψ into mean and perturbation components: ψ = Ψ(z) + ψ′, and substitute it into Eq. (1.10). The equation is then linearized about the background state, i.e. terms with product of perturbation components are neglected. Similar treatment is also done on Eq. (1.4). The conventional way of solving the resulting linear equation is by applying the method of normal modes. In this method we assume the perturbation streamfunction and perturbation density to have the from ψ ′ (x, z, t) = ϕ(z)eiα(x−ct) and ρ ′ (x, z, t) = ρ̂(z)eiα(x−ct) respectively, where i is the imaginary unit, α is the real wavenumber, and c = cr + ici is the complex phase speed of the mode. The quantities αci and cr respectively denote the exponential growth rate and phase speed of the mode. When ci > 0 the mode is unstable. Substitution of normal mode type perturbations into linearized Eqs. (1.10) and (1.4) yield the following eigenvalue problem-set: [(c− U)(α2 −D2)−D2U ]ϕ = Jρ̂− i αRe (α2 −D2)2ϕ (1.12) (U − c)ρ̂ = ϕDρ̄ (1.13) where D ≡ d/dz. The above two equations can be combined to produce the stratified, viscous linear stability equation: [(c− U)(α2 −D2)−D2U ]ϕ = J Dρ̄ (U − c)ϕ− i αRe (α2 −D2)2ϕ (1.14) 6 1.2. Hydrodynamic stability theory If the flow is assumed inviscid (Re → ∞), Eq. (1.14) produces the Taylor- Goldstein equation (Taylor, 1931; Goldstein, 1931): D2ϕ+ [ −J Dρ̄ (U − c)2 − D2U (U − c) − α 2 ] ϕ = 0 (1.15) If the flow is homogeneous and viscous, Eq. (1.14) produces the Orr-Sommerfield equation (Drazin & Reid, 2004): iαRe[(c− U)(α2 −D2)−D2U ]ϕ = (α2 −D2)2ϕ (1.16) 1.2.2 Non-modal stability In the above section, we have introduced the concept of normal-mode based linear stability theory. Although this theory is extremely useful, it provides limited insight into the physical mechanism(s) responsible for hydrodynamic instability. The answer to why an infinitesimal perturbation vigorously grows from a stable background flow is provided in the form of non-intuitive mathe- matical theorems - Rayleigh’s inflection point theorem and Fjørtoft’s extension for the case of homogeneous flows, and Miles-Howard criterion for stratified flows (Drazin & Reid, 2004). Since linear instability is the first step towards understanding the more complicated and highly elusive non-linear processes like chaos and turbulence, it is helpful to formulate alternative theories which are able to provide intuitive explanations. Another shortcoming of the normal mode stability theory is the normal-mode assumption itself. The extensive 7 1.2. Hydrodynamic stability theory work by Farrell (1984), Trefethen et al. (1993), Schmid & Henningson (2001) and others have shown that shear allows rapid non-modal transient growth due to non-orthogonal interaction between the modes. Farrell & Ioannou (1996) developed the Generalized stability theory for linear dynamical systems, and showed the process of obtaining the optimal non-modal growth from a singular value decomposition of the propagator matrix of the linear dynamical system. Theories have been proposed to provide an intuitive understanding of the hydrodynamic instability process. Probably the first mechanistic picture of stratified shear instabilities was provided by Holmboe (1962). Using idealized velocity and density profiles, Holmboe postulated that the resonant interaction between stable propagating waves, each existing at a discontinuity in the back- ground flow profile (density profile discontinuity produces gravity waves and vorticity profile discontinuity produces vorticity waves), yields exponentially growing instabilities. He was able to show that Kelvin-Helmholtz instability (Rayleigh, 1880) is the result of the interaction between two vorticity waves (also known as Rayleigh waves). Moreover, Holmboe also found a new type of instability, later came to be known as the “Holmboe instability”, produced by the interaction between vorticity and gravity waves. Bretherton, a contempo- rary of Holmboe, proposed a similar theory to explain mid-latitude cyclogenesis (Bretherton, 1966). He hypothesized that the cyclones form due to a baro- clinic instability mechanism caused by the interaction between two Rossby edge waves (vorticity waves in a rotating frame of reference), one existing at the earth’s surface and the other located at the atmospheric tropopause. 8 1.3. Overview The theory proposed by Holmboe and Bretherton has been refined and re- interpreted over the years, see Cairns (1979); Hoskins et al. (1985); Caulfield (1994); Baines & Mitsudera (1994); Heifetz et al. (1999); Carpenter et al. (2013). As reviewed in Carpenter et al. (2013), resonant interaction between two edge waves in an idealized homogeneous or stratified shear layer occurs when these waves attain a phase-locked state, i.e. they are at rest relative to each other. Maintaining this phase-locked configuration, the waves grow equally at an exponential rate. Hoskins et al. (1985) summarized the Rossby edge wave induced shear instability mechanism in one sentence: ‘The induced velocity field of each Rossby wave keeps the other in step, and makes the other grow.’ 1.3 Overview In this thesis we have investigated various aspects of two dimensional shear flow instabilities. Chapter 2 is dedicated to providing a physical understand- ing of shear instabilities by looking into the problem from wave interaction perspective. The physical reason behind the growth of normal-mode type (i.e. exponentially growing) instabilities is obtained from this analysis. The inves- tigation also provides valuable insight into the non-modal stability theory. The idea proposed in Chapter 2 is extended to the non-linear regime in Chapter 3 for the case of Kelvin-Helmholtz instability. It has been shown that the non-linear interaction between two counter-propagating vorticity waves 9 1.3. Overview produce elliptical vortex patches connected via thin braids. The vortex and braid dynamics have been investigated. In Chapter 4 we have followed the conventional normal-mode approach to investigate how channel flows can be unconditionally stabilized/destabilized. In a parameter space of non-dimensional wall velocity and non-dimensional cross flow velocity, we have investigated how these two parameters influence the stability of a channel flow. It has been shown that the two effects compensate each other. Moreover, a small range of cross-flows have been shown to exist for which the channel flow is unconditionally linearly stable even in the absence of wall motion. 10 Chapter 2 Understanding the mechanism of shear instabilities from wave interaction perspective 1 Before turning to those moral and mental aspects of the matter which present the greatest difficulties, let the inquirer begin by mas- tering more elementary problems. — Sherlock Holmes in A Study in Scarlet. 2.1 Introduction In this Chapter we study shear instabilities in terms of interacting interfacial waves. Holmboe was probably the first to postulate that the resonant inter- action between stable progressive waves, each existing at a discontinuity in the background flow profile (can be vorticity or density discontinuity), yields 1A version of this chapter has been submitted for publication. A. Guha and G. A. Lawrence (2013) A wave interaction approach to studying non-modal homogeneous and stratified shear instabilities. 11 2.1. Introduction exponentially growing instabilities. Holmboe showed that Kelvin-Helmholtz instability (Rayleigh, 1880) is the result of the interaction between two vor- ticity waves (also known as Rayleigh waves). Moreover, Holmboe also found a new type of instability, later came to be known as the “Holmboe instabil- ity”, produced by the interaction between vorticity and gravity waves. Holm- boe’s hypothesis has been refined and re-interpreted over the years, see Cairns (1979); Hoskins et al. (1985); Caulfield (1994); Baines & Mitsudera (1994); Heifetz et al. (1999); Carpenter et al. (2013). In recent years, Heifetz and co-authors (Heifetz et al., 1999, 2004; Heifetz & Methven, 2005) have extensively studied the interaction between Rossby edge waves. By not limiting the Rossby edge waves to be of the normal-mode type, they were able to obtain non-modal instability and transient growth mecha- nisms. While Heifetz et al. (2004) derived the governing equations using the Hamiltonian approach, Heifetz & Methven (2005) used the streamfunction- vorticity approach. Their successful attempt has motivated us to formulate a generalized interfacial wave interaction technique for studying homogeneous and stratified shear instabilities. Without forcing the wave to be of the normal- mode type, and furthermore, without assuming any particular type of wave- form (e.g. gravity wave or vorticity wave), we have formulated the governing equations for wave interaction. Our equations are derived from the linearized kinematic and dynamic (for stratified flows) conditions. Unlike Heifetz et al. (2004) and Heifetz & Methven (2005), our key variables are the vertical dis- placement and the vertical velocity of the wave. The choice of variables, along 12 2.2. Linear wave(s) at an interface with the derivation methodology, provide a deeper and probably more intuitive understanding of the wave interaction process. 2.2 Linear wave(s) at an interface In the present study we consider multi-layered flows with constant density and vorticity in each layer. This configuration makes the equations of mo- tion for perturbations within a layer to be the same as that in an irrotational background flow. The interface between two adjacent layers signifies a discon- tinuity in vorticity or density. The former is a vorticity interface, while the latter is a density interface. Let such an interface existing at a location z = zi be perturbed by an infinitesimal vertical displacement ηi, given as follows: ηi = Aηi(t) cos[αx+ φηi(t)] (2.1) This displacement manifests itself in the form of stable, progressive wave(s), the amplitude and phase of which are Aηi and φηi respectively. While a vortic- ity interface produces a vorticity wave, two oppositely traveling gravity waves are produced in the case of a density interface. We have assumed the interfa- cial displacement (or the wave) to be monochromatic, having a wavenumber α. Moreover, the interface satisfies the kinematic condition - a particle initially on the interface will remain there forever. The linearized kinematic condition 13 2.2. Linear wave(s) at an interface is given by ∂ηi ∂t + Ui ∂ηi ∂x = wi (2.2) where Ui ≡ U(zi) is the background velocity and wi is the vertical velocity at the interface. We prescribe the latter to be as follows: wi = Awi(t) cos[αx+ φwi(t)] (2.3) Here Awi is the amplitude and φwi is the phase of wi. The interfacial deforma- tion sets up an irrotational velocity field away from the interface (Holmboe, 1962; Caulfield, 1994): 4ψ′ = 0 when z 6= zi (2.4) When the perturbation streamfunction of the form ψ ′ (x, z, t) = ϕ(z)ei(αx+θψ(t)) is substituted in the above equation we get ∂2ϕ ∂z2 − α2ϕ = 0 (2.5) The above equation yields ϕ = e−α|z−zi|ϕi. The vertical velocity w = −∂ψ∂x is then given by w = e−α|z−zi|wi (2.6) Substituting Eqs. (2.1) and (2.3) in Eq. (2.2), we obtain Ȧηi cos (αx+ φηi)−Aηi ( αUi + φ̇ηi ) sin (αx+ φηi) = Awi cos (αx+ φwi) (2.7) 14 2.2. Linear wave(s) at an interface Here φηi , φwi ∈ [−pi, pi]. The frequency and the growth rate of a wave are respectively defined as σi ≡ −φ̇ηi and γi ≡ Ȧηi/Aηi (overdot denotes d/dt). Using these definitions in Eq. (2.7), we get σi = αUi −ωi sin (4φii) (2.8) γi = ωi cos (4φii) (2.9) where4φii ≡ φwi−φηi . In order to get Eqs. (2.8)-(2.9) from Eq. (2.7), we write the R.H.S. of Eq. (2.7) as follows: Awi cos (αx+ φwi) = Awi cos (αx+ φηi +4φii). The cosine function is expanded using a standard trigonometric identity. Fi- nally we collect the coefficients of sin(αx) and cos(αx). Eq. (2.8) shows that the frequency of a wave consists of two components - (i) the Doppler shift αUi due to the background velocity, and (ii) the intrinsic frequency −ωi sin (4φii), where ωi ≡ Awi/Aηi . The phase-speed ci ≡ σi/α of the wave is found to be ci = Ui − ωi α sin (4φii) (2.10) The last term (including the negative sign) denotes the intrinsic phase-speed. Noting that a wave in isolation cannot grow or decay on its own, Eq. (2.9) demands that |4φii| = pi/2. Therefore for a stable wave, the vertical velocity field at the interface has to be in quadrature with the interfacial deformation. According to Eqs. (2.8) and (2.10), the quadrature condition makes the mag- 15 2.2. Linear wave(s) at an interface nitude of the intrinsic frequency and the intrinsic phase-speed to be ωi and ωi/α respectively. The intrinsic direction of motion of the wave, however, is determined by 4φii. For waves moving to the left relative to the interfacial velocity Ui, 4φii = pi/2. Similarly for right moving waves, 4φii = −pi/2. When such a stable, progressive wave is acted upon by external influence(s) (e.g. when another wave interacts with the given wave, as detailed in §2.3), the quadrature condition is no longer satisfied, i.e. |4φii| 6= pi/2. Therefore, the wave may grow (γi > 0) or decay (γi < 0), and its intrinsic frequency and phase-speed may change. In our analyses we will consider two types of progressive interfacial waves - vorticity waves and internal gravity waves. 2.2.1 Vorticity waves Vorticity waves, also known as Rayleigh waves, exist at a vorticity interface (i.e. regions involving a sharp change in vorticity). Such interfaces are a common feature in the atmosphere and oceans. In rotating frame, the analogue of vorticity wave is the Rossby edge wave which exists at a sharp transition in the potential vorticity. When Rossby edge waves propagate in a direction opposite to the background flow, they are called “counter-propagating Rossby waves” or CRWs (Heifetz et al., 1999). In order to evaluate the frequency ωi of vorticity waves, let us consider a 16 2.2. Linear wave(s) at an interface velocity profile having the form U (z) =  Ui z ≥ zi Sz z ≤ zi (2.11) Here the constant S = Ui/zi is the vorticity, or the shear in the region z ≤ zi (Carpenter et al., 2013). In Eq. (2.11), the vorticity dU/dz is discontinuous at z = zi. This condition supports a vorticity wave. A deformation ηi of the interface adds vorticity S to the upper layer and removes it from the lower layer, creating a vorticity imbalance in the horizontal direction. This imbalance produces a velocity field which causes the wave to propagate in the horizontal direction. The horizontal component ui of the perturbation velocity field set up by the interfacial deformation undergoes a jump at the interface, the value of which can be determined from Stokes’ Theorem (Appendix 5.3.3): u+i − u−i = Sηi (2.12) By taking an x derivative of Eq. (2.12) and invoking the continuity relation, we get −∂w + i ∂z + ∂w−i ∂z = S ∂ηi ∂x (2.13) By substituting Eqs. (2.1) and (2.3) in Eq. (2.13), we obtain ωi = −S 2 sin (αx+ φηi) cos (αx+ φwi) = S 2 sin (4φii) (2.14) 17 2.2. Linear wave(s) at an interface The fact that ωi is always positive demands 4φii = pi 2 sgn(S) (2.15) where sgn( ) is the sign function. From Eq. (2.14), the intrinsic frequency of a vorticity wave is found to be −S/2. The phase-speed ci can be evaluated by substituting Eq. (2.14) into Eq. (2.10): ci = Ui − S 2α (2.16) If S > 0, the vorticity wave moves to the left relative to the background flow. Alternative derivations of the frequency and phase-speed of a vorticity wave can be found in several references, e.g. Caulfield (1994), Sutherland (2010) or Carpenter et al. (2013). 2.2.2 Interfacial internal gravity waves Interfacial gravity waves exist at a density interface, i.e. regions involving sharp change in density. The most common example of interfacial gravity wave is the surface wave existing at the interface of air and water. Here we will be considering interfacial internal gravity waves (hereafter, gravity waves) only. Such waves exist in density stratified flows having a thin density interface (pycnocline). Since most natural water bodies like lakes, estuaries and oceans are density stratified, gravity waves are ubiquitous. 18 2.2. Linear wave(s) at an interface In the case of gravity waves, the intrinsic frequency −ωi sin (4φii) can be evaluated by considering the dynamic condition. The latter implies that the pressure at the density interface must be continuous. Let the density of upper and lower fluids be ρ1 and ρ2 respectively. Then the linearized dynamic condition at the interface after some simplification becomes (Eq. (3.13) of Caulfield (1994)) : ∂ψi ∂t + Ui ∂ψi ∂x = g ′ 2α ∂ηi ∂x (2.17) Here g ′ ≡ g(ρ2 − ρ1)/ρ0 is the reduced gravity and ρ0 is the reference density. Under Boussinesq approximation ρ0 ≈ ρ1 ≈ ρ2. By taking an x derivative of Eq. (2.17) and using the streamfunction relation {ui, wi} = {−∂ψi/∂z, ∂ψi/∂x}, we get ∂wi ∂t + Ui ∂wi ∂x = g ′ 2α ∂2ηi ∂x2 (2.18) Substitution of Eqs. (2.1) and (2.3) in Eq. (2.18) yields Ȧwi cos (αx+ φwi)− Awiφ̇wi sin (αx+ φwi) −αUiAwi sin (αx+ φwi) = − g ′ α 2 Aηi cos (αx+ φηi) (2.19) The quantity φ̇wi = φ̇ηi = −σ = −αci. On substituting this relation in Eq. (2.19), we obtain ωi = g ′ 2 sin (4φii) (Ui − ci) (2.20) 19 2.2. Linear wave(s) at an interface Since ωi is a positive quantity, Eq. (2.20) demands 4φii = pi 2 sgn(Ui − ci) (2.21) An important aspect of Eq. (2.20) is that it has been derived independent of the kinematic condition. The presence of single or multiple interfaces does not alter the expression in Eq. (2.20), implying that this equation provides a generalized description of ωi. Inclusion of kinematic condition yields an expression for ωi which is simpler but problem specific. For example, when a single interface is present, inclusion of kinematic condition in Eq. (2.20) produces the well known expression for gravity wave frequency. Substituting Eq. (2.8) (this equation has been derived from the kinematic condition for a single interface, i.e. Eq. (2.2)) in Eq. (2.20) and considering only the positive value, we obtain the dispersion relation for gravity waves: ωi = √ g′α 2 (2.22) Moreover Eq. (2.9) requires |4φii| = pi/2. Substituting it along with Eq. (2.20) in Eq. (2.10)) produces the well known expression for the phase-speed of a gravity wave: ci = Ui ± √ g′ 2α (2.23) The above equation shows that each density interface supports two gravity waves, one moving to the left and the other to the right relative to the back- 20 2.3. Interaction between two linear interfacial waves ground velocity Ui. Alternative approaches to deriving the frequency and phase-speed of a gravity wave can be found in Caulfield (1994); Sutherland (2010); Carpenter et al. (2013)). 2.3 Interaction between two linear interfacial waves Let us now consider a system with two interfaces, one at z = z1 and the other one at z = z2. The linearized kinematic condition at each of these interfaces is then given by: ∂η1 ∂t + U1 ∂η1 ∂x = w1 + e −α|z1−z2|w2 (2.24) ∂η2 ∂t + U2 ∂η2 ∂x = e−α|z1−z2|w1 + w2 (2.25) It has been implicitly assumed that both waves have the same wavenumber. The R.H.S. of Eqs. (2.24)-(2.25) reveal the subtle effect of wave interaction, and can be understood as follows. The effect of w1 extends away from the interface z1, hence it can be felt by a wave existing at another location, say z2. Therefore the vertical velocity of the wave at z2 gets modified - it becomes the linear superposition of its own vertical velocity w2 and the component of w1 existing at z2. This phenomenon is also known as “action-at-a-distance”, see Heifetz & Methven (2005). 21 2.3. Interaction between two linear interfacial waves On substituting Eqs. (2.1) and (2.3) in Eqs. (2.24)-(2.25), we get Ȧη1 cos (αx+ φη1)− Aη1 ( αU1 + φ̇η1 ) sin (αx+ φη1) = Aw1 cos (αx+ φw1) + e −α|z1−z2|Aw2 cos (αx+ φw2) (2.26) Ȧη2 cos (αx+ φη2)− Aη2 ( αU2 + φ̇η2 ) sin (αx+ φη2) = e−α|z1−z2|Aw1 cos (αx+ φw1) + Aw2 cos (αx+ φw2) (2.27) Proceeding in a manner similar to §2.2, the growth rate γi and phase-speed ci of each wave is found to be γ1 = Aw1 Aη1 cos (4φ11) + Aw2 Aη1 e−α|z1−z2| cos (4φ12) (2.28) c1 = U1 − 1 α [ Aw1 Aη1 sin (4φ11) + Aw2 Aη1 e−α|z1−z2| sin (4φ12) ] (2.29) γ2 = Aw2 Aη2 cos (4φ22) + Aw1 Aη2 e−α|z1−z2| cos (4φ21) (2.30) c2 = U2 − 1 α [ Aw2 Aη2 sin (4φ22) + Aw1 Aη2 e−α|z1−z2| sin (4φ21) ] (2.31) Here 4φij ≡ φwj − φηi . When α |z1 − z2| → ∞, the two waves decouple, and we recover Eqs. (2.9)-(2.10) for each wave. As argued in §2.2, a wave in isolation cannot grow or decay on its own. Therefore, the first term in each of Eq. (2.28) and Eq. (2.30) should be equal to zero, implying |4φii| = pi/2. In all our analyses, we will be considering a system with a left moving top wave (4φ11 = pi/2) and a right moving bottom wave (4φ22 = −pi/2), the wave motion being relative to the background velocity at the corresponding 22 2.3. Interaction between two linear interfacial waves interface. Let the phase shift between the bottom and top waves be Φ ≡ φη2−φη1 . Therefore4φ12 = Φ−pi/2 and4φ21 = pi/2−Φ. Defining amplitude ratio R ≡ Aη1/Aη2 , we re-write Eqs. (2.28)-(2.31) to obtain γ1 = ω2 R e−α|z1−z2| sin Φ (2.32) c1 = U1 − 1 α [ ω1 − ω2 R e−α|z1−z2| cos Φ ] (2.33) γ2 = Rω1e −α|z1−z2| sin Φ (2.34) c2 = U2 + 1 α [ ω2 −Rω1e−α|z1−z2| cos Φ ] (2.35) Eqs. (2.32)-(2.35) describe the linear hydrodynamic stability of the system. Unlike the conventional linear stability analysis, we did not impose normal- mode type perturbations (they only account for exponentially growing insta- bilities) in our derivation. Therefore the equation set provides a non-modal description of hydrodynamic stability in idealized (multi-layered) shear flows. We refer to this theory as the “wave interaction theory (WIT)”. WIT is only applicable to those hydrodynamic stability problems where the discrete spec- trum dynamics is of interest and the continuous spectrum can be neglected. A schematic description of the process of wave interaction is illustrated in Fig. 2.1. Subtracting Eq. (2.34) from Eq. (2.32) and Eq. (2.35) from Eq. (2.33), we 23 2.3. Interaction between two linear interfacial waves Figure 2.1: Schematic of the interfacial wave interaction mechanism. The deformation and associated vertical velocity of each wave is shown by the same color. Interaction imposes an additional vertical velocity (shown by different color). The horizontal arrow associated with a wave indicates the intrinsic wave propagation direction. Both the waves are counter-propagating (move against the background velocity at that location). find dR dt = R (γ1 − γ2) = ( ω2 −R2ω1 ) e−α|z1−z2| sin Φ (2.36) dΦ dt = α (c1 − c2) = α (U1 − U2)−[ ω1 +ω2 − ( Rω1 + ω2 R ) e−α|z1−z2| cos Φ ] (2.37) The two parameters have the following range of values: R ∈ (0, ∞) and Φ ∈ [−pi, pi]. Eqs. (2.36)-(2.37) represent a two dimensional, autonomous, non- linear dynamical system. The two equilibrium points of the system, found by imposing the conditions dR/dt = 0 in Eq. (2.36) and dΦ/dt = 0 in Eq. (2.37), 24 2.3. Interaction between two linear interfacial waves are given by (R, Φ) = (RNM ,ΦNM) and (RNM ,−ΦNM) (2.38) where RNM = √ ω2 ω1 (2.39) ΦNM = ± cos−1 [{ ω1 +ω2 − α (U1 − U2) 2 √ ω1ω2 } eα|z1−z2| ] (2.40) Eq. (2.40) reveals that the equilibrium points exist only if ∣∣∣∣{ω1 +ω2 − α (U1 − U2)2√ω1ω2 } eα|z1−z2| ∣∣∣∣ ≤ 1 (2.41) The linear behavior of the dynamical system around the equilibrium points is of interest. To understand this behavior, we evaluate the Jacobian matrix, J at the equilibrium points: J (RNM ,±ΦNM) = −2√ω1ω2e−α|z1−z2|  sin (±ΦNM) 0 0 sin (±ΦNM)  (2.42) Eq. (2.42) shows that the two eigenvalues corresponding to each equilibrium point are equal. Further analysis reveals that every vector at the equilibrium point is an eigenvector. The equilibrium point (RNM ,ΦNM) produces neg- ative eigenvalues, while the eigenvalues corresponding to (RNM ,−ΦNM) are 25 2.3. Interaction between two linear interfacial waves positive. Hence the dynamical system represented by Eqs. (2.36)-(2.37) is of “source-sink” type. In terms of the classical normal-mode analysis, each equilibrium point corresponds to a normal-mode of the discrete spectrum - (RNM ,ΦNM) corresponds to the growing normal-mode (signifying exponential growth) and (RNM ,−ΦNM) corresponds to the decaying normal-mode (signi- fying exponential decay). Normal-mode type instabilities can exist only if the condition in Eq. (2.41) is satisfied. Therefore Eq. (2.41) denotes the condition for exponentially growing instabilities in idealized, homogeneous and stratified shear layers. WIT allows understanding hydrodynamic instability from two different per- spectives - wave interaction and dynamical systems. According to the former, exponentially growing instabilities signify resonant interaction between the two waves. From dynamical systems point of view, resonance implies “equilibrium condition” (dΦ/dt = 0 and dR/dt = 0). Wave interaction interpretation of each of the two components of the equilibrium condition are as follows: (a) Phase-Locking : Reduction in the phase-speed of each wave occurs through the interaction mechanism - the vertical velocity field produced by the distant wave acts so as to diminish the phase-speed of the given wave. Furthermore, if the waves are “counter-propagating” (meaning, the direction of the intrin- sic phase-speed, −ωi sin (4φii) /α, is opposite to the background flow), the background flow causes an additional reduction in the phase-speed. Both wave interaction and counter-propagation work in tandem until the two waves 26 2.3. Interaction between two linear interfacial waves (a) (b) Figure 2.2: The figure shows that any initial condition (R0,Φ0) finally yields the resonant configuration (RNM ,ΦNM), provided Eq. (2.41) is satisfied. The case depicted here is Kelvin-Helmholtz instability (interaction between two vorticity waves) corresponding to α = 0.4. Any other shear instability will show qualitatively similar characteristics. (a) Φ versus t corresponding to Φ0 = −pi,−pi/2, 0,ΦNM , pi/2 and pi. The value of R0 is held constant, and is equal to 2. (b) R versus t corresponding to R0 = 1/2, 1(RNM), 2 and 5. The value of Φ0 is held constant, and is equal to −pi/2. 27 2.3. Interaction between two linear interfacial waves “phase-lock”, i.e. they are stationary relative to each other. In other words, this means dΦ/dt = 0. (b) Mutual Growth: The phase shift at the phase-locked state, ΦNM , is a unique angle producing the resonant configuration. This configuration causes the two waves to grow equally (i.e. γ1 = γ2) via interaction. Eq. (2.36) shows that equal growth rate implies dR/dt = 0. Furthermore, Eqs. (2.32) and (2.34) imply γ1 = γ2 = constant, meaning that the wave amplitudes grow at an exponential rate3. This exponential growth explains why the equilibrium point (RNM ,ΦNM) corresponds to the growing normal-mode of the discrete spectrum. WIT shows that the individual waves grow when 0 ≤ Φ ≤ pi and decay when −pi ≤ Φ ≤ 0. The largest possible growth, also known as the optimal growth, occurs when Φ = pi/2. These results are in accordance with the analysis of Heifetz & Methven (2005). Using the Generalized stability theory of Farrell & Ioannou (1996), Heifetz & Methven (2005) showed that the optimal growth in a barotropic shear layer occurs when the two counter-propagating Rossby waves are in quadrature. WIT also reveals that the left moving top wave and the right moving bottom wave eventually phase-lock (which then leads to mutual growth), pro- vided the condition in Eq. (2.41) is satisfied. Any arbitrary initial condition 3There are systems where phase-locking does not produce exponential growth. For ex- ample, stable barotropic and baroclinic modes result from the phase-locking between deep water surface gravity and internal gravity waves; see Chapter 7 of Kundu & Cohen (2004), Pg. 259-261. Shear is always absent in such systems. 28 2.3. Interaction between two linear interfacial waves (say R = R0 and Φ = Φ0) finally leads to phase-locking, and is evident from Fig. 2.2. As mentioned already, Eqs. (2.32)-(2.37) describe the non-modal in- stability process. Non-modal instability signifies non-orthogonal interaction between the two wave modes, and is the entire process occurring prior to the phase-locking event. Phase-locking is the final or steady state configuration, and corresponds to the the growing normal-mode (RNM ,ΦNM) of the discrete spectrum. The fact that ΦNM signifies the growing normal-mode configuration implies 0 ≤ ΦNM ≤ pi. A misconception might arise for phase shifts in the range −pi ≤ Φ ≤ 0. In this case, the reader might form an impression that the instability will not appear (because both the waves are decaying according to Eqs. (2.32) and (2.34)). However, the wave decay process is temporary. The two waves continuously adjust Φ so as to enter the growing zone 0 ≤ Φ ≤ pi. After reaching this zone, the two waves still continue to adjust Φ until the steady- state value (i.e. the resonant configuration), ΦNM , is reached. This fact can be better understood by considering the case Φ0 = −pi/2 in Fig. 2.2(a). Although initially −pi ≤ Φ ≤ 0, the value of Φ eventually enters the growing range and finally attains the steady-state value. 29 2.4. Homogeneous and stratified shear instabilities 2.4 Homogeneous and stratified shear instabilities 2.4.1 The Kelvin-Helmholtz instability Let us consider a piecewise linear velocity profile U (z) =  U1 z ≥ z1 Sz z2 ≤ z ≤ z1 U2 z ≤ z2 (2.43) This profile is a prototype of barotropic shear layers occurring in many geo- physical and astrophysical flows, see Chapter 3. It supports two vorticity waves, one at z1 and the other at z2. The shear S = (U1 − U2)/(z1 − z2). We non-dimensionalize the problem by choosing a length scale h = (z1 − z2)/2 and a velocity scale ∆U = (U1 −U2)/2. In a reference frame moving with the mean flow Ū = (U1 + U2)/2, the non-dimensional velocity profile becomes U (z) =  1 z ≥ 1 z −1 ≤ z ≤ 1 −1 z ≤ −1 (2.44) where both U and z are now non-dimensional quantities. This profile, along with the vorticity waves, is shown in Fig. 2.3(a). The top wave is left moving 30 2.4. Homogeneous and stratified shear instabilities (a) (b) Figure 2.3: (a) The setting leading to the Kelvin-Helmholtz instability. The velocity profile in Eq. (2.44) is shown on the left, while the vorticity waves (marked by “V”) are shown on the right. (b) Linear stability diagram of the Kelvin-Helmholtz instability (γ denotes the modal growth rate). 31 2.4. Homogeneous and stratified shear instabilities while the bottom wave is right moving. Both the waves counter-propagate, i.e. move in a direction opposite to the background flow. The wave interaction and subsequent instability mechanism can be understood from WIT, see §2.3. The classical normal-mode based linear stability analysis of the profile in Eq. (2.44) was first performed by Rayleigh (1880). He showed that if the non- dimensional wavenumber α is in the range 0 ≤ α ≤ 0.64, the flow is unstable; see Fig. 2.3(b). Thus, the piecewise linear profile and the ensuing instability are often referred to as the “Rayleigh’s shear layer” and “Rayleigh’s shear instability” respectively. However we will address the latter as the “Kelvin- Helmholtz instability (KH)”, following the wider acceptance of this terminol- ogy in the stratified shear layer community (Carpenter et al., 2013). The non-modal analysis of the piecewise linear profile was performed in detail by Heifetz et al. (1999); Heifetz & Methven (2005). Following the footsteps of Bretherton (1966) and Hoskins et al. (1985), Heifetz and co-authors were able to put forward a comprehensive mechanistic picture of KH (in rotating frame) in terms of counter-propagating Rossby wave interactions. By using the Gen- eralized Stability Theory (Farrell & Ioannou, 1996), Heifetz & Methven (2005) showed how wave interaction leads to optimal growth in shear layers. Here we study the KH problem in terms of WIT, i.e. Eqs. (2.32)-(2.35). Since the two waves involved in the KH problem are vorticity waves, we sub- stitute Eq. (2.14) in the WIT equation-set and after non-dimensionalization 32 2.4. Homogeneous and stratified shear instabilities we obtain γ1 = 1 2R e−2α sin Φ (2.45) c1 = 1− 1 2α [ 1− 1 R e−2α cos Φ ] (2.46) γ2 = R 2 e−2α sin Φ (2.47) c2 = −1 + 1 2α [ 1−Re−2α cos Φ] (2.48) Eqs. (2.45)-(2.48) are isomorphic to Eqs. (14a)-(14d) of Heifetz et al. (1999) and homomorphic to Eqs. (7a)-(7d) of Davies & Bishop (1994). These two referenced equation-sets describe edge wave interactions in two different types of rotating physical systems. While the one described by Heifetz et al. (1999) shows how CRW interactions lead to barotropic shear instability, the equation- set formulated by Davies & Bishop (1994) shows how baroclinic instability is produced through the interaction of temperature edge waves of the Eady model. Furthermore, Heifetz et al. (1999) showed that their set of equations is homomorphic to that of Davies & Bishop (1994). Eqs. (2.45)-(2.48) demonstrate how wave interaction causes amplitude growth and phase-speed modification of the individual vorticity waves, thereby lead- ing to KH. The fact that the wave interaction modifies the phase-speed of a vorticity wave can be understood by comparing Eq. (2.46) and Eq. (2.48) with the non-dimensional form of Eq. (2.16) (non-dimensionalization means substituting S = 1 and Ui = 1 or −1 in Eq. (2.16)). 33 2.4. Homogeneous and stratified shear instabilities The generalized non-linear dynamical system given by Eqs. (2.36)-(2.37) in this case translates to dR dt = 1 2 ( 1−R2) e−2α sin Φ (2.49) dΦ dt = (2α− 1) + 1 2 ( R + 1 R ) e−2α cos Φ (2.50) The equilibrium points of this system are (RNM ,±ΦNM), where RNM = 1 (2.51) ΦNM = cos −1 [(1− 2α) e2α] (2.52) The phase portrait is shown in Fig. 2.4. It confirms that the dynamical system is indeed of source-sink type, as predicted in §2.3. When the inter- action between the two waves is weak, the waves do not reach the resonant configuration, hence exponentially growing instability is never produced. The resulting dynamical system only produces periodic orbits; see Fig. 2.5. The necessary and sufficient condition for instability expressed via Eq. (2.41) in this case translates to −1 ≤ (1− 2α) e2α ≤ 1 implying 0 ≤ α ≤ 0.64 (2.53) The range of unstable wavenumbers obtained from the above equation corrob- orates Rayleigh’s normal-mode analysis. Rayleigh also found the wavenumber of maximum growth (also known as 34 2.4. Homogeneous and stratified shear instabilities Figure 2.4: Phase portrait of Kelvin-Helmholtz instability corresponding to α = 0.4. The system has two equilibrium points - one unstable (◦) and the other stable (•). Φ is the phase difference between the lower and upper waves, while R represents the ratio of the upper wave amplitude to the lower wave amplitude. the critical wavenumber) to be αcrit = 0.4. This value can be verified through WIT by imposing the normal-mode condition and maximizing γ1 or γ2 with respect to α. The fact that KH develops into a standing wave instability can be verified by applying the normal-mode condition in Eqs. (2.46) and (2.48). Performing the necessary steps we find c1 = c2 = 0, i.e. the waves have become stationary after phase-locking. In this configuration, the waves start to grow exponen- tially. Hence the shear layer grows in size. The growth process eventually becomes non-linear, and as shown in Chapter 3, the shear layer modifies into elliptical patches of constant vorticity. 35 2.4. Homogeneous and stratified shear instabilities Figure 2.5: Phase portrait when the two interacting vorticity waves do not produce KH (α = 2). 2.4.2 The Taylor instability Let us consider a uniform shear layer with two density interfaces U(z) = Sz and ρ (z) =  ρ0 − ∆ρ2 z ≥ z1 ρ0 z2 ≤ z ≤ z1 ρ0 + ∆ρ 2 z ≤ z2 (2.54) The shear S is constant. We choose ∆ρ/2 as the density scale, h = (z1− z2)/2 as the length scale, and thereby non-dimensionalize Eq. (2.54). The physical state of the system is determined by the competition between the density stratification and the shear, the non-dimensional measure of which is given by the Bulk Richardson number J = g ′ /(hS2), where g ′ = g(∆ρ/2)/ρ0 is the reduced gravity, and ρ0 is the reference density. The dimensionless velocity 36 2.4. Homogeneous and stratified shear instabilities (a) (b) Figure 2.6: (a) The setting leading to the Taylor instability. The velocity and density profiles in Eq. (2.55) are shown on the left, while the gravity waves “G” are shown on the right. (b) Linear stability diagram of the Taylor instability. The contours represent the growth rate. 37 2.4. Homogeneous and stratified shear instabilities and density profiles therefore become U(z) = z and ρ (z) =  −1 z ≥ 1 0 −1 ≤ z ≤ 1 1 z ≤ −1 (2.55) This flow configuration is shown in Fig. 2.6(a). Contrary to the conven- tional notion that gravitationally stable density stratified flows are usually stable, Taylor (1931) put forward the flow given by Eq. (2.55) and showed it to be linearly unstable. The interplay between the background shear and the gravity waves existing at the density interfaces produce the destabilizing ef- fect. The ensuing instability is thus known as the “Taylor instability”. Taylor found that for each value of J , there exists a band of unstable wavenumbers (and vice-versa), shown in Fig. 2.6(b). This unstable range is given by (see Eq. (2.154) of Sutherland (2010)) 2α 1 + e−2α ≤ J ≤ 2α 1− e−2α (2.56) Caulfield (1994), and more recently Carpenter et al. (2013), have described Taylor instability in terms of wave interactions. As discussed in §2.2.2, each density interface (located at z = 1 and z = −1) supports two gravity waves. The interaction between the left moving gravity wave at the upper interface and right moving gravity wave at the lower interface leads to Taylor instability. 38 2.4. Homogeneous and stratified shear instabilities To understand Taylor instability in terms of WIT, we substitute Eq. (2.20) in Eqs. (2.32)-(2.35). After non-dimensionalization and some algebra, we ob- tain γ1 = J 2R(1 + c2) e−2α sin Φ (2.57) c1 = 1− √ J 2α ( 1− β R e−2α cos Φ ) (2.58) γ2 = JR 2(1− c1)e −2α sin Φ (2.59) c2 = −1 + √ J 2α ( 1− R β e−2α cos Φ ) (2.60) Here β = ω2/ω1 = (1− c1)/(1 + c2), and by definition is a positive quantity. From Eqs. (2.58) and (2.60) we construct a quadratic equation for β: β2 + βe−2α cos Φ ( 1 R −R ) − 1 = 0 (2.61) Of the two roots, only the positive one is relevant. The coupled nature of Eqs. (2.58) and (2.60) makes it more complicated than the KH problem. The non-linear dynamical system in this case is given by dR dt = J 2 ( 1 1 + c2 − R 2 1− c1 ) e−2α sin Φ (2.62) dΦ dt = 2α− √ Jα 2 ( 1− β R e−2α cos Φ ) − √ Jα 2 ( 1− R β e−2α cos Φ ) (2.63) 39 2.4. Homogeneous and stratified shear instabilities At phase-locking R = RNM = √ β. Substituting this value in Eq. (2.61) gives β = 1. Therefore RNM = 1 and c1 = c2 = 0 at resonance. This implies that Taylor instability, like KH, also evolves into a standing wave instability. Although this fact is previously known, WIT demonstrates why this is the case. The non-linear structure of Taylor instability is similar to KH. Caulfield et al. (1995) has experimentally shown that this instability evolves into standing billows. The phase-shift ΦNM is evaluated from Eq. (2.63): ΦNM = cos −1 [( 1− 2α J ) e2α ] (2.64) The condition for Taylor instability is given by −1 ≤ ( 1− 2α J ) e2α ≤ 1 implying 2α 1 + e−2α ≤ J ≤ 2α 1− e−2α (2.65) The latter result corroborates the classical normal-mode result given in Eq. (2.56). 2.4.3 The Holmboe instability Let us consider the following velocity and density profiles U(z) =  U1 z ≥ z1 Sz z ≤ z1 and ρ (z) =  ρ0 z ≥ z2 ρ0 +4ρ z ≤ z2 (2.66) 40 2.4. Homogeneous and stratified shear instabilities Figure 2.7: Phase portrait of Taylor instability corresponding to an unstable combination of α and J . Here α = 0.2 and J = 0.7264. We non-dimensionalize Eq. (2.66) exactly like the Taylor problem, which gives us the dimensionless velocity and density profiles: U(z) =  1 z ≥ 1 z z ≤ 1 and ρ (z) =  0 z ≥ 0 2 z ≤ 0 (2.67) The vorticity interface at the top supports a vorticity wave, while the density interface at the bottom supports two gravity waves. The interaction between the left moving vorticity wave at the upper interface and the right moving gravity wave at the lower interface leads to an instability mechanism, known as the “Holmboe instability”. The corresponding flow setting is shown in Fig. 2.8(a). 41 2.4. Homogeneous and stratified shear instabilities (a) (b) Figure 2.8: (a) The setting leading to the Holmboe instability. The velocity and density profiles in Eq. (2.67) are shown on the left, while the vorticity wave “V” and the gravity wave “G” are shown on the right. (b) Linear stability diagram of the Holmboe instability. The contours represent the growth rate. 42 2.4. Homogeneous and stratified shear instabilities Holmboe (1962) was the first to consider the instability mechanism re- sulting from the interaction between vorticity and gravity waves. In his actual problem, Holmboe considered a flow setting more complicated than Eq. (2.67). His problem consisted of a velocity profile given by Eq. (2.44), however the density profile is the same as that in Eq. (2.67). Holmboe performed a linear stability analysis and showed that in addition to the conventional KH mode, there is another mode of instability - the Holmboe mode. Unlike the KH mode, the Holmboe mode is characterized by traveling waves. Presence of this unsta- ble mode reveals that stable density stratification can also have a destabilizing influence. This aspect of Holmboe instability is much like the Taylor instabil- ity. Recent non-modal analysis by Constantinou & Ioannou (2011) has shown that Holmboe instability is susceptible to substantial transient growths. Such growths especially occur for parameter values for which there is no instability but are close to the stability boundary. Analyzing the “authentic” Holmboe instability in terms of WIT implies considering the interaction of three waves - two vorticity waves and a gravity wave. An extended version of WIT can handle this problem, however this will not be considered in this work4. Baines & Mitsudera (1994) simplified Holm- boe’s problem by introducing the profile in Eq. (2.67). This allows studying the interaction of a vorticity and a gravity wave, and is therefore suitable for this work. Linear stability analysis shows that corresponding to each value of 4Stratified shear layer instabilities resulting from the interaction of multiple waves have been addressed by Caulfield (1994), however he limited his study to the normal-mode wave- form. 43 2.4. Homogeneous and stratified shear instabilities Figure 2.9: Phase portrait of Holmboe instability corresponding to α = 1 and J = 0.5. J , there exists a band of unstable wavenumbers. This is shown in Fig. 2.8(b). The stability boundary has been evaluated in Appendix 5.3.3. In order to understand Holmboe instability in terms of WIT, we substi- tute Eq. (2.14) and Eq. (2.20) in Eqs. (2.32)-(2.35). After performing non- dimensionalization and some algebra, we obtain γ1 = J Rc2 e−α sin Φ (2.68) c1 = 1− 1 α ( 1 2 − J Rc2 e−α cos Φ ) (2.69) γ2 = R 2 e−α sin Φ (2.70) c2 = 1 4α ( −Re−α cos Φ + √ R2e−2α cos2 Φ + 16αJ ) (2.71) 44 2.4. Homogeneous and stratified shear instabilities Like the Taylor case, this equation-set is also of coupled type. The non-linear dynamical system in this case is given by: dR dt = ( 4αJ −Re−α cos Φ +√R2e−2α cos2 Φ + 16αJ − R2 2 ) e−α sin Φ (2.72) dΦ dt = α− 1 2 ( 1−Re−α cos Φ)+ 4αJ −Re−α cos Φ +√R2e−2α cos2 Φ + 16αJ ( e−α cos Φ R − 1 ) (2.73) The equilibrium points of this system are (RNM ,±ΦNM), where RNM = √√√√1− 2α +√32αJ + (1− 2α)2 2 (2.74) ΦNM = cos −1 [( R2NM + 1− 2α 2RNM ) eα ] (2.75) The condition for Holmboe instability is found to be −1 ≤ ( R2NM + 1− 2α 2RNM ) eα ≤ 1 (2.76) This provides the range of J leading to Holmboe instability, and is as follows: 1 2A ( −B − √ B2 − 4AC ) ≤ J ≤ 1 2A ( −B + √ B2 − 4AC ) (2.77) 45 2.5. Conclusion where A = 16α2 B = −α [8 (2α− 1)2 + 36 (2α− 1) e−2α + 27e−4α] C = ( 2α− 1 + e−2α) (2α− 1)3 Eq. (2.77) corroborates the normal-mode result given in Appendix 5.3.3. The phase portrait of Holmboe instability, corresponding to an unstable combination of α and J , is shown in Fig. 2.9. This phase portrait is slightly different from Taylor and KH, because RNM 6= 1 in this case. Another feature of Holmboe instability is that, unlike Taylor and KH cases, its phase-speed is non-zero at the equilibrium condition. This phase-speed is found to be c1 = c2 = 2J R2NM = 4J 1− 2α + √ 32αJ + (1− 2α)2 (2.78) In the limit of large α and J , the two phase-locked waves move with unit speed to the right. 2.5 Conclusion This chapter is devoted to investigating the wave interaction problem in a generalized sense. The governing equations of hydrodynamic instability in ho- mogeneous and stratified shear layers have been derived without imposing the 46 2.5. Conclusion wave type, or the normal-mode waveform. We refer to this equation-set as the “wave interaction theory (WIT)”. Using WIT we have shown in Fig. 2.2 that two counter-propagating linear interfacial waves, having arbitrary initial amplitudes and phases, eventually resonates, provided they satisfy the con- dition for resonant interaction (Eq. (2.41)). Resonance makes the interfacial waves, and therefore the shear layer, to grow at an exponential rate - leading to the normal-mode instability. The condition for resonant interaction provides the range of unstable wavenumbers causing exponential growth. By consider- ing three different types of shear instabilities - Kelvin-Helmholtz, Taylor and Holmboe, we have shown that the resonant condition in each case matches the predictions of the canonical normal-mode based linear stability theory. In addition to perceiving the shear instability problem in terms of two interacting linear waves, WIT also provides an alternate perspective - under- standing shear instability in terms of dynamical systems. According to WIT, an unstable shear layer represents a non-linear dynamical system of the source- sink type, source and sink being the two equilibrium points. The equilibrium condition of the dynamical system is analogous to the resonant condition of the two-interacting-waves system. In terms of canonical linear stability the- ory, the source and the sink respectively correspond to the decaying and the growing normal-modes of the discrete spectrum. The most important aspect of WIT is that it provides a non-modal de- scription of idealized shear instabilities. Non-modal instability signifies non- orthogonal interaction between the two wave modes, and is the entire process 47 2.5. Conclusion occurring prior to resonance. Rapid transient growth is a key feature of non- modal instability process. WIT shows that optimal growth occurs when the two waves are in quadrature. The phenomenon of transient growth has been studied very briefly in this thesis; thorough research is required in future for its detailed understanding. 48 Chapter 3 Non-linear Kelvin-Helmholtz instability in a piecewise linear shear layer 1 The interplay of ideas and the oblique uses of knowledge are often of extraordinary interest. — Sherlock Holmes in The Valley of Fear. 3.1 Introduction Shear layers are ubiquitous in atmospheres and oceans. These layers can be- come hydrodynamically unstable, giving rise to Kelvin-Helmholtz instability (KH). The non-linear manifestation of KH is usually in the form of spiraling billows, whose breaking generates turbulence in geophysical flows. In theoretical and numerical studies, the hyperbolic tangent velocity profile is often used to model smooth barotropic shear layers (Hazel, 1972). Initially 1A version of this chapter has been published. A. Guha, M. Rahmani and G. A. Lawrence (2013) Evolution of a barotropic shear layer into elliptical vortices, Physical Review E. 49 3.1. Introduction interested in understanding the long time evolution of KH emanating from the hyperbolic tangent profile, we performed a direct numerical simulation (DNS); see Fig. 3.1(a). The flow re-laminarizes once the KH billow completely breaks down into small scales via turbulent processes. At this stage, the thickness of the shear layer has approximately quadrupled and, more importantly, the profile has almost become piecewise linear (see Fig. 3.1(b)). As shown in Fig. 3.1(c), the new base flow can be well described by Eq. (2.44), which serves as the initial condition for the subsequent instability processes. Detailed study on the linear instability arising from this base flow has been performed in Chapter 2.4. The study has shown that this velocity profile produces two counter- propagating vortcity waves (see Fig. 3.1(d)), which interacts to produce the subsequent Kelvin-Helmholtz instability mechanism. The underlying theory has been referred to as the wave interaction theory (WIT). The first non-linear analysis of a piecewise linear shear layer was performed by Pozrikidis & Higdon (1985). They used a boundary integral method known as contour dynamics to show that the piecewise linear shear profile evolves into nearly elliptical patches of constant vorticity - Kirchhoff vortices. We hypothesize that the initial shear layer profile determines the asymptotic form of the ensuing KH - hyperbolic tangent shear layers give rise to spiral billows, while piecewise linear shear layers produce Kirchhoff vortices. The spiraling billow form of KH has been thoroughly investigated in the past. In fact, the spiral billow shape has become the signature of KH (Smyth & Moum, 2012); but little is known about the non-linear evolution of the 50 3.2. Non-linear simulation piecewise linear shear layer. This is because the piecewise linear profile is usu- ally considered to be of little practical relevance, hence its usage is restricted to theoretical studies - mainly as an approximation of smooth shear layers (Pozrikidis, 1997; Carpenter et al., 2013). On the contrary, our DNS result in Fig. 3.1(a,b) indicates that the quasi piecewise linear profile is also likely to occur in nature. The fact that this profile produces elliptical vortices similar to those observed in geophysical and astrophysical flows, e.g. meddies in At- lantic ocean (Chérubin et al., 2003), stratospheric polar vortices (Waugh & Randel, 1999), Great Red Spot and other Jovian vortices (Morales-Jubeŕıas et al., 2002), Neptune’s Great Dark Spot (Polivani et al., 1990), and coherent vortices in the atmosphere of Uranus (Liu & Schneider, 2010), has motivated us to investigate further. 3.2 Non-linear simulation Over the last few decades, sophisticated computational techniques have been developed for precise understanding of turbulent processes. Two such tech- niques worth mentioning are direct numerical simulation (DNS) and vortex methods. Many geophysical and astrophysical flows can be assumed homogeneous, incompressible and quasi-inviscid. In such flows, vorticity plays a major role in driving non-linear processes like chaos and turbulence (Saffman, 1995). Vortex methods are especially useful under such circumstances; they numerically solve 51 3.2. Non-linear simulation Figure 3.1: (a) 3D DNS performed to capture the complete turbulent dis- sipation of a KH billow ensuing from a hyperbolic tangent velocity profile. False color is added to aid visualization. (b) Spanwise averaged mean velocity profile corresponding to each instant. (c) The magenta line is the continuous velocity profile obtained from Eq. (3.7), while the thick black line below it is the piecewise linear profile from Eq. (2.44). (d) CVWs (exaggerated) existing at the vorticity discontinuities. 52 3.2. Non-linear simulation the inviscid and incompressible Navier-Stokes equations (Euler equations). An example of one such 2D vortex method is contour dynamics (Deem & Zabusky, 1978). 3.2.1 Contour dynamics High Reynolds number flows have a tendency to develop finite-area vortex regions or “vortex patches” with steep sides (Deem & Zabusky, 1978; Pullin, 1992). If the flow is incompressible and 2D, the area as well as the vorticity of a vortex patch are conserved quantities. The vortex patch boundary is referred to as a “contour”, and is both a vorticity jump and a material surface (Saffman, 1995). A small perturbation on this contour sets up a vorticity wave (Deem & Zabusky, 1978). Substantial simplification is possible for constant vorticity patches; the governing 2D Euler equations can be reduced to a 1D boundary integral. This provides a significant advantage for computing the vortex patch evolution by only solving the contour motion. This methodology is known as the contour dynamics (CD) (Pullin, 1992). The contour motion is obtained by tracing the patch boundary with N Lagrangian markers. The evolution of the i- th Lagrangian marker is obtained by solving the following integro-differential equation: dxi dt = − Ω 4pi ∫ C ln [( 4x′i )2 + ( 4y′i )2] dx ′ (3.1) 53 3.2. Non-linear simulation Figure 3.2: Contour dynamics simulation of a piecewise linear shear layer by Pozrikidis & Higdon (1985) for (a) α = 0.639, (b) α = 0.5 and (c) α = 0.0625. where Ω = 1 is the patch vorticity, x = [x, y]T , 4x′i = xi − x′ , 4y′i = yi − y′ and C is the contour around one patch5. Observing that a piecewise linear shear layer can be represented by a hori- zontally periodic patch of constant vorticity, Pozrikidis & Higdon (1985) mod- ified Eq. (3.1) for a periodic domain (Pozrikidis, 1997) : dxi dt = − Ω 4pi ∫ C ln [ cosh ( α4y′i ) − cos ( α4x′i )] dx ′ (3.2) On perturbing the contour with sinusoidal disturbances, the shear layer rolls up producing nearly elliptical patches of constant vorticity (Pozrikidis & Higdon, 1985). 5Unlike Chapter 2, the cross streamwise coordinate is denoted by ‘y’ instead of ‘z’. This is primarily to stress that there is no density variation along this coordinate. 54 3.2. Non-linear simulation We interpret the rolling up of the piecewise linear shear layer in terms of the interaction between two CVWs. To better understand this point, let us first consider the simple case of a circular patch of constant vorticity. A small perturbation on its contour produces a linearly stable vorticity wave (also known as “Kelvin wave”) (Deem & Zabusky, 1978; Mitchell & Rossi, 2008). Similarly, perturbation on a periodic, piecewise linear shear layer produces two counter-propagating vorticity waves, one at the upper interface and the other at the lower interface. These two CVWs interact with each other causing the shear layer to evolve into elliptical vortex patches. Since WIT represents the linear dynamics of the shear layer evolution process, CD can be regarded as a non-linear extension of WIT. We assume the pair of CVWs to be initially of the form: η+ = a0 cos (αx− Φ) at y = 1 (3.3) η− = a0 cos (αx) at y = −1 (3.4) where the perturbation amplitude a0  1 and the phase shift Φ ∈ [−pi, pi]. Eqs. (3.3)-(3.4) are similar to the initial perturbation profiles assumed by Pozrikidis & Higdon (1985), except that they interpreted these equations as contour perturbations. Also our choice of α and Φ are different from Pozrikidis & Higdon (1985). We set α = αcrit = 0.4 (i.e. the fastest growing mode) to ensure that the resulting non-linear structure is most likely to be realized in nature; whereas, Pozrikidis and Higdon used several values of α, not including 55 3.2. Non-linear simulation Figure 3.3: Time evolution of Kelvin-Helmholtz instability - comparison be- tween DNS and CD. 56 3.2. Non-linear simulation αcrit. Corresponding to each α, Pozrikidis & Higdon (1985) observed the shear layer to produce a unique non-linear structure; see Fig. 3.2. To ensure exponential growth from time t = 0, we set Φ = θcritNM = 0.353pi; whereas, in Pozrikidis and Higdon, Φ = 0 or pi. The angle θcritNM is θNM corresponding to αcrit. Recall from Chapter 3 that θNM is the phase-locking angle (see Eq. (2.52)), and is given by: θNM = cos −1 [(1− 2α) e2α] In summary, WIT provides the initial conditions (i.e. αcrit and θ crit NM) for mod- eling the non-linear evolution. We follow the numerical method outlined in Pozrikidis (1997) to solve Eq. (3.2). The contour is approximated by a polygonal line which connects the successive marker points. The integral over a segment that does not contain the i-th marker point is easy to compute by using standard numerical inte- gration methods, e.g. trapezoidal rule. However the integral over the adjacent segments Sj of the i-th marker point, corresponding to j = i−1 and i, exhibits logarithmic singularity and is computed as follows: ∫ Sj ln [ cosh ( α4y′i ) − cos ( α4x′i )] dx ′ = ∫ Sj ln [ cosh ( α4y′i )− cos (α4x′i)(4x′i)2 + (4y′i)2 ] dx ′ + ∫ Sj ln [( 4x′i )2 + ( 4y′i )2] dx ′ (3.5) 57 3.2. Non-linear simulation The first integral on the R.H.S. is non-singular and is computed by using trapezoidal rule. The singularity has been shifted to the second integral, which can fortunately be integrated analytically: ∫ Sj ln [( 4x′i )2 + ( 4y′i )2] dx ′ = (xj+1 − xj) ( ln |xj+1 − xj|2 − 2 ) (3.6) We consider a domain one wavelength (λcrit ≡ 2pi/αcrit = 5pi) long, and solve Eq. (3.2) using central differencing for space derivatives and 4th order Runge-Kutta for time. Each wave is initially represented by 400 points. During its evolution process, an adaptive point insertion-deletion algorithm is used to check if the neighboring points are within a desired distance. 3.2.2 Direct numerical simulation Previous studies have used CD mainly as a tool for qualitative understand- ing of problems involving inviscid vortical flows (Pullin, 1992). In order to demonstrate the quantitative capabilities of CD, we validate our CD simula- tion against a pseudo-spectral DNS. The DNS code uses full Fourier transforms in the horizontal direction, and half-range sine or cosine Fourier transforms in the vertical direction in order to convert the set of partial differential equa- tions (Navier-Stokes equations) into ordinary differential equations. Time in- tegration is performed using a third-order Adams-Bashforth method. Detailed description of this code can be found in Winters et al. (Winters et al., 2004). We consider a domain of length λcrit in the horizontal direction and nine 58 3.2. Non-linear simulation times the initial shear layer thickness in the vertical direction. The horizontal boundary condition is periodic while the vertical boundary condition is no- flux free-slip. We perform a 2D simulation at Reynolds number Re = 10, 000 (Re = 1/ν, where ν is the fluid viscosity). This Reynolds number is high enough to mimic quasi-inviscid flow conditions. To simulate the smallest scales of motion, we resolve our domain using 2880 points in horizontal and 3456 points in vertical direction. Since non-differentiable profiles like Eq. (2.44) are subject to Gibbs phe- nomena, we use a smooth velocity profile that resembles the piecewise linear profile very closely; see Fig. 3.1(c). This velocity profile is derived from a vorticity distribution having the form Ω (y) = 1 2 [ 1− tanh ( y2 − 1  )] (3.7) Integrating Eq. (3.7), the velocity profile is obtained directly: U = ∫ Ωdy. The linear stability characteristic of this profile matches with that of the piecewise linear profile almost exactly, and has the same αcrit. Equating the total cir- culation of DNS with CD yields  = 0.100. The vorticity field in DNS is perturbed to match the initial wave amplitude growth in CD. 59 3.3. Results and discussion Figure 3.4: Temporal variation of the wave amplitude a. The straight blue line is the prediction from linear theory, the black line corresponds to CD while the magenta line represents DNS results. The green and red lines in the inset respectively show the variation of the ellipse aspect ratio r and the angular rotation rate ω with time. These variations are obtained by solving Eqs. (3.9)-(3.10). The black markers indicate the corresponding data points measured from DNS. 3.3 Results and discussion 3.3.1 Pre-saturation and saturation phases The non-linear evolution of KH is illustrated in Fig. 3.3. It shows the vorticity field from DNS and contour lines from CD. The vorticity of the region enclosed by the contour lines is conserved in CD. However the presence of viscosity makes conservation of vorticity invalid in DNS. The implementation of high Re DNS minimizes the viscous effects, making DNS comparable to CD. The basic premise behind our simulations is the conservation of total circulation Γ = ΩA (where the vorticity Ω = 1 and A is the shear layer area), which 60 3.3. Results and discussion comes from Kelvin’s Circulation theorem (Saffman, 1995). Its corollary is the conservation of shear layer area - a quantity that remains fixed at its initial value A = 2λcrit. Fig. 3.4 shows the time evolution of the wave amplitude a = |max(η+) − 1| = |min(η−) + 1|. The maximum shear layer thickness, H = 2(1 + a) evolves in a fashion similar to the wave amplitude a. We find the growth to be exponential, at least for t ≤ 20. CVW interaction causes the shear layer to grow non-linearly. This phenomenon leads to the roll-up and formation of the elliptical core vortex (Kirchhoff vortex). The evolution process is shown in Fig. 3.3. The part of the shear layer between the crest of the lower wave and the trough of the upper wave (see Fig. 3.1(d)) gives rise to the elliptical core, the initial length of which is given by linit = ( 1 + θcritNM pi ) λcrit 2 (3.8) The flow saturates (i.e. the amplitude reaches a maxima) at tsat = 34. For t ≥ tsat, approximately 80% of Γ is concentrated in the core. H also reaches a maxima at saturation, and has the value Hmax = 8.7. The fully formed elliptical cores are connected by thin filaments of fluid called braids. These braids wrap around the rotating cores, producing a complex spiral like structure, as shown in Fig. 3.3. 61 3.3. Results and discussion 3.3.2 Early post-saturation phase After saturation, the core rotates with an angular velocity ω, causing the wave amplitude, a, to oscillate with a time period Tamp ≈ 13; see Fig. 3.4. The core also nutates, i.e. the core aspect ratio r (defined as the ratio between the ellipse major axis and the minor axis) undergoes a periodic oscillation. This nutation phenomenon is apparent in both Figs. 3.3 and 3.4. Nutation To better understand the nutation process, we consider the simple model pro- posed by Kida (1981). An isolated Kirchhoff vortex rotates in the presence of a constant background strain-rate γ. The velocity field associated with this strain-rate is given by us = γσ, ws = −γξ where σ and ξ are the principal axes with the origin at the centre of the ellipse. In our case, this velocity field mimics the leading order straining effect induced by the rotation of other Kirchhoff vortices. Note that the periodic boundary condition takes into ac- count the effects of other Kirchhoff vortices. Let the clockwise angle between σ and the ellipse major axis be θ at any instant. Then θ and r evolve as follows (Kida, 1981): ω ≡ dθ dt = −γ ( r2 + 1 r2 − 1 ) sin (2θ) + Ωr (r + 1)2 (3.9) dr dt = 2γr cos (2θ) (3.10) Eq. (3.10) implies that the nutation is caused by strain. It also reveals that r 62 3.3. Results and discussion reaches maxima at θ = ±pi/4. Simultaneously, Fig. 3.3 shows that the core nutates with a maximum value of r along the x axis and a minimum along the y axis. Therefore the σ axis must make an angle of pi/4 with the x axis. The angle made by the braid with the x axis at the stagnation point(s) is also pi/4. This is because the braid aligns itself with the streamlines. To investigate why the σ axis makes an angle of pi/4 with the x axis, we consider an ideal problem where an infinite number of Kirchhoff vortices, each of circulation Γcore = Γ (note Γ = 2λcrit), are placed along the x axis with a constant spacing λcrit between their centers. The effect of a vortex patch at a point outside its own boundary is exactly the same as that of a point vortex of equivalent strength placed at the center of the patch (Pozrikidis, 1997). Therefore we replace all the Kirchhoff vortices with point vortices of strength Γcore. This provides a simplistic understanding of the mechanism by which the rotation of distant vortex patches strain a given patch. We find this ideal strain-rate to be γ′ = Γcore 2piλ2crit ∞∑ n=−∞,n6=0 n−2 = 0.067 (3.11) The principal axes of the strain field produced by this infinite array of point vortices make angles of ±pi/4 with the x axis. Hence, this ideal strain field and the strain field of our actual problem have the same orientation. Before comparing the magnitudes of these two fields, it is important to note that the presence of braids complicate the actual problem by making the strain- rate magnitude vary spatially. We simplify the analysis by assuming a strain 63 3.3. Results and discussion field of constant magnitude acting on the elliptical core, thereby reducing the problem to the Kida problem described by Eqs. (3.9)-(3.10). DNS is used to supplement the analysis by providing the values of r and θ wherever necessary. By applying this methodology, the magnitude of the actual strain-rate is found to be γ = 0.073, which is quite close to the ideal value of 0.067 obtained from Eq. (3.11). We also capture the evolution of r and ω by solving Eqs. (3.9)-(3.10); see the inset in Fig. 3.4. The initial values are obtained from DNS, and γ = 0.073. The spike in ω at t ≈ 45 is caused by r → 1. The figure shows that Kida’s model compares well with the DNS. In DNS we fit an ellipse on the Kirchhoff vortex using a least-square based algorithm (Hendel, 2008). Points on the surface of the vortex are selected manually and then the best-fit ellipse is generated. The code supplies the length of the ellipse axes, as well as the orientation. The nutation period is found to be Tnut ≈ 13, while the period of core rotation is Tcore = 2pi/ω̄ ≈ 26 (overbar denotes average). Tcore ≈ 2Tnut is because one full rotation corresponds to passing the coordinate axes twice. Likewise, Tcore ≈ 2Tamp because the braids are connected to the two ends of the core. Small length scale production The smallest length scales are found to occur in the braid region adjacent to the core; see Fig. 3.3. This is due to the straining effect of the rotating elliptical 64 3.3. Results and discussion t = 67 Figure 3.5: Formation of winding filaments around the elliptical vortex during the late post-saturation phase. core which causes the braid region in its vicinity to thin exponentially fast. In real flows, when a fluid element becomes sufficiently thin, the balance between the strain-rate and the viscous dissipation determines the small length scales. The order of magnitude of the core rotation induced strain-rate, γlocal, can be obtained by replacing the core with a point vortex of equivalent strength and located at the ellipse centre: γlocal ∼ Γcore 2pil2 = 1 2 (3.12) where l is the characteristic length of the core. Notice that this local strain- rate is one order of magnitude greater than the background strain-rate γ or γ′. 65 3.4. Conclusion The smallest length scale appearing in a 2D “turbulent” flow is L2D ∼ Re−1/2 (Davidson, 2004). This length scale is a 2D analogue of the Taylor microscale occurring in 3D turbulent flows. In order to estimate the time when L2D appears in our flow, we formulate a braid evolution equation similar to Eq. (2.8) of Corcos and Sherman (Corcos & Sherman, 1976): δ2 (t∗) = δ2 (0) e−2γlocalt ∗ + pi 2γlocalRe ( 1− e−2γlocalt∗) (3.13) where δ (t∗) is the braid thickness adjacent to the core at time t∗ = t − tsat. We estimate δ (0) from our DNS and solve Eq. (3.13). L2D is found to appear soon after saturation, around t∗ ≈ 4, implying that 2D transitional flows like KH can give rise to “turbulent” features at a very early stage. 3.4 Conclusion When a piecewise linear shear layer becomes unstable, it evolves into a series of elliptical vortices of constant vorticity (Kirchhoff vortices) connected by thin braids. The interaction between two counter-propagating vorticity waves is the driving mechanism behind this instability process. Although this wave interaction perspective was known previously, linearized approximations forced the analysis to be valid only in the linear regime. By finding and exploiting the link between two quite different theories, namely wave interaction theory and contour dynamics, we are able to extend the analysis to the fully non-linear 66 3.4. Conclusion domain. The production of Kirchhoff vortices shows that KH arising from a piece- wise linear shear layer is very different from the classical spiraling billow type KH ensuing from a smooth shear layer. The characteristics of this little known KH have been investigated. The rotation and nutation of the Kirchhoff vor- tices are found to be consistent with the predictions of Kida. The time period of rotation of these vortices is twice the period of nutation and the period of maximum shear layer height oscillation. The braids connecting the Kirchhoff vortices thin exponentially fast to a length scale which is the 2D equivalent of Taylor microscale. Elliptical vortical structures, similar to those found in our simulations, are quite common in nature, especially in regions with quasi piecewise linear shear. Examples of such vortices include meddies in the Atlantic ocean, stratospheric polar vortices, and vortices in the gas giant planets like Jupiter, Neptune and Uranus. Our analysis may motivate further investigation of their formation and evolution. 67 Chapter 4 Linear stability of a Plane Couette-Poiseuille flow in the presence of cross-flow 1 Each fact is suggestive in itself. Together they have a cumulative force. — Sherlock Holmes in The Bruce-Partington Plans. 4.1 Introduction From the perspective of applications in technology, Poiseuille flow of viscous fluid along a duct is undoubtedly one of the most important flows studied as it underpins the field of hydraulics. Instability and subsequent transition from laminar flow marks a paradigm shift in the dominant transport mechanisms of mass, momentum and heat, and it is for this reason that the subject remains of enduring interest, even after more than 100 years of study. In this chapter we focus on two methods for affecting the linear stability of plane Poiseuille 1A version of this chapter has been published. A. Guha and I.A. Frigaard (2010) On the stability of plane Couette-Poiseuille flow with uniform crossflow, Journal of Fluid Mechanics. 68 4.1. Introduction (PP) flow. The first method consists of introducing a Couette component to the flow, by translation of one of the walls. The second method consists of introducing a cross-flow, e.g. via injection through a porous wall. While both effects have been studied individually to some extent, there are fewer studies of the two effects combined, which is the main focus here. PP flow is linearly unstable when the critical Reynolds number exceeds Rec ≈ 5772, (Reynolds number based on the center-line axial velocity and the half-width of the channel; see Orszag (1971)), whereas the plane-Couette (PC) flow is absolutely stable with respect to infinitesimal amplitude disturbances, Rec = ∞; see Romanov (1973). Superimposing PP and PC flows, we may ask if a small Couette component can affect the stability of the PP flow. Stability of plane Couette-Poiseuille (PCP) flow was first studied by Potter (1966) and later by Hains (1967), Reynolds & Potter (1967) and Cowley & Smith (1985). The results are typically understood with respect to a Reynolds number that is based on the maximal velocity of the Poiseuille component, say Rep, and the ratio of wall velocity to maximal velocity of the Poiseuille component, denoted k. For small Couette components, k, it is possible to observe some destabilization of the flow, (depending on the wavenumber), but as soon as k > 0.3 a strong stabilization of the flow sets in. As the velocity ratio k exceeds 0.7, the neutral stability curve completely vanishes and the flow becomes unconditionally linearly stable, i.e. Rec → ∞. The term “cut- off” velocity has been used to describe this stabilization; see Reynolds & Potter (1967). 69 4.1. Introduction The stability of PP flow with cross-flow was first analysed by Hains (1971) and Sheppard (1972), both of whom have shown that a modest amount of cross-flow produces significant increase of the critical Reynolds number. These results are however slightly problematic to interpret in absolute terms, since at a fixed pressure gradient along the channel, increasing the cross-flow decreases the velocity along the channel, (hence effectively the Reynolds number). This difficulty was noted by Fransson & Alfredsson (2003), who used the maxi- mal channel velocity as their velocity scale (instead of that based on the PP flow without cross-flow), and thus separated the effects of base velocity mag- nitude from those of the base velocity distribution. Using this velocity scale in their Reynolds number Re, they showed regimes of both stabilization and destabilization as the cross-flow Reynolds number was increased. For example, for Re = 6000 and wavenumber α = 1, Fransson & Alfredsson (2003) have shown that the cross-flow was stabilizing up to a cross-flow Reynolds num- ber Reinj ≈ 3.4, and then starts destabilizing before re-stabilizing again at Reinj ≈ 635. The initial regime of stabilization is the one corresponding to the earlier results. In the present study we focus on the combination of cross-flow and Cou- ette component. Our motivation stems from a desire to understand how the two mechanisms interact, since in terms of technological application different mechanical configurations may be more or less amenable to cross-flow and/or wall motion. This means that there is value in knowing when one effect may compensate for the other in stabilizing (or destabilizing) a given flow. To the 70 4.2. Stability of plane Couette-Poiseuille flow with cross-flow best of our knowledge, stability of PCP flow with cross-flow has only been studied in any generality by Hains (1971). In considering the base flow for PCP flow with cross-flow, (which is parameterized by Reinj and k), the rela- tion kReinj = 4, defines an interesting paradigm in which the base velocity in the axial direction is linear. These Couette-like flows have been studied by Nicoud & Angilella (1997) for increasing Reinj. They found a critical value of Reinj ≈ 24, below which which no instability occurs, (we have translated their critical value of 48 into the Reinj that we use). Therefore, we observe that understanding of cross-flow PCP flows is far from complete. We aim to contribute to this understanding. The 3D linear stability of PCP flow with cross-flow is amenable to Squires transformation, so that the linear instability occurs first for 2D (spanwise- independent) perturbations. It is these perturbations that we study here. Our aim is to demarcate clearly in the (Reinj, k)-plane, regions of unconditional stability, i.e. where there is a cut-off wall velocity or injection velocity. We also wish to understand the underlying linear stability mechanisms as Reinj and k are varied. 4.2 Stability of plane Couette-Poiseuille flow with cross-flow The base flow considered in this work is a plane Couette-Poiseuille flow (PCP) with imposed uniform cross-flow. This flow is two-dimensional, viscous, incom- 71 4.2. Stability of plane Couette-Poiseuille flow with cross-flow pressible and fully developed in the streamwise direction, x. The imposed base velocity in the y-direction 6, v , is constant and equal to the injection/suction velocity Vinj. Since v is constant, the x-component of velocity, u depends only on y. The flow domain is bounded by walls at y = ±h, and is driven in the x-direction by a constant pressure gradient and by translation of the upper wall, at speed Uc. The x-component of velocity, u(y), is found from the x-momentum equation, which simplifies to: Vinj ∂u ∂y = −1 ρ ∂p ∂x + ν ∂2u ∂y2 (4.1) where ρ is the density, ν = µ/ρ is the kinematic viscosity, and µ is the dynamic viscosity. The boundary conditions at y = ±h are: u(−h) = 0, u(h) = Uc (4.2) To scale the problem we scale all lengths with h, hence (x, y) = (x/h, y/h) (i.e. henceforth x and y are non-dimensional variables) For the velocity scale two choices are common. First, the imposed pressure gradient defines a “Poiseuille” velocity scale: Up = −h 2 2µ ∂p ∂x , (4.3) which is equivalent to the maximum velocity of the plane Poiseuille flow, driven by the pressure gradient alone. Second, we may take the maximum velocity, 6Unlike Chapter 2, the cross streamwise coordinate is denoted by ‘y’ instead of ‘z’. This is primarily to stress that there is no density variation along this coordinate. 72 4.2. Stability of plane Couette-Poiseuille flow with cross-flow which we need to compute. Up is the choice of Potter (1966), and thus allows one to compare directly with the studies of PCP flows. In the absence of a cross-flow, the maximal velocity is not actually very sensitive to the wall velocity Uc, at least for Uc < Up, which covers the range over which the flow stabilizes. However, in the case of a strong cross-flow the x-velocity is reduced significantly below Up, which therefore loses its meaning. Consequently we adopt the second choice and scale with the maximal velocity, Umax. This choice retains physical meaning in the base velocity, but does introduce algebraic complexity. The solution is found from (4.1)–(4.2) after detailed but straightforward algebra: u(y) Up = [ 4 coshReinj − kReinje−Reinj + [kReinj − 4]eReinjy 2Reinj sinhReinj + 2y Reinj ] , (4.4) where k and Reinj are defined by: k = Uc Up , (4.5) Reinj = Vinjh ν . (4.6) These two dimensionless parameters uniquely define the dimensionless base flow. The parameter k is the velocity ratio of Couette to Poiseuille velocities, which is useful as it allows direct comparison with earlier results on stabi- lization of PCP flows without cross-flow. The parameter Reinj is simply a 73 4.2. Stability of plane Couette-Poiseuille flow with cross-flow Reynolds number based on the injection velocity. Primarily here we consider the ranges: k ∈ [0, 1] and Reinj ≥ 0. For relatively weak crossflow velocities, the velocity component u(y) has a single maximum at a value of y = ymax defined by: eReinjymax = sinhReinj Reinj 4 4− kReinj (4.7) The maximal velocity Umax, is then evaluated from (4.4). Since, sinhReinj ≥ Reinj, we can see that ymax > 0 for k ≥ 0 and Reinj > 0. Both the injection cross-flow and Couette component act to skew the velocity profile towards the upper wall. For stronger cross-flow velocities, (or sufficiently large k), the maximal velocity occurs at the upper wall, i.e. Umax = Uc. The division between weak and strong cross-flows, taking into account also the Couette component, is defined by the line kReinj = 4 [ 1− sinhReinj ReinjeReinj ] (4.8) 74 4.2. Stability of plane Couette-Poiseuille flow with cross-flow The dimensionless base velocity is given by: u(y) =  4 coshReinj − kReinje−Reinj − [4− kReinj]eReinjy + 4y sinhReinj 4 coshReinj − kReinje−Reinj − [4− kReinj]eReinjymax + 4ymax sinhReinj , kReinj ≤ 4 [ 1− sinhReinj Reinje Reinj ] , 1 k [ 4 coshReinj − kReinje−Reinj − [4− kReinj]eReinjy 2Reinj sinhReinj + 2y Reinj ] , kReinj > 4 [ 1− sinhReinj Reinje Reinj ] . (4.9) It can be verified that in the limit Reinj → 0, with k fixed, the classical form of PCP base velocity profile is retrieved: u(y) ∼ 1− y 2 + k 2 (1 + y) +Reinj[ 1 3 (y − y3)− k 4 (1− y2)] 1 + k 2 + k 2 16 −Reinj[k6 − k 3 64 ] (4.10) as Reinj → 0, with k ≤ 4[1− sinhReinj/(ReinjeReinj)]/Reinj ∼ 4[1 +Reinj/3]. Examples of the base velocity profile are given in Fig. 4.1, for for k = 0.5 and different values of Reinj. Observe that for Reinj = 8, when kReinj = 4, the velocity profile is linear. This flow has been termed a “generalised Couette” flow by Nicoud & Angilella (1997). We shall denote differentiation with respect to y by the operator D. The first and second derivatives of the base flow, Du andD2u respectively, influence the stability of the flow. We find that D2u has sign determined by (kReinj−4), and increases in absolute value exponentially towards the upper wall. For kReinj < 4, the velocity is concave, and is convex otherwise. Since D 2u does not change sign, the maximal absolute value of the first derivative is found 75 4.2. Stability of plane Couette-Poiseuille flow with cross-flow Figure 4.1: Mean Velocity Distribution for k = 0.5 and Reinj = 0, 1, 4, 8, 15 and 30 (Reinj = 30 marked with a ) either at the upper or lower wall, y = ±1. The maximal velocity gradients are found at the lower wall for small Reinj and also for a range of Reinj close to kReinj = 4, but otherwise are found at y = 1; see Fig. 4.2(a). At large Reinj the maximal velocity increases almost linearly: |Du|max = |Du(y = 1)| = 1 k [ [kReinj − 4]eReinj 2 sinhReinj + 2 Reinj ] ∼ Reinj − 4 k + 2 kReinj +O(e−2Reinj). (4.11) Figure 4.2(b) shows examples of the profiles of D2u. We observe that 76 4.2. Stability of plane Couette-Poiseuille flow with cross-flow (a) (b) Figure 4.2: (a) Maximal velocity gradient, |Du|max, plotted against Reinj for k = 0.35, 0.5, 0.65, (k = 0.65 marked with a ). The thick line indicates where the maximum is attained at y = −1; otherwise at y = 1. (b) Variation of D2u with y for k = 0.5 for : Reinj = 0, (); Reinj = 4, (◦); Reinj = 8, (×); Reinj = 12, (). D2u ≈ 0 over a large range of y, close to the lower wall, whenever a significant amount of cross-flow is present, i.e. Reinj & 1. 4.2.1 Streamwise Reynolds number The base base flow is fully defined by the parameters k and Reinj, as discussed above. In addition, the transient flow and associated stability problem will depend on the streamwise Reynolds number, Re, which we define in terms of Umax, i.e. Re = Umaxh ν . (4.12) 77 4.2. Stability of plane Couette-Poiseuille flow with cross-flow 4.2.2 The stability problem The base flow is two-dimensional, but since v = Reinj is constant, the 3D linear stability equations are only modified by the addition of a constant convective term: Reinj ∂ ∂y u′, where u′ = (u′, v′, w′) denotes the linear perturbation. The classical Squire transformation can therefore be applied to the temporal problem, showing that for any unstable 3D linear disturbance there exists an unstable 2D linear disturbance at lower Re; see Squire (1933). It suffices to consider only 2D disturbances and we adopt the normal mode approach introduced in Chapter 1. However the inclusion of cross-flow term adds an extra inertial component ReinjD(α 2−D2)φ (where φ is the amplitude of the perturbation streamfunction) to the conventional Orr-Sommerfeld (O-S) equation (see Eq. 1.16): iαRe[(c− u)(α2 −D2)−D2u]φ−ReinjD(α2 −D2)φ = (α2 −D2)2φ, (4.13) The boundary conditions in this case are given by φ(±1) = Dφ(±1) = 0. (4.14) Note that Reinj in Eq. (4.13) also influences stability via the base velocity profile u(y). Equations (4.13)–(4.14) constitute the eigenvalue problem. The 78 4.2. Stability of plane Couette-Poiseuille flow with cross-flow eigenvalue c is parameterised by the 4 dimensionless groups (α,Re,Reinj, k) and the condition of marginal stability is: ci(α,Re,Reinj, k) = 0 (4.15) We attempt to characterise the stability of (4.13)–(4.14) for positive (α, Re, Reinj, k). We may note that the limit Re → ∞ for finite Reinj, reduces (4.13) to the Rayleigh equation. Since D2u is of one sign only, there are no inflection points and hence no purely inviscid instability. This suggests that the instabilities of (4.13)–(4.14) will be viscous in nature. Addition of the constant cross-flow terms does not fundamentally alter the O-S problem, and we expect a discrete spectrum. To find the spectrum of (4.13)–(4.14) we use a spectral approach, representing φ by a truncated sum of Chebyshev polynomials: φ = N∑ n=0 anTn (y) for y ∈ [−1, 1] (4.16) where N is the order of the truncated polynomial, an is the coefficient of the n-th Chebyshev polynomial, Tn (y). This method is described for example in Schmid & Henningson (2001) and is widely used. The discretised problem is coded and solved in Matlab. The accuracy of the code has been checked against the results of Mack (1976) for the Blasius boundary layer, with various results for PP flow in Schmid & Henningson (2001), with the PCP flow results of Potter (1966), and finally against results for PP flow with cross-flow; see 79 4.2. Stability of plane Couette-Poiseuille flow with cross-flow Sheppard (1972), Fransson & Alfredsson (2003). The results are accurate up to three, four and five significant places when validated against Potter (1966), Mack (1976) and Fransson & Alfredsson (2003), respectively. All the numerical results given below have been computed with N = 120. On using 200 collocation points, the growth rates changed only in the fourth significant place in the worst case. 4.2.3 Characteristic effects of varying k and Reinj Before starting a systematic analysis of (4.13)–(4.14), we briefly show some example results that illustrate the characteristic effects of varying the Couette component, k, and the cross-flow component, Reinj. These examples also serve to establish the framework of analysis used later in this chapter. With reference to PP flow, Potter (1966) first observed that the stability is increased by adding a Couette component while Fransson & Alfredsson (2003) showed that cross-flow can stabilize or destabilize PP flow. Eigenspectra Setting (α,Re) = (1, 6000), we investigate variations in the eigenspectrum of (4.13)–(4.14). According to a classification proposed by Mack (1976), the spectrum of PP flow spectra may be divided into 3 distinct families: A, P and S. Family A exhibits low phase velocity and corresponds to the modes concentrated near the fixed walls. Family P represents phase velocities, cr, close to the maximum velocity in the channel. Family S corresponds to the 80 4.2. Stability of plane Couette-Poiseuille flow with cross-flow mean modes and has phase velocity cr close to the mean velocity. In Figs 4.3(a) & (b), we track the eigenmodes as k and Reinj, respectively, are varied from zero. The initial condition (denoted by ) represents the PP flow. Referring to Fig. 4.3(a), (where Reinj = 0), addition of the Couette com- ponent increases the mean velocity: the S modes shift from cr = 0.6667 at k = 0 (PP flow) to cr = 0.7513 at k = 1. The family of P modes is also shifted to the right. The A modes are associated with both walls and as k increases we see a splitting of the family, with the upper wall modes moving to the right as k is increased. The least stable mode is a wall mode associated with the lower wall, which we observe stabilizes monotonically as k is increased. Figure 4.3(b) shows the effects of increasing Reinj, (holding k = 0). The least stable A mode of PP flow initially stabilizes and then destabilizes with increasing Reinj. This behavior has also been observed by Fransson & Alfredsson (2003). For large Reinj the A, P, and S families have disappeared, instead leaving two distinct families of modes. It appears that each of the A, P, and S families splits, with some modes entering each of the two families (this alternate splitting is most evident for the S modes). As observed by Nicoud & Angilella (1997), the phase speed no longer lies in the range of the axial velocity. This does not violate the conditions on cr, given by Joseph (1968) and Joseph (1969), since these conditions are derived for parallel flows only. 81 4.2. Stability of plane Couette-Poiseuille flow with cross-flow (a) (b) Figure 4.3: Eigenspectrum of (α,Re) = (1, 6000) by varying k and Reinj. 40 least stable modes are considered. (a) Effect of increasing k from 0 to 1 in steps of 0.01, keeping Reinj = 0. (b) Effect of increasing Reinj from 0 to 100 in steps of 0.05, keeping k = 0 (PP flow). The symbols in (a) and (b) are similar and are denoted as follows: k = 0 or Reinj = 0 by (), k = 1 or Reinj = 100 by (◦) and intermediate k or Reinj by (.). Note that the PP flow spectrum is represented by the  in both figures, and shows the vertical family of S-modes, the branch of A-modes (diagonally upwards from centre to left) and branch of P-modes (diagonally upwards from centre to right). 82 4.2. Stability of plane Couette-Poiseuille flow with cross-flow (a) (b) Figure 4.4: (a) Effect of increasing Reinj on the stability of PCP flow, for (α,Re) = (1, 6000) and different values of k = 0 (), 0.5 (◦), 1 (×). (b) Maxi- mal growth rate for increasing Reinj at different Re, (k = 0.5 and the step in values of Re between curves is 104). Increasing Reinj Next, we illustrate the qualitative effects of increasing Reinj at fixed (α,Re, k), in Fig. 4.4(a). We again fix α = 1 and Re = 6000, and show the variation of the least stable eigenvalue, for k = 0, 0.5, 1. Our results for k = 0 (PP flow) may be compared directly with those of Fransson & Alfredsson (2003). We observe that as Reinj increases we have an initial range of stabilization (ci,crit decreasing), followed by a range of destabilization (ci,crit increasing), and finally again stabilization at large Reinj, (ci,crit decreasing). Qualitatively, we have observed these same three ranges of decreasing/increasing ci,crit, as Reinj increases, for all numerical results that we have computed, and this provides a convenient framework within which to describe our results. For fixed (α,Re, k), the case Reinj = 0 may either be stable or unsta- ble, in which cases there are respectively two or three marginal stability val- 83 4.2. Stability of plane Couette-Poiseuille flow with cross-flow ues of Reinj. We denote these marginal values of Reinj by: Reinj,1, Reinj,2, Reinj,3, noting that in the case that Reinj = 0 is stable Reinj,1 is absent. More clearly, Reinj,2 will always represent a transition from stable to unstable, while Reinj,1 & Reinj,3 denote transitions from unstable to stable. The PCP flows for k = 0.5 and 1 are stable for (α,Re) = (1, 6000) in the absence of cross-flow, Reinj = 0. For a larger Re, k = 0.5 is unstable at Reinj = 0, but k = 1 remains stable for all (α,Re). Figure 4.4(b) shows the maximal growth rate γ, for increasing Reinj at dif- ferent Re, with k = 0.5. The maximal growth rate is computed over wavenum- bers α ∈ [0, 1]: γ = max α∈[0,1] {αci}, (4.17) which often captures the largest growth rates over all α. We observe that the first marginal value Reinj,1 increases with Re, but appears to converge towards a finite value as Re → ∞. The second marginal value of Reinj,2 appears independent of Re, (at least numerically). For k = 0.5 we have Reinj,2 ≈ 24.7. Nicoud & Angilella (1997) have observed a similar behaviour in studying the generalised Couette flow, (for which the constraint, kReinj = 4, is always satisfied). They have found Reinj,2 ≈ 24 (note that Nicoud and Angilella use the full channel width as their length-scale, and therefore report Reinj,2 ≈ 48, in their variables). In contrast, the third marginal value, Reinj,3, is strongly dependent on Re. For example, for k = 0.5, the values corresponding to Re = 10000 and 100000 are Reinj,3 ≈ 83 and Reinj,3 ≈ 287 respectively. 84 4.2. Stability of plane Couette-Poiseuille flow with cross-flow (a) (b) Figure 4.5: Maximal growth rate versus Reinj at Re = 40000: (a) Reinj,2 & Reinj,3 for k = 0 (◦) to 1 (); (b) Reinj,1 for k = 0 (◦) to 0.6 (). Step size is 0.2 in both figures. Increasing k Figure 4.5 explores the effects of increasing the Couette component k, on γ and on the marginal values of Reinj. Figure 4.5(a) indicates that the sensitivity of Reinj,2 to k is also not extreme: we have found that this transition occurs within the range ∼ 22 − 25 for k ∈ [0, 1]. For each value of k examined, we also observe numerically a similar independence of Reinj,2 to Re as seen earlier in Fig. 4.4(b) for k = 0.5. The 3rd marginal value, Reinj,3, is strongly dependent on k. For example, at Re = 40000, Reinj,3(k = 1) ≈ 120 and Reinj,3(k = 0.8) ≈ 135. In general, increasing k shifts Reinj,3 to the left, thereby decreasing the span of the unstable region. Increasing k also decreases the maximum value of γ. Figure 4.5(b) looks at the first transition, Reinj,1 at Re = 40000. Potter 85 4.3. PCP flows and the effects of small Reinj (1966) was the first to observe that for PCP flows (i.e. Reinj = 0), a gradual increase in the wall velocity results in crossing a “cut-off” value of k, say k1, such that for k > k1 the flow is unconditionally linearly stable. It has already been pointed out from the results of Fig. 4.4(b) that Reinj,1 is finite as Re→∞. In addition, the results in Fig. 4.5(b) indicate that Reinj,1 decreases with k at a finite Re. Hence, it can be inferred that as Re → ∞, the cut-off wall velocity, k1 = k1(Reinj) must decrease with Reinj. 4.3 PCP flows and the effects of small Reinj Having developed a broad picture of the different transitions occurring in the flow, we now focus in depth at each range of Reinj, to understand the stability mechanisms in play. We start with the range of small Reinj. PCP flows without cross-flow are stable to inviscid modes, but viscosity admits additional modes, i.e. the Tollmien-Schlichting (TS) waves, which may destabilize, according to the value of k. When αR 1 with c ∼ O(1), viscous effects occur in thin oscillatory layers: (i) adjacent to the walls, (of thickness ∼ (αRe)−1/2), and (ii) close to the critical point(s), yc, where u(yc) = cr,crit are found, (of thickness ∼ (αRe)−1/3). It is in the critical layers that we see peaks in the distribution of energy production, implying transfer from the base flow. Potter (1966) put forward the argument that for a dimensionless wall velocity that exceeds cr,crit, the critical layer near the moving wall will vanish and there remains only one critical layer, near the fixed lower wall. The thickness of this 86 4.3. PCP flows and the effects of small Reinj second layer increases with wall velocity, thereby favouring stabilization. This mechanism appears to correctly describe the long wavelength pertur- bations, (at Reinj = 0), which are found to be the least stable for k ∼ O(1). Indeed Cowley & Smith (1985) developed a long wavelength analysis (α ∼ Re), in order to predict the cut-off value k1(Reinj = 0) ≈ 0.7. For values k ∼ O(1), PCP flows have only a single neutral stability curve (NSC). However, Cow- ley & Smith (1985) noted that for smaller k, multiple neutral stability curves could exist, and at shorter wavelengths. For example, when 0 ≤ k ≤ Re−2/7 there is one NSC, when Re−2/7 ≤ k ≤ Re−2/13 there are three NSC’s, and when Re−2/13 ≤ k  1 there are two NSC’s; see Cowley & Smith (1985). Thus to understand the effect of cross-flow in PCP flows, the different regimes of k need to be considered separately. For Reinj ≈ 0, we expect the stability behaviour to be close to that of the PCP flow without cross-flow. Intuitively we expect cross-flow to stabilize, and so study the range cr,crit < k ≤ k1(Reinj = 0). We examine the NSC’s obtained from the O-S equation corresponding to k = 0.5, under different values of Reinj; see Fig. 4.6(a). As expected, increasing Reinj results in a progressively larger critical Re = Recrit. We also observe that both the upper and the lower branches are oriented at an angle of 45 degrees, (i.e. α ∼ Re−1), at high values of Re. On fixing Reinj and increasing k we have found that for successively large k the upper and lower branches move together as Recrit increases, eventually coalescing at k = k1(Reinj). This mechanism is identical with that observed by Cowley & Smith (1985), suggesting the applicability of 87 4.3. PCP flows and the effects of small Reinj (a) (b) Figure 4.6: Critical values for k = 0.5: (a) neutral stability curves for Reinj = 0 (×), 0.3 (◦) and 0.53 (); (b) variation in cr,crit with Reinj. Reinj αcrit Recrit cr,crit 0 0.3851 22600 0.2344 0.1 0.3576 22538 0.2370 0.2 0.3275 22924 0.2394 0.3 0.2950 23986 0.2415 0.4 0.2550 26321 0.2433 0.5 0.2000 31656 0.2452 0.6 0.1200 51115 0.2461 Table 4.1: Critical values for k = 0.5 and increasing Reinj. a long wavelength approximation in order to predict k1(Reinj). Figure 4.6(b) plots the values of cr at criticality, as Reinj is varied, also for k = 0.5. The critical values are tabulated in Table 4.1. The dependence is initially linear. We observe that k > cr,crit over the computed range. 88 4.3. PCP flows and the effects of small Reinj (a) (b) Figure 4.7: (a) Long wave NSC’s showing the dependence of λ on k for Reinj = 0 (dashed line), 0.3 (dash-dot), 0.5 (solid), 0.7 (dash-dot-dot) and 1 (long dash); (b) ci,crit versus k for λ = 2.5 × 10−5, and Reinj = 0 (dashed line), 0.3 (dash-dot), 1 (solid), 1.2 (dash-dot-dot) and 1.3 (long dash). 4.3.1 Long wavelength approximation We follow the long wavelength distinguished limit approach of Cowley & Smith (1985), taking α → 0 and Re → ∞, with λ = (αRe)−1 fixed. The product αRe is fixed along the upper and lower branches of the NSC. Thus, as the two branches of the NSC coalesce, in the (k, λ)-plane we observe k → k1(Reinj). In the long wavelength limit, equation 4.13 becomes: iλ [ D4 −ReinjD3 ] φ+ (u− c)D2φ− (D2u)φ = 0, (4.18) with boundary conditions (4.14). Figure 4.7(a) shows the NSC obtained from (4.18), plotted in the (k, λ)- plane for various Reinj. The cut-off value k1(Reinj) is the maximal value of k 89 4.3. PCP flows and the effects of small Reinj Reinj,1 k1 cr,crit k̃(k1, Reinj,1) 0 0.70 0.2331 0.5070 0.3 0.60 0.2431 0.4657 0.5 0.54 0.2455 0.4386 0.7 0.48 0.2472 0.4085 1.0 0.38 0.2358 0.3489 1.29 0.19 0.1556 0.1939 Table 4.2: Cut-off values, k1, and wavespeed cr,crit, for increasing Reinj. on each of these curves. These values are listed in Table 4.2. We also list the dimensionless wall speeds at cut-off, i.e. k̃(k1, Reinj). We observe that the cut- off wall speed decreases with Reinj. This is in agreement with the concluding remarks of §4.2.3. Figure 4.7(b) shows ci for the least stable eigenvalue of the long wavelength problem, for fixed λ = 2.5× 10−5 and different values of Reinj, as k is varied. When Reinj ≥ 1.3, we find that ci,crit ≤ 0, ∀ k ∈ (cr,crit, k1(0)], implying that there are no neutral or unstable long wavelength perturbations in this range of k, (i.e. at least until we approach the second transition at Reinj,2). Thus, in this initial range of say Reinj . 1.3, provided that k > cr,crit, we can talk equally of a cut-off value for k or for Reinj. 4.3.2 Effects of asymmetry of the velocity profile We observe that Reinj enters the stability problem in two distinct ways. The first one represents the direct contribution of the additional third order inertial term, ReinjD(α 2 − D2)φ, in the O-S equation (4.13). For the second one, 90 4.3. PCP flows and the effects of small Reinj (a) (b) Figure 4.8: Eigenspectrum for (k, α, Re) = (0.5, 0.2, 31656) (a) Reinj = 0.5 (Critical Conditions) and (b) Reinj = 23.5. Symbol ◦ indicates the eigen- spectrum from the O-S equation while  indicates the spectrum obtained by neglecting the additional cross-flow inertial term. Reinj influences the base velocity profile. To explore which of these effects is dominant, we show in Fig. 4.8 the spectra of (4.13)–(4.14) obtained with and without the term, ReinjD(α 2 −D2)φ, included in the computation. The critical parameters corresponding to Reinj = 0.5, in Table 4.1, are chosen and fixed for this comparison. Figure 4.8(a) shows the two spectra at Reinj = 0.5, which are near identical, completely overlapping on the figure. This suggests that at smaller Reinj, the effects of cross-flow manifest completely via the base flow velocity profile. Figure 4.8(b) shows a similar comparative study at a larger value of Reinj, closer to Reinj,2. In this figure we see a distinct difference between the spectra. The additional third order term is apparently responsible for the splitting of the A, P, and S families, illustrated earlier in Fig. 4.3(b). 91 4.3. PCP flows and the effects of small Reinj In Fig. 4.9, we plot k1 against Reinj(= Reinj,1). A linear dependence is evi- dent. The slope of the line is approximately −1/3. The flow is unconditionally linearly stable above the line and conditionally unstable otherwise. For small values of Reinj, we have seen in Fig. 4.1 that the principal effect is to skew the velocity profile towards the upper wall. A similar asymmetric skewing of the velocity profile is also induced in an annular Couette-Poiseuille (ACP) flow, through geometric means by varying the radius ratio, η, (defined as the radius of the outer stationary cylinder to the radius of inner moving cylinder). ACP flow has been studied extensively by Sadeghi & Higgins (1991), and we superimpose their results on ours, in Fig. 4.9. The comparison is striking. We believe there are 2 features of Fig. 4.9 that are unusual and worthy of note. Unsurprising is of course the identical limits Reinj = 0 = (η − 1). Note that Reinj → 0 is the PCP flow, and η → 1 represents the narrow gap limit of ACP, which is also the PCP flow. The first feature is the very similar linear decay in critical k = k1(Reinj), from the PCP values. It can be argued along the lines of Mott & Joseph (1968), that for a fixed Couette component (k), increasing the cross-flow for the PCP flow, or the radius ratio in the ACP flow of Sadeghi & Higgins (1991), skews the velocity profile more towards the moving boundary, thus increasing asymmetry and thereby stability. Since it has been already observed in Fig. 4.8(a) that for small Reinj the influence of injection on the eigenspectrum is through the velocity profile only, we do expect stabilization. However, when (η − 1) and Reinj are of O(1), we can see no obvious quantitative relation between these 92 4.3. PCP flows and the effects of small Reinj Figure 4.9: k1 as a function of Reinj,1 (shown by ) and the radius ratio, η (shown by ) in ACP flow (Sadeghi & Higgins (1991)) flows and even the stability operators are quite different. The second noteworthy feature of Fig. 4.9 is that there is a minimum value of k1 (k1,min) below which it is not possible to produce unconditional stability by applying (modest) cross-flow. This minimum value is found when k1 → cr,crit. We have found approximately that k1,min = 0.19 and the corresponding Reinj,1 = 1.29. This is very similar to Sadeghi & Higgins (1991), who found that the critical layer near the moving wall of ACP flows remained up to cr,crit ≈ 0.18. Linear energy budget considerations The strong analogy with the ACP results of Sadeghi & Higgins (1991) suggests that a similar mechanism may be responsible for the stabilization and cut- off behaviour. To investigate this we examine the linear energy equation, 93 4.3. PCP flows and the effects of small Reinj derived in modal form from the Reynolds-Orr energy equation. This yields the following two identities: ci = 〈(φrDφi − φiDφr)Du〉 − 1 αRe [I22 + 2α 2I21 + α 4I20 ] I21 + α 2I20 , (4.19) cr = 〈(α2|φ|2 + |Dφ|2)u〉+ Reinj αRe 〈α2(φrDφi − φiDφr) + (DφrD2φi −DφiD2φr)〉 I21 + α 2I20 (4.20) where Ik = Ik(φ) is the semi-norm defined by: Ik = [∫ 1 −1 |Dkφ|2 dy ]1/2 , k = 0, 1, 2, and where 〈f〉 = ∫ 1 −1 f(y) dy. Before proceeding further, we note that Reinj only appears indirectly in (4.19), reinforcing the assertion that for order unity Reinj, the principle contribution to stability of injection is via the mean flow. Indeed, in the long wavelength limits of cut-off k that we have studied, we have found values λ = (αRe)−1 . 10−4 for instability. Thus, in (4.20) the term directly involving Reinj has minimal effect on cr, explaining the observations in Fig. 4.8(a). The identity (4.19) can also be interpreted as an energy equation, in form: d dt 〈T1〉 = 〈T2〉 − 1 Re 〈T3〉 (4.21) 94 4.3. PCP flows and the effects of small Reinj where T1 = 0.5 (|Dφ|2 + α2 |φ|2) , d dt T1 = αciT1, (4.22) T2 = 0.5ατDu, τ = φrDφi − φiDφr, (4.23) T3 = 0.5( ∣∣D2φ∣∣2 + 2α2 |Dφ|2 + α4 |φ|2). (4.24) The left-hand side of (4.21) represents the temporal variation of the spatially averaged (one wavelength) kinetic energy. The first term on the right-hand side of (4.21) is the exchange of energy between the base flow and the distur- bance. The last term, ( 1 Re 〈T3〉 ) , represents the rate of viscous dissipation. At criticality, the two terms on the right-hand side balances each other, but the spatial distributions of T2 and T3/Re indicate where the energy is generated and dissipated in the channel. Sadeghi & Higgins (1991) extensively utilised this linear energy approach in studying the effect of k on stability of ACP flow. They found that increase in the value of k − cr,crit decreases the Reynolds stress (τ) near the moving wall until it becomes negative, hence stabilizing. The critical layer near the moving wall vanishes for k > cr,crit and as k increases the Reynolds stress becomes progressively negative within the critical layer at the fixed wall, but this behavior is destabilizing since the velocity gradient is negative there for ACP flow. Figures 4.10(a)-(d) examine the distribution of T2 and T3/Re for the least stable eigenmode for the parameters listed in Table 4.1, i.e. we fix k = 0.5 and 95 4.3. PCP flows and the effects of small Reinj (a) (b) (c) (d) Figure 4.10: Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to criticality at Reinj= (a) 0, (b) 0.2, (c) 0.4 and (d) 0.6. In all the cases, k = 0.5. Dash-dot-dot line with symbol  represents T2, dashed line with filled 4 represents 1ReT3 and solid vertical line represents the location of the critical layer. increase Reinj up to Reinj = Reinj,1 ≈ 0.6. The critical layer is marked with a vertical line. We observe that both the rate of energy transfer and the rate of viscous dissipation decrease with the cross-flow. Without cross-flow, T2 is positive and negative respectively in the lower (injection) and upper (suction) halves of the domain. Increasing the cross-flow decreases both the positive (near injection wall) and negative (near suction wall) peaks. The location of the critical layer also moves away from the injection wall due to the skewing of the velocity profile. When Reinj ≈ Reinj,1, 〈T2〉 and 1Re 〈T3〉 not only equalize 96 4.4. Intermediate Reinj and short wavelength instabilities but (since φ has been normalised), will have magnitudes O(α−1) since αRe = constant at cut-off; (see also Sadeghi & Higgins (1991)). This reduced energy budget as Reinj ≈ Reinj,1. This is the primary reason for the cut-off. 4.3.3 Summary For the range of small to order unity Reinj with k ≥ cr,crit, the flow instability is dominated by long wavelength perturbations. This instability mchanism exhibits a cut-off phenomenon characterised by a near linear boundary in the (Reinj, k)-plane. The initial cut-off mechanism is very similar to that for ACP, as studied by Sadeghi & Higgins (1991), combining skewing of the velocity profile, shifting of the critical layer and decay of the net perturbation energy. 4.4 Intermediate Reinj and short wavelength instabilities We now consider the range 0 ≤ k ≤ cr,crit, in which the critical layer at the upper wall is still present. We investigate its stability characteristics by adding cross-flow of intermediate strength (0 ≤ Reinj . 21), avoiding for the moment the second transition. It is intuitive that the presence of the critical layer will affect the stability behavior. To verify this we have studied the two extremities of the range of k considered, i.e. k = 0 (PP flow) and k = 0.18. The respective NSCs are shown in Fig. 4.11. It is evident that the presence of the critical 97 4.4. Intermediate Reinj and short wavelength instabilities (a) (b) Figure 4.11: Neutral Stability Curves (NSCs) for (a) k = 0 and (b) 0.18 at different Reinj. The symbols indicate different values of Reinj and are as follows: × → Reinj = 0, ◦ → Reinj = 6 in (a) and 4 in (b), → Reinj = 12 in (a) and 8 in (b) layers render shorter wavelength modes unstable. Yet, it is also observed that with Reinj in this intermediate range, the stability increases dramatically. We have been unable to make any advance analytically in this range of Reinj, and therefore have proceeded numerically. First we note that when we have considered k & 0.19 for the range of 1.3 < Reinj < 21, we have found that the least stable modes are long wavelength modes and that these are linearly stable. Thus, k & 0.19 appears to represent an absolute cut-off in this range of Reinj. For smaller k we have seen that the NSC’s occur with wavenumbers that are O(1) and apparently increasing with Reinj. Unlike the long wavelength problem, the asymptotic behaviour along the branches of the NSC’s is not easily treated. At fixed large Re, we are able to compute numerically a cut-off value of k for increasing Reinj, i.e. k = k1(Reinj, R). These cut-off curves do 98 4.4. Intermediate Reinj and short wavelength instabilities k1 Reinj,1 αcrit 0 20.8 3.5227 0.0225 20.6 3.7458 0.0450 20.0 3.9381 0.0675 18.6 3.9831 0.0900 15.6 3.4146 0.1125 14.4 3.2112 0.1350 13.4 3.2112 0.1575 12.6 3.2112 0.1800 11.8 3.2112 Table 4.3: Cut-off values evaluated for shorter wavelength instabilities for Re = 106. lie below k ∼ 0.19, but are not wholly independent of Re, at least within the range of Re up to which our numerical code is reliable, i.e. it is quite possible that these asymptote to a cut-off curve as Re → ∞, but we cannot reliably evaluate this limit numerically. As an example of this numerical cut-off, (at Re = 106), we have computed the cut-off values Reinj,1, as listed in Table 4.3 and shown in Fig. 4.12(a). For the range 1.4 < Reinj < 11.8, the cut-off is close to k ∼ 0.19. Although we see that the unstable wavenumbers increase with Reinj in Fig. 4.11, note that asymptotically as α→∞ the short wavelengths are stable. To see this, from (4.19) we bound 〈(φrDφi − φiDφr)Du〉 ≤ |Du|maxI0I1 ≤ 0.5|Du|max[αI20 + I21/α], 99 4.4. Intermediate Reinj and short wavelength instabilities Figure 4.12: Shorter wavelength cut-off showing k1 as a function ofReinj,1. The flow is linearly stable for Re ≤ 106 above the curve. The values in Table 4.3 are marked by . so that ci < 0 provided that: Re < |Du|max 2α2 , (4.25) (and better bounds are certainly possible). In Table 4.3 we note that the maximal critical wavenumber is in fact attained at an intermediate Reinj. 100 4.4. Intermediate Reinj and short wavelength instabilities 4.4.1 Behaviour of preferred modes for intermediate Reinj. In our preliminary results, (§4.2.3), we saw that at fixed values of (Re, k, α), increasing the Reinj led to regimes of stabilization, then destabilization, and then finally stabilization. For k ≥ cr,crit, only long wavelengths appear unsta- ble and how the cut-off values of k and Reinj vary in this regime are illustrated in Fig. 4.9. For the lower range of k, our results are primarily numerical, indi- cating a cut-off value k ≈ 0.19 for 1.3 . Reinj . 11.8 and then with decaying cut-off k for 11.8 . Reinj . 20.8, as illustrated in Fig. 4.12. Therefore, we have linear stability as we cross some cut-off frontier, k > k1(Reinj) in the (Reinj, k)-plane, (alternatively for Reinj > Reinj,1). We now consider what happens to the certain eigenmodes (preferred modes) as we extend the injection cross-flow up until the second critical Reinj. Our analysis up to now suggests that the behaviour may be different depending on whether we consider small or moderate k. In Fig. 4.13, we have plotted the locations of certain eigenmodes as Reinj is increased, by keeping the Reynolds number Re constant at 106. This gives us some idea of how cut-off behaviour changes with Reinj. Although the “preferred modes” are simply those we have selected, we implicitly mean modes that are involved in the transition from stable to unstable as one of our dimensionless parameters is varied (here Reinj), i.e. at some point a preferred mode becomes the least stable mode and then unstable. 101 4.4. Intermediate Reinj and short wavelength instabilities Figure 4.13(a) shows two eigenmodes corresponding to k = 0, (PP flow). A least stable long wavelength mode is tracked for α = 0.001, denoted by ‘A’. This mode is stable at Reinj = 0 and its stability increases further as Reinj increases up to around 1.7. However, further increases in Reinj destabilize this mode progressively until it becomes unstable at Reinj = 25. In the inset of Fig. 4.13(a) we have also plotted the least stable short wavelength mode at α = 3.5227. Such modes become unstable only under the influence of cross- flow of intermediate strength. This particular mode, (denoted by ‘B’), starts becoming unstable approximately when Reinj > 15, but recovers stability later for Reinj ≥ 20.8. This behavior is a direct consequence of the trajectory of the NSCs observed in Fig. 4.11(a). The preferred mode ‘B’ is the critical mode at cut-off, (see Table 4.3). Thus PP flow with cross- flow is unconditionally linearly stable in the range 20.8 ≤ Reinj . 25. For larger k, the stability behavior is primarily governed by the long wave- length modes, as shown in Fig. 4.13(b) for k = 0.5. The least stable mode corresponding to α = 0.01 is unstable for Reinj = 0, denoted mode ‘C’. This viscous mode becomes stable when Reinj increases to 0.6, which is indeed the cut-off value, i.e. Reinj,1. This is expected, according to Table 4.2. Mode ‘C’ is weakly damped and its stability increases for Reinj . 3, after which it starts destabilizing. The mechanism of this destabilization can probably be analysed along the lines of resonant interactions of the Tollmien-Schlichting (T-S) waves; see Baines et al. (1996). To show this interaction, we have traced the locus, (for Reinj = [7, 30]), of the least stable inviscid short wavelength mode ‘D’, at 102 4.4. Intermediate Reinj and short wavelength instabilities (a) (b) Figure 4.13: Behavior of preferred modes (belonging to different wavelengths and denoted by alphabets ‘A’-‘D’) under the influence of cross-flow with Re = 106. Symbols  and ◦ respectively imply the starting and the ending position of the preferred mode in the ci, cr plane, whereas the dots (‘.’) trace the locus. The difference in Reinj between consecutive dots is 0.1. (a) k = 0. Mode ‘A’ has α = 0.001 and is traced for Reinj=[0,25]. Mode ‘B’ has α = 3.5227 and is traced for Reinj=[15,21] (shown in the inset), the position at Reinj = 15 is marked by ‘∗’. (b) k = 0.5. Mode ‘C’ has α = 0.01 and is traced for Reinj=[0,30]. Mode ‘D’ has α = 2.5 and is traced for Reinj=[7,30]. α = 2.5. This mode, being inviscid, remains stable but has ci very close to zero as Reinj increases. The wave speed cr decreases continuously with Reinj for mode ‘D’. The resonant interaction takes place when its wave speed matches with that of mode ‘C’, which signals the destabilization of mode ‘C’. This destabilization continues until mode ‘C’ becomes unstable when Reinj & 30. In Fig. 4.14 we show examples of the streamfunction for the preferred modes, corresponding to various k and Reinj in the transitions of Fig. 4.13. For the long wavelength mode ‘C’, Figs. 4.14(a)-(c) show that strong Reinj appears to skew the streamlines towards the lower wall. The same is true for the long wavelength mode ‘A’ under strong injection; see Fig. 4.14(f). On the 103 4.4. Intermediate Reinj and short wavelength instabilities Figure 4.14: Isovalues of the normalised perturbation stream functions (ψ′) for the preferred modes at Re = 106 under different Reinj. The streamwise extent of the domain is one wavelength. Corresponding to k = 0.5, the long wavelength mode ‘C’ is shown for (a) Reinj = 0.1 (unstable), (b) Reinj = 1 (stable) and (c) Reinj = 30 (unstable). Corresponding to k = 0, ψ ′ for two different preferred modes, viz. ‘A’ and ‘B’ are shown. The shorter wavelength mode ‘B’ (α = 3.5227) is shown for (d) Reinj = 15 (unstable) and (e) Reinj = 21 (stable). The longer wavelength mode ‘A’ (α = 0.001) is shown for (f) Reinj = 25 (unstable). 104 4.5. Stability and instability at large Reinj. contrary, Figs. 4.14(d)-(e) show that for large Reinj, the streamlines of the shorter wavelength mode ‘B’ are skewed and localised towards the upper wall. 4.5 Stability and instability at large Reinj. We turn now to the transition to instability at Reinj,2 and then later to sta- bilizing effects at very large Reinj. As observed in §4.2.3, the transition at Reinj,2 appears to be independent of streamwise Reynolds number Re (see Fig. 4.4(b)) and occurs for all k. Although there is sensitivity to k, it is not very significant. Values of Reinj,2 are found for all k ∈ [0, 1] and are in a fairly tight range of Reinj ∼ 22− 25. As suggested in the previous subsection, although instability at moderate Reinj may be either short wavelength or long wavelength, according to (k − cr,crit), as we approach Reinj,2 from below it is the long wavelengths that are unstable. Figure 4.15a shows the neutral stability curves corresponding to PP flow for Reinj just above Reinj,2. The NSC’s are nested with decreasing Reinj and as we approach Reinj,2 the upper and lower branches of the NSC are seen to coalesce. The slope of the two branches suggests that α ∼ R−1 in the limit of cut-off, and hence the previous long wavelength approximation, leading equation (4.18), should be effective for predicting the cut-off in the (k, λ)-plane; (recall λ = (αR)−1). Figure 4.15b shows the NSC’s obtained from long wave approximation. The cut-off velocity k2 is the maximum value of k encountered along the NSC 105 4.5. Stability and instability at large Reinj. (a) (b) Figure 4.15: (a) NSC of PP flow (k = 0) when Reinj → Re−inj,2. The different values of Reinj are 22.5 (dashed line with ×), 23 (dash-dot line with ) and 24 (dash-dot-dot line with ◦). Near cut-off, αRe is constant along the upper and lower branches. (b) Long wave NSCs showing the dependence of log10 λ on k. The different values of Reinj are 22.4 (), 23.5 (◦), 24 () and 25 (×). Cut-off is achieved over the en- tire range of k, i.e. [0, 1]. for a given Reinj = Reinj,2. Unlike Fig. 4.7, the entire range of k becomes unconditionally stable. For Reinj,1 < Reinj < 22.2, ci < 0 ∀ k ∈ [0, 1]. Another significant difference with Fig. 4.7 and the results of Cowley & Smith (1985) is that “bifurcation from infinity” is not observed as k → 0. This is possibly because the curves bifurcate from infinity for negative values of k, but we have not studied this range. Finally, we mention that for cross-flow rates slightly greater than Reinj,2, the Recrit is relatively low for the entire range of k. For example, Reinj,2 ≈ 23.8 for k = 0.5, (implying Recrit → ∞ as Reinj → Re−inj,2). Increasing Reinj to 25 decreases Recrit to around 6000. Thus, on crossing Reinj,2 we find a dramatic decrease in the flow stability. 106 4.5. Stability and instability at large Reinj. (a) (b) Figure 4.16: (a) Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to criticality at Reinj = 25. Dash- dot-dot line with symbol  represents T2, dashed line with filled 4 represents 1 Re T3 and solid vertical line represents the location of the critical layer. (b) Reynolds Stress τ distribution at criticality for Reinj = 0 (denoted by  sym- bol), Reinj = 0.6 (denoted by ×) and Reinj = 25 (denoted by ◦). The location of the critical layers are shown by solid lines with corrosponding symbols. 4.5.1 Linear energy balance at Reinj,2 An interesting feature of transition at Reinj,2 is the independence with re- spect to Re. With reference to the energy equation (4.21), this insensitiv- ity implies that in this range |T2| is much larger than the viscous dissipa- tion, 1 Re T3. In other words, at criticality ci = 0 is achieved by a balance of energy production and dissipation within T2, more so than via balance with the viscous dissipation. Figure 4.16(a) investigates the energy budget at criticality for k = 0.5 at Reinj = 25. The critical parameters are ob- served to be (αcrit, Recrit) = (0.31, 6000). This implies that crossing the cut-off Reinj,2, there is a transition from unconditional stability (Recrit →∞) to high 107 4.5. Stability and instability at large Reinj. instability(Recrit = 6000). Comparing with Fig. 4.10 (which shows energy distribution corresponding to criticality for k = 0.5 and Reinj ≤ Reinj,1) , it is obvious that T2 has a higher amplitude while the viscous dissipation 1 Re T3 is weaker. This behavior is due to the generation of larger Reynolds stresses τ , as Reinj increases, as illustrated in Fig. 4.16(b). The dominance of T2 over the viscous dissipation suggests that the critical layers have little to do with instability in this range. Note that τ is small in the critical layer, which has now moved towards the channel centre, and hence T2 is also small. Referring to Fig. 4.2(b), the vanishing vorticity gradient (D2u) found in the bulk of the flow domain at high values of Reinj removes/diminishes the singular effects associated with the critical layer. The growth of τ is probably not responsible for the spreading of the spec- trum along the real axis, that we have observed in Fig. 4.3(b). Equation (4.20) may be rewritten as: cr = 〈(α2|φ|2 + |Dφ|2)u〉+ Reinj αRe 〈α2τ − φrD3φi + φiD3φr〉 I21 + α 2I20 . (4.26) The first term leads simply to values of cr in the range of u. The second term does contain α2τ , i.e. longitudinal gradients of the Reynolds stresses. However, note that even for the shorter wavelengths we have α ∼ O(1), and if we consider long wavelengths, we have typically found instability only for 108 4.5. Stability and instability at large Reinj. (a) (b) (c) (d) (e) (f) Figure 4.17: Distribution of energy production (T2) and dissipation ( 1 Re T3) terms across the domain corresponding to mode ‘C’ at Reinj= (a) 0, (b) 0.6, (c) 1, (d) 3, (e) 10 and (f) 30. Dash-dot-dot line represents T2, solid line represents 1 Re T3. αRe  1. Thus, even for these larger Reinj, the term involving α2τ is likely to be insignificant. The extension of cr beyond the usual bounds of the base flow velocity is therefore due to the 3rd derivative terms in (4.26), which cannot be bounded by the denominator. Interestingly therefore, the larger values of cr, which indicate less regular eigenmodes, also lead to larger viscous dissipation, and hence more stable modes. This explains the shape of the spectrum in Fig. 4.3(b). In Fig. 4.13(b), we tracked the behaviour of mode ‘C’ as Reinj increased. 109 4.5. Stability and instability at large Reinj. This mode becomes unstable for Reinj ≥ Reinj,2, implying that it governs the transition behavior. In Fig. 4.17 we show the evolution of the energy balance terms for this mode as Reinj increases from zero. This mode is stable from Reinj = 0.6 to 30. The cut-off achieved at Reinj = 0.6 is primarily due to the increased viscous dissipation at both walls. This phenomenon continues until Reinj ≈ 3, see Fig. 4.17(d). At this point, the ‘viscous hump’ observed near the lower wall gets amplified. This mechanism is probably due to the resonant interaction between mode ‘C’ and an (approximately) neutrally stable inviscid mode ,for example mode ‘D’. Further increase in Reinj thins out the viscous layer at the suction wall faster than that at the injection wall. Suction negates both the exchange as well as the dissipation of energy, and the viscous hump is localised within the lower half of the channel, i.e. injection side. 1 Re T3 reduces faster than T2 and finally the mode becomes unstable when Reinj increases to 30. The condition at this point is 〈T2〉 > 1Re 〈T3〉; see Fig. 4.17(f). It is interesting to observe that the mode becomes unstable when the viscous hump reaches the centre of the channel. Further increase of Reinj results in a gradual reduction of T2 and the mode becomes stable again. The mean perturbation kinetic energy q(y) distribution provides further insight into the instability mechanism. It is defined modally to be: q = 1 4 (|Dφ|2 + α2 |φ|2) (4.27) Figure 4.18 shows the mean perturbation kinetic energy profiles for mode ‘C’ 110 4.5. Stability and instability at large Reinj. at different Reinj. Each distribution of q has been normalised by its maximum value. Without any cross-flow, the amount of energy in the two halves of the domain are comparable, the suction half having ∼ 43% of the energy, (note k = 0.5). Increasing cross-flow up to Reinj ≈ Reinj,1 = 0.6 increases the secondary peak until the cut-off is achieved. The energy in the suction half at this point is 46.7%. The primary peak moves toward the lower wall but cannot reach it because of the no-slip conditions. At Reinj > 1, the primary peak starts moving away from the suction wall. At Reinj = 3, the mode is at its maximum stability (see Fig. 4.13(b)). At this point, the perturbation energy is highly localised within the lower 1 8 th of the channel, along with a small secondary peak at the upper quarter. Further increase of Reinj to 10 causes the secondary peak to vanish; the energy content in the suction half being only ∼ 7.6%. The resonant interactions of T-S waves result in the development of a secondary peak from the primary peak itself. During this process, the secondary peak slowly separates from the primary peak and moves in the direction of the upper wall. For Reinj = 30, the perturbation reaches the channel centre and the mode becomes unstable. The amount of energy in the suction half increases to 18.1%. For even higher values of Reinj, for example 45, the upper half holds ∼ 32% of the energy. Thus it appears that the onset of the cut-off at Reinj,1 occurs when the secondary peak holds maximum energy. Increasing injection decays this peak until it reaches a minimum and then starts to grow out from the primary 111 4.5. Stability and instability at large Reinj. Figure 4.18: Non-dimensional mean perturbation kinetic energy profiles for mode ‘C’ at different Reinj. Solid lines with symbols denote the unstable modes. For each Reinj, q has been scaled by its maximum value. peak. The end of the cut-off regime, marked by Reinj > Reinj,2, occurs when the secondary peak reaches the channel centre and holds sufficient energy. 4.5.2 Eventual stabilization at Reinj,3 We have not studied in detail the eventual stabilization of the flow at very large Reinj (i.e. Reinj ∼ Reinj,3), but we believe the energetics of this stabi- lization are due to a decay in the energy production. This can be seen most clearly from the identity (4.19), which is in the same form as that for any parallel shear flow, i.e. cross-flow only influences (4.20) directly. Joseph has used this expression to derive general bounds that depend on |Du|max, and 112 4.5. Stability and instability at large Reinj. various functional inequalities; see Joseph (1968, 1969). For example, we have linear stability provided that: αRe|Du|max < max(ξ2pi + 23/2α3, ξ2pi + α2pi) (4.28) where ξ = 2.36502 is the least eigenvalue of a vibrating rod with clamped ends at y = ±1. The condition (4.28) evidently holds for the flows we consider, but is very conservative and especially so in the limit of large Reinj. This conservatism at large Reinj stems directly from the simplistic treatment of Du in bounding the energy production term: 〈(φrDφi − φiDφr)Du〉 < |Du|maxI0I1. With reference to Figs. 4.1 & 4.2 and to (4.11), we see that at large Reinj the base velocity profile consists of a thin layer near the upper suction wall, within which Du ∼ |Du|max ∼ Reinj, which has thickness of O(Re−1inj). Away from this thin boundary layer, the velocity gradients are of size Du ∼ 2(kReinj)−1 + O(Reinje −Reinj(1−y)). Note however, that within this suction layer, we have φ ∼ (1− y)2 due to the boundary conditions on the perturbation. Therefore, 113 4.6. Summary taking a nominal suction layer boundary at y = ys, we may estimate as follows: 〈(φrDφi − φiDφr)Du〉 = ∫ 1−ys −1 (φrDφi − φiDφr)Du dy + + ∫ 1 1−ys (φrDφi − φiDφr)Du dy ≤ 2 kReinj ∫ 1−ys −1 |φrDφi − φiDφr| dy + O(|Du|max(1− ys)4) ≤ 2 kReinj I0I1 +O(Re −3 inj) (4.29) Following Joseph (1969), this leads directly to the bound 2αRe kReinj . max(ξ2pi + 23/2α3, ξ2pi + α2pi), (4.30) sufficient for linear stability at large Reinj, (with asymptotically kReinj & 4 required). In other words, at large Reinj, the energy production T2 will decay like (kReinj) −1 at leading order, so that the viscous dissipation need only be of this order to stabilize the flow. 4.6 Summary To summarise, we have presented a detailed analysis of linear stability and instability in the (Reinj, k)-plane, for PCP flow with cross-flow. The most complete analysis concerns the important range of low Reinj and modest k. In this range we have demonstrated that the stabilization mechanism, due to 114 4.6. Summary Figure 4.19: Variation of k with Reinj. The filled  symbols show the long wavelength cut-off achieved for 0.7 ≥ k ≥ 0.19. The filled ◦ symbols show the shorter wavelength cut-off for 0.19 > k ≥ 0 evaluated numerically for Re = 106. The filled  symbols imply the second long wavelength cut-off. The shaded region depicts the entire zone of unconditional linear stability. either injection or wall motion, is essentially the same. Long wavelengths dom- inate. Skewing of the velocity profile shifts the critical layer and at the same time the energy production is diminished until viscous dissipation dominates at cut-off. In Fig. 4.9, we have also shown an interesting quantitative analogy with the cut-off behavior of ACP flows; see Sadeghi & Higgins (1991). This lower range of Reinj and modest k is probably that which is most important practically. Essentially, this range allows one to compensate cross- flow by wall-motion and vice-versa, achieving unconditional linear stability via either mechanism. With reference to Fig. 4.1, it is the range of Reinj in which 115 4.6. Summary the cross-flow and wall motion are modifications of a base Poiseuille flow. Due to the scaling, the peak velocity is always 1, but at larger Reinj with modest k the Poiseuille component is completely dominated by cross-flow and wall motion. Globally, the cut-off regimes in the (Reinj, k)-plane are as illustrated in Fig. 4.19. The shaded area shows the region of unconditional linear stability. In the intermediate range of approximately 1.3 ≤ Reinj ≤ 20.8 values of k & 0.19 are dominated by long wavelengths and are stable. Below this value, we are able to compute numerical cut-off curves for fixed Re. With the limits of our computations, we cannot determine if these cut-off curves asymptote to an unconditional cut-off curve as Re→∞. There appears to be a short band of unconditional linear stability for all computed values of k around approximately 20.8 ≤ Reinj ≤ 22, before the destabilization occurs at larger Reinj = Reinj,2. Since this band can make PP flow unconditionally stable, it could be effectively used in applications where wall motion is not feasible, e.g. cross-flow filtration, medical dialysis. From the practical perspective, it is worth noting that the transition across Reinj,2, is from unconditional stability to critical values of Re which are relatively modest (e.g. in the range 103 − 104) just a short distance beyond Reinj,2. Assuming that the PP flow is linearly unstable, this means stabilization can be achieved with cross-flow velocities of the order of 1% of the mean axial flow velocity. This destabilization at Reinj,2 is again a long wavelength mechanism, which we have analysed using the long wavelength approximation of Cowley & Smith 116 4.6. Summary (1985). A possible cause of this instability has been found to be resonant interactions of the T-S waves. Study of the linear energetics of the upper limit, Reinj,2, has shown that neither viscous dissipation, nor the involvement of a critical layer are significant. Rather, the balance of energy production and dissipation within T2 keeps the mode neutrally stable. Energy analysis of the preferred mode ‘C’ has revealed that the precursor of the transition to instability from unconditional stability is the amplification of disturbances near the injection wall. The mean perturbation kinetic energy has also been analyzed. It has shown that the lower limit occurs when the secondary peak holds maximum energy. Increasing injection decreases the secondary peak until it reaches a minimum and then it starts to grow from the primary peak. When the secondary peak reaches the channel centre and holds a sufficient amount of energy, the unconditional stability mechanism breaks down. The final stabilization occurring at large Reinj ≥ Reinj,3 has been analyzed using linear energy bounds. By careful treatment of the energy production term, we are able to show that the energy production terms decreases asymp- totically like Re−1inj as Reinj → ∞. We believe that this mechanism leads to the eventual domination of the viscous dissipation at large enough Reinj. In terms of the spatial structure of the perturbations, we note that the stabilization at small and moderate Reinj are both long wavelength phenomena for which the approximation of Cowley & Smith has been shown effective. Implicitly therefore, the critical wavenumbers scale like Re−1 in these limits. For the shorter wavelength instabilities we have not analysed the asymptotic 117 4.6. Summary behaviour of the wavenumber with Re. A more detailed look at the spatial structure of certain eigenmodes has been presented in Fig. 4.14. This shows a skewing of the streamline recirculatory regimes towards the lower wall for long wavelengths as Reinj is increased, and towards the upper wall at shorter wavelengths as Reinj is increased. 118 Chapter 5 Conclusions Education never ends Watson. It is a series of lessons with the greatest for the last. — Sherlock Holmes in The Red Circle. 5.1 Summary In this thesis we have investigated various aspects of two dimensional shear flow instabilities. In Chapter 2 we have investigated the physical reason behind the exponential amplification of an infinitesimal disturbance in an otherwise stable shear flow. By considering idealized shear layers, we have shown that hydrodynamic instability occurs because two (or more) linear interfacial waves, having arbitrary initial amplitudes and phases, interact with each other in such a way that they eventually attain a resonant configuration. The latter provides the condition for (idealized) shear instabilities; it furnishes with the range of unstable wavenumbers. This generalized wave interaction approach (referred to as wave interaction theory (WIT)) has been validated against three different types of shear instabilities, Kelvin-Helmholtz (KH), Holmboe, and Taylor. WIT also provides a non-modal description of idealized shear instabilities. 119 5.1. Summary Non-modal instability signifies non-orthogonal interaction between the two wave modes, and is the entire process occurring prior to resonance. Rapid transient growth is a key feature of non-modal instability process. WIT shows that optimal growth occurs when the two waves are in quadrature. In Chapter 3, the idea proposed in Chapter 2 has been extended to the non-linear regime for the case of KH. It has been shown that the non-linear interaction between two counter-propagating vorticity waves produce elliptical vortices. Contrary to the common notion that KH manifests itself in the form of spiral billow structures, our simulation has shown that elliptical vortices are another asymptotic form of KH. The dynamics of such elliptical vortices (rotation and nutation) as well as the attached thin braids have been investi- gated. Simple models are shown to provide a leading order description of the vortex and braid dynamics. It has been hypothesized that elliptical vortices in geophysical flows arise due to this phenomenon. In Chapter 4 we have followed the conventional normal-mode approach to investigate how channel flows can be unconditionally stabilized/destabilized. A channel flow can be unconditionally stabilized by making one of the walls to move (i.e. adding ”Couette” component). Previous studies have shown that adding small amount of cross-flow (suction from one wall and injection from the other) produces a stabilizing influence, but little is known about the influence of moderate or large cross-flows. In a parameter space of non-dimensional wall velocity and non-dimensional cross flow velocity, we have investigated how these two parameters influence the stability of a channel flow. It has been 120 5.2. Main contributions shown that the two effects compensate each other (in the parameter space considered here). Most importantly, it has been found that there exists a small band of cross-flows for which the channel flow is unconditionally linearly stable even in the absence of wall motion. 5.2 Main contributions The contributions of each Chapter are as follows: 5.2.1 Chapter 2 Here we have been able to provide a mechanistic understanding of shear insta- bilities in terms of interfacial wave interactions. Moreover we have showed that in the case of idealized shear instabilities, normal modes of spectral analysis, resonant condition of wave mechanics, and equilibrium points of dynamical systems can be brought under one umbrella. We have formulated a condition for idealized, homogeneous or stratified shear instabilities. Our formulation provides a non-modal description of instabilities, hence transient growth pro- cesses can be well understood in the light of this theory. 5.2.2 Chapter 3 The main contribution of this work is to bridge contour dynamics with wave interaction theory. We have shown that contour dynamics is the non-linear extension of WIT, and using it we have studied the non-linear evolution of a 121 5.2. Main contributions piecewise linear shear layer. It is found that this shear layer produces elliptical vortex. The dynamics of this vortex is well predicted by the theory proposed by Kida. By comparing our results with a number of geophysical vortices, e.g. meddies, Jupiter’s Great Red Spot, Neptune’s Great Dark Spot, we have conjectured that elliptical shaped geophysical vortices probably arise due to the above-mentioned mechanism. 5.2.3 Chapter 4 The main contribution of this chapter is finding the range of non-dimensional wall velocities and non-dimensional injection Reynolds numbers for which the two effects compensate each other. This finding can help developing a robust compensatory design of flow stabilization using either mechanism. The most important contribution of this work is finding a range of non-dimensional in- jection Reynolds numbers for which the flow is unconditionally stable even in the absence of wall motions. Practically speaking, wall motions are not a feasible option for flow stabilization. Suction-injection on the other hand is comparatively easier to implement. The study can therefore be valuable for various engineering applications where stabilizing the flow is crucial. 122 5.3. Future research 5.3 Future research 5.3.1 Chapter 2 Wave interaction theory proposed in Chapter 2 is only valid in the linear regime. In future, this theory needs to be extended into the non-linear regime to understand finite amplitude modal and non-modal instability mechanisms. Although we have extended the KH instability into non-linear domain in Chap- ter 3, Holmboe and Taylor instabilities need to be considered in future. Although we have limited our study to the interaction between two waves, the Taylor and Holmboe profiles actually involve multiple wave interactions, which we have neglected. Each density interface in these profiles supports two gravity waves, out of which only the counter-propagating wave has been considered. Noting that two wave interactions are sufficient to produce the normal mode characteristics of Holmboe and Taylor instabilities, it can be argued that the inclusion of co-propagating gravity wave would have been unnecessary. However this wave might have some effect during the initial interaction (non-modal) stages, which needs to be studied in future. 5.3.2 Chapter 3 To better understand the formation of geophysical vortices, more realistic ve- locity profiles need to be studied. In such flows, the dynamics of the elliptical vortices will be affected by strain-rate as well as background shear. Another interesting aspect worth studying is the nature of the ensuing KH when the 123 5.3. Future research velocity profile is gradually varied from piecewise linear to hyperbolic tangent. Moreover, our analysis is mainly restricted to the early post-saturation phase. We have briefly studied the late post-saturation phase when Kelvin waves ap- pear on the elliptical core. Future study may provide valuable insight into the sustainability of geophysical vortices. 5.3.3 Chapter 4 Future experimental and numerical research needs to be undertaken to shed light into the stabilization of the channel flow by the application of cross-flow. Linear stability analysis has provided us with a range of injection Reynolds number for which the flow is unconditionally stable. Since viscous normal mode stability theory does not correspond very well with experiments or com- putations, It is unlikely to get the exact predicted range in practice. However linear theory is known to provide an estimate, which can be favorably utilized while designing experiments or performing numerical analysis. 124 Bibliography Baines, P.G., Majumdar, S.J. & Mitsudera, H. 1996 The mechanics of the tollmien-schlichting wave. Journal of Fluid Mechanics 312 (107), 107– 124. Baines, P.G. & Mitsudera, H. 1994 On the mechanism of shear flow instabilities. J. Fluid Mech. 276, 327–342. Bretherton, F. P. 1966 Baroclinic instability and the short wavelength cut-off in terms of potential vorticity. Q. J. Roy. Meteor. Soc. 92 (393), 335–345. Cairns, R.A. 1979 The role of negative energy waves in some instabilities of parallel flows. J. Fluid Mech. 92, 1–14. Carpenter, Jeffrey R., Tedford, Edmund W., Heifetz, Eyal & Lawrence, Gregory A. 2013 Instability in stratified shear flow: Review of a physical interpretation based on interacting waves. Applied Mechanics Reviews 64 (6), 061001. Caulfield, C.P. 1994 Multiple linear instability of layered stratified shear flow. J. Fluid Mech. 258, 255–285. Caulfield, C.P., Peltier, W.R., Yoshida, S. & Ohtani, M. 1995 An experimental investigation of the instability of a shear flow with multilayered density stratification. Phys. Fluids 7, 3028–3041. Chérubin, L. M., Serra, N. & Ambar, I. 2003 Low-frequency variabil- ity of the mediterranean undercurrent downstream of portimão canyon. J. Geophys. Res. 108 (C3), 3058. Constantinou, Navid C. & Ioannou, Petros J. 2011 Optimal excitation of two dimensional holmboe instabilities. Phys. Fluids 23 (7), 074102. 125 Bibliography Corcos, G. M. & Sherman, F. S. 1976 Vorticity concentration and the dynamics of unstable free shear layers. J. Fluid Mech. 73 (02), 241–264. Cowley, SJ & Smith, FT 1985 On the stability of poiseuille-couette flow: a bifurcation from infinity. Journal of Fluid Mechanics 156, 83–100. Davidson, P. A. 2004 Turbulence: an introduction for scientists and engi- neers . Oxford University Press, USA. Davies, HC & Bishop, CH 1994 Eady edge waves and rapid development. J. Atmos. Sci. 51 (13), 1930–1946. Deem, Gary S. & Zabusky, Norman J. 1978 Vortex waves: Stationary ”v states,” interactions, recurrence, and breaking. Phys. Rev. Lett. 40, 859–862. Drazin, P.G. & Reid, W.H. 2004 Hydrodynamic Stability , 2nd edn. Cam- bridge University Press. Farrell, B. 1984 Modal and non-modal baroclinic waves. J. Atmos. Sci. 41 (4), 668–673. Farrell, B.F. & Ioannou, P.J. 1996 Generalized stability theory. part i: Autonomous operators. J. Atmos. Sci. 53 (14), 2025–2040. Fransson, J.H.M. & Alfredsson, P.H. 2003 On the hydrodynamic sta- bility of channel flow with cross flow. Physics of fluids 15, 436–441. Goldstein, S. 1931 On the stability of superposed streams of fluids of dif- ferent densities. Proc. R. Soc. Lond. A 132, 524–548. Green, MA 2006 Turbulent channel structures, princeton university art of science gallery. Hains, FD 1967 Stability of plane couette-poiseuille flow. Physics of Fluids 10, 2079–2080. Hains, FD 1971 Stability of plane couette-poiseuille flow with uniform cross- flow. Physics of Fluids 14 (8), 1620–1623. Hazel, P. 1972 Numerical studies of the stability of inviscid stratified shear flows. J. Fluid Mech. 51, 39–61. 126 Bibliography Heifetz, E., Bishop, CH & Alpert, P. 1999 Counter-propagating Rossby waves in the barotropic Rayleigh model of shear instability. Q. J. R. Mete- orol. Soc. 125 (560), 2835–2853. Heifetz, E., Bishop, CH, Hoskins, BJ & Methven, J. 2004 The counter-propagating rossby-wave perspective on baroclinic instability. i: Mathematical basis. Q. J. Roy. Meteor. Soc 130 (596), 211–231. Heifetz, Eyal & Methven, John 2005 Relating optimal growth to coun- terpropagating Rossby waves in shear instability. Phys. Fluids 17 (6), 064107. Hendel, Tal 2008 http://www.mathworks.com/matlabcentral/fileexchange/22423- ellipse-fit. Holmboe, J. 1962 On the behavior of symmetric waves in stratified shear layers. Geofys. Publ. 24, 67–112. Hoskins, B.J., McIntyre, M.E. & Robertson, A.W. 1985 On the use and significance of isentropic potential vorticity maps. Q. J. Roy. Meteor. Soc. 111 (470), 877–946. Joseph, D.D. 1968 Eigenvalue bounds for the orr-sommerfeld equation. J. Fluid Mech 33 (part 3), 617–621. Joseph, D.D. 1969 Eigenvalue bounds for the orr-sommerfeld equation. part 2. J. Fluid Mech 36 (part 4), 721–734. Kida, Shigeo 1981 Motion of an elliptic vortex in a uniform shear flow. J. Phys. Soc. Jpn. 50 (10), 3517–3520. Kundu, PK & Cohen, IM 2004 Fluid Mechanics . Elsevier, Boston. Liu, J. & Schneider, T. 2010 Mechanisms of jet formation on the giant planets. J. Atmos. Sci. 67, 3652–3672. Mack, L.M. 1976 A numerical study of the temporal eigenvalue spectrum of the blasius boundary layer. Journal of Fluid Mechanics 73 (3), 497–520. Mitchell, T. B. & Rossi, L. F. 2008 The evolution of kirchhoff elliptic vortices. Phys. Fluids 20 (5), 054103. 127 Bibliography Morales-Jubeŕıas, R., Sánchez-Lavega, A., Lecacheux, J. & Co- las, F. 2002 A comparative study of jovian anticyclone properties from a six-year (1994–2000) survey. Icarus 157 (1), 76–90. Mott, J. & Joseph, D.D. 1968 Stability of parallel flow between concentric cylinders. Physics of Fluids 11 (10), 2065–2073. Nicoud, F. & Angilella, JR 1997 Effects of uniform injection at the wall on the stability of couette-like flows. Physical Review E 56 (3), 3000–3009. Orszag, S.A. 1971 Accurate solution of the orr-sommerfeld stability equa- tion. J. Fluid Mech 50 (4), 689–703. Polivani, L. M., Wisdom, J., DeJong, E. & Ingersoll, A. P. 1990 Simple dynamical models of neptune’s great dark spot. Science 249, 1393– 1398. Potter, M.C. 1966 Stability of plane couette-poiseuille flow. J. Fluid Mech 24 (4), 609–619. Pozrikidis, C. 1997 Introduction to Theoretical and Computational Fluid Dynamics , first edition edn. Oxford University Press. Pozrikidis, C. & Higdon, J. J. L. 1985 Nonlinear kelvin–helmholtz insta- bility of a finite vortex layer. J. Fluid Mech. 157, 225–263. Pullin, D. I. 1992 Contour dynamics methods. Annu. Rev. Fluid Mech. 24 (1), 89–115. Rayleigh, J.W.S. 1880 On the stability, or instability, of certain fluid mo- tions. Proc. Lond. Math. Soc. 12, 57–70. Reynolds, WC & Potter, M.C. 1967 Finite-amplitude instability of par- allel shear flows. J. Fluid Mech 27 (3), 465–492. Romanov, VA 1973 Stability of plane-parallel couette flow. Functional Anal- ysis and its Applications 7 (2), 137–146. Sadeghi, V.M. & Higgins, B.G. 1991 Stability of sliding couette–poiseuille flow in an annulus subject to axisymmetric and asymmetric disturbances. Physics of Fluids A: Fluid Dynamics 3 (9), 2092–2104. 128 Bibliography Saffman, P. G. 1995 Vortex Dynamics , first edition edn. Cambridge Uni- versity Press. Schmid, P.J. & Henningson, D.S. 2001 Stability and transition in shear flows , , vol. 142. Springer Verlag. Sheppard, D.M. 1972 Hydrodynamic stability of the flow between parallel porous walls. Physics of Fluids 15 (2), 241–244. Smyth, W. D. & Moum, J. N. 2012 Ocean mixing by kelvin-helmholtz instability. Oceanography 25, 140–149. Squire, HB 1933 On the stability for three-dimensional disturbances of vis- cous fluid flow between parallel walls. Proceedings of the Royal Society of London. Series A 142 (847), 621–628. Sutherland, B.R. 2010 Internal gravity waves . Cambridge University Press. Taylor, G.I. 1931 Effect of variation in density on the stability of superposed streams of fluid. Proc. R. Soc. Lond. A 132, 499–523. Trefethen, L.N., Trefethen, A.E., Reddy, S.C. & Driscoll, T.A. 1993 Hydrodynamic stability without eigenvalues. Science 261, 578–584. Waugh, D.W. & Randel, W.J. 1999 Climatology of arctic and antarctic polar vortices using elliptical diagnostics. J. Atmos. Sci. 56 (11), 1594–1613. Winters, K. B., MacKinnon, J. A. & Mills, B. 2004 A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Ocean. Tech. 21, 69–94. 129 Appendix Stokes’ Theorem applied on a vorticity interface Stokes’ theorem relates the surface integral of the curl of a vector (in our case this vector is velocity) field ~u over a surface A to the line integral of the vector field over its boundary δA: ∮ δA ~u.d~l = x A (∇× ~u) .d ~A (5.1) Holmboe (1962) used this theorem to relate the interfacial displacement ηi with the difference in velocity perturbation (u+i − u−i ) produced at a vorticity interface; see Eq. (2.12). This equation is referred to as “Eq. (3.2)” in his paper. However, the relevant steps required to derive this equation has not been provided. In order to understand how Eq. (2.12) is obtained, we first graphically describe the problem in Fig. 5.1. The background velocity is such that the flow is irrotational when z > zi, and has a constant vorticity, say S, when z ≤ zi. When the interface is disturbed by an infinitesimal displacement ηi (solid black curve in Fig. 5.1), the velocity field also changes slightly - the perturbation velocity in the upper layer (z > zi) becomes u + i and that in the lower layer (z ≤ zi) becomes u−i . 130 Appendix Figure 5.1: Schematic of a vorticity interface - left half shows unperturbed velocity field, while the right half depicts infinitesimal interfacial displacement. Let us consider a circuit A-B-C-D, see Fig. 5.1. Applying Stokes’ theorem, we obtain (u+i − u−i )4x = S.A (5.2) where A = ηi4x is the area of A-B-C-D, and S = ∇ × ~u is the vorticity in this area. Therefore we obtain u+i − u−i = Sηi (5.3) which is basically Eq. (2.12). 131 Appendix Normal mode form of Holmboe instability Both interfaces in the Holmboe profile (Eq. (2.67)) individually satisfy the kinematic condition: ∂η1 ∂t = ∂ ∂x ( e−αψ2 + 1− 2α 2α η1 ) (5.4) ∂η2 ∂t = ∂ ∂x ( ψ2 + e−α 2α η1 ) (5.5) where ψ2 is the stream function perturbation at the lower interface. This interface being a density interface also satisfies the dynamic condition: ∂ψ2 ∂x = J α ∂η2 ∂x (5.6) We assume the perturbations to be of normal-mode form: η1 = η̂1e iα(x−ct), η2 = η̂2e iα(x−ct), and ψ2 = ψ̂2eiα(x−ct). Here the wave speed c is generally complex. Defining ς̂ = [ ψ̂2 η̂2 η̂1 ]T , we obtain the following eigenvalue problem: (M + cI) ς̂ = 0 (5.7) where M =  0 J/α 0 1 0 e−α/(2α) e−α 0 (1− 2α)/(2α)  (5.8) 132 Appendix Eq. (5.7) generates the following characteristic polynomial: c3 + ( 1− 2α 2α ) c2 − J α c− J ( 1− 2α 2α2 ) + J e−2α 2α2 = 0 (5.9) This equation produces complex conjugate roots only when the discriminant is negative. Since the presence of complex roots signify normal-mode instability, negative values of the discriminant is of our interest. The discriminant (D) in this case is given by: D = 16α2J2 − αJ [8 (2α− 1)2 + 36e−2α (2α− 1) + 27e−4α] − (1− 2α)3 (2α− 1 + e−2α) (5.10) Imposing the condition D < 0, we find 1 2A ( −B − √ B2 − 4AC ) ≤ J ≤ 1 2A ( −B + √ B2 − 4AC ) (5.11) where A = 16α2 B = −α [8 (2α− 1)2 + 36 (2α− 1) e−2α + 27e−4α] C = ( 2α− 1 + e−2α) (2α− 1)3 Thus Holmboe instability occurs only when the condition in Eq. (5.11) is satisfied. 133