NON-SYMMETRIC HOLMBOE WAVES By Susan Patricia Haigh B. Math. (Pure & Applied Mathematics) University of Waterloo, 1988 M. A. Sc. (Aerospace) University of Toronto, 1990 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1995 © Susan Patricia Haigh, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of The University of British Columbia Vancouver, Canada Date DE-6 (2188) 3 Oc:- Abstract When two flows of different velocity and density meet, a shear layer with a density gradient is formed. Under certain conditions this flow can be unstable. A statically stable stratified shear flow in which the density interface is much thinner than the shear layer thickness can be linearly unstable to two modes with equal growth rates and equal and opposite phase speeds. The superposition of these two modes is called a Rolmboe instability. This instability is oniy possible when the flow is symmetric about the center of the shear layer. We examine the effect of breaking this symmetry by allowing the center of the density interface to be displaced with respect to the center of the shear layer. There are three major components to this study: linear stability analysis, nonlinear numerical simulations, and comparison with laboratory experiments. Linear stability analysis is used to examine the effect of the density interface offset on the overall stability of the flow. Both inviscid theory with piecewise linear background velocity and density profiles, and viscous theory with smooth background profiles are used. As in previous studies, it was found that the growth rate of one mode increases and that of the other mode decreases as the density interface displacement is increased. The precise behaviour depends on the relative thickness of the density interface with respect to the shear layer thickness. For inviscid theory with piecewise linear background profiles, we show that the initial perturbations must be two-dimensional. When the effects of viscosity and diffusion are included, however, it may be possible that the weaker mode is initially three-dimensional. Detailed analysis of the energy transfer in the linear regime indicates that when the background flow loses its symmetry, it is the mode with the larger growth rate that is primarily responsible for the extraction of energy from the mean flow. 11 Two-dimensional numerical simulations are used to examine the nonlinear develop ment of “non-symmetric” Holmboe instabilities. We start by pertiirbing the flow with the stronger mode predicted by linear theory. We then examine the response of the flow to weaker mode. Finally, we impose both modes. By comparing the development of these three flows we are able to study the interaction between the two modes and the effect of initial conditions on the development of instabilities. Although the initial development of instabilities depends on the initial conditions, this dependence weakens as the density interface offset is increased. Also, preliminary results indicate that long term behaviour of the perturbations are independent of initial conditions. Results of the numerical simulations are compared to both symmetric and nonsymmetric Holmboe instabilities that have been observed in laboratory experiments. Since we are unable to compute the flow at the high Prandtl number of the salt stratified experimental flows, we run the simulation for a range of increasing Prandtl numbers to determine how instabilities of thermally stratified flows (with low Prandtl number) and salt stratified flows differ. Results of these simulations indicate that differences between the experimental and numerical flows can be attributed to the thicker density interface and lower Prandtl number used in the numerical simulations. When the density interface is thicker and the Prandtl number lower, the waves or billows formed by the instabilities are not as sharply defined as in the experiments. 111 Table of Contents Abstract ii Table of Contents iv List of Tables vii List of Figures viii Acknowledgements xi 1 Introduction 1 2 Background 8 2.1 3 Governing Equations . . 8 2.1.1 Small Density Variations 2.1.2 Stream Function Representation 11 2.1.3 Nondimensional Equations 11 9 . . . 2.2 Linear Stability Analysis 13 2.3 Literature Review 16 Linear Stability Analysis 23 3.1 Derivation of Linear Equations 23 3.2 Piecewise Linear Profiles 24 3.2.1 Results: Piecewise Linear Profiles 26 3.2.2 Three-Dimensional Perturbations: Piecewise Linear Profiles iv . 28 3.3 4 Smooth Profiles 31 3.3.1 Numerical Method 32 3.3.2 Results: Smooth Profiles 3.3.3 Three-Dimensional Instabilities: Smooth Profiles 42 3.3.4 Eigenfunctions 43 Nonlinear Numerical Simulations 4.1 4.2 6 34 73 Numerical Method 4.1.1 5 . Numerical Stability Results 4.2.1 Methodology 4.2.2 Parameter Values 4.2.3 Right Mode 4.2.4 Left Mode 4.2.5 Both Modes Comparison with Laboratory Experiments 128 5.1 Symmetric Case 129 5.2 Non-Symmetric Case 131 Summary and Conclusions 145 6.1 Linear Stability Analysis 145 6.1.1 Inviscid Results 146 6.1.2 Viscous Results 146 6.2 Nonlinear Simulations 149 6.3 Comparison with Laboratory Experiments 151 6.4 Future Work 151 V List of Symbols 153 Bibliography 158 A Dispersion Relation for Piecewise Linear Profiles 164 B Maximum Growth Rates for Holmboe’s Equation 168 B.1 Ilolmboe Instabilities 169 B.2 Kelvin-Helmholtz Instabilities 170 C Energy Equations 172 C.1 Kinetic Energy 172 C.2 Potential Energy 176 D Period of Perturbation Kinetic Energy vi 178 List of Tables 3.1 Values of Prandtl number for water 36 4.1 Parameters used for the numerical simulations 84 5.1 Parameters of flows from various laboratory experiments vii 130 List of Figures 1.1 Schematic of a two layer flow with a thin density interface 6 1.2 Two-layer flow studied by Holmboe [23] 7 2.1 Two-layer flow considered by Helmholtz [22] and Kelvin [25] 22 3.1 Background flow for piecewise linear velocity and density profiles 51 3.2 Stability diagrams for inviscid case 52 3.3 Growth rates of most unstable modes for inviscid case 53 3.4 Phase speeds of most unstable modes for inviscid case 54 3.5 Linear stability diagrams for small J when 55 3.6 oc/’J for Kelvin-Helmholtz and Holmboe waves 56 3.7 c//T for e 57 3.8 Background flow for hyperbolic tangent velocity and density profiles 3.9 Maximum growth rate in Holmboe regime with varying R = = 0 0, 0.25, and 0.5 58 . 59 3.10 Linear stability diagrams for smooth profiles: varying R 60 3.11 Maximum growth rates for smooth profiles: varying R 61 3.12 Linear stability diagrams for smooth profiles: varying H 62 3.13 Maximum growth rates for smooth profiles: varying H 63 3.14 Smooth profiles: varying Re 64 3.15 Eigenfunctions and correlations when = 0 65 3.16 Eigenfunctions and correlations when e = 0.5 66 3.17 Evolution of perturbation kinetic energy and correlations with e 3.18 Evolution of perturbation kinetic energy and correlations with viii = 0. 0.25 . . . 67 69 3.19 Evolution of perturbation kinetic energy and correlations with = 0.5 71 4.1 Growth rates from linear stability analysis 100 4.2 Absolute stability region of fourth order Runge-Kutta scheme 101 4.3 Perturbation kinetic energy for right mode 102 4.4 Perturbation kinetic energy for left mode 4.5 Perturbation kinetic energy for both modes 4.6 Effect of grid size when e 4.7 Density and vorticity for right mode only with = 0 106 4.8 Density and vorticity for right mode only with e = 0.25 107 4.9 Density and vorticity for right mode only with e 0.5 108 = . . 103 . 104 0.5 105 4.10 Linear growth rates from nonlinear simulations 109 4.11 Positions of waves for right mode only 110 4.12 Average phase speeds from nonlinear simulations 111 4.13 Energy budgets for right mode oniy with e = 0 4.14 Energy budgets for right mode only with e = 0.5 . . . 112 . 113 4.15 Evolution of mean profiles for right mode only 114 4.16 Density and vorticity for left mode only with = 0.25 4.17 Density and vorticity for left mode only with e = 0.5 4.18 Positions of waves for left mode only . . . . 115 116 117 4.19 Density and vorticity for both modes with 0 118 4.20 Density and vorticity for both modes with e = 0.1 119 4.21 Density and vorticity for both modes with e = 0.25 120 4.22 Density and vorticity for both modes with = 0.5 121 4.23 Positions of waves for both modes 122 4.24 Energy budgets for both modes with ix = 0 123 4.25 Wave positions and energy budgets in linear regime for e = 0 . 4.26 Wave positions and energy budgets in nonlinear regime for e 4.27 Energy budgets for both modes with E = . 0 0.5 . . 124 125 126 4.28 Evolution of mean density profiles for both modes . . 127 . 5.1 Experiments of Zhu [72] 5.2 Results from numerical simulations for J 5.3 Experiments of Guez & Lawrence [20] 139 5.4 Experiments of Lawrence et al. [38] 140 5.5 Results from numerical simulations for J = 0.03 with R = 3 5.6 Results from numerical simulations for J = 0.03 with R = 6 5.7 Results from numerical simulations for J = 0.03 with R = 10 5.8 Results from numerical simulations for J = 0.06 137 B.1 Holmboe’s [23] stability diagram = 0.3 138 . . . . . . 141 142 143 144 171 x Acknowledgements This thesis would not have been completed without the constant guidance and encour agement from my supervisor, Dr. Gregory Lawrence. I am indebted to Dr. Brian Wetton for keeping me numerically honest (at least to order Ax ). In addition, I would like to 2 thank the other members of my supervisory committee for their helpful suggestions in the preparation of this document. The work in this thesis is based on previous work by Dr. Bill Smyth who was generous in his advice and who promptly responded to my pleas for help. I also greatly benefited from my acquaintances with Dr. Coim Caulfield, Dr. Craig Stevens, and Dr. Noboru Yonemitsu. Although I often found dealing with computers and their software frustrating, the following people make the experience less overwhelming: Bob Bajwa, Djun Kim, Dave Moulton, Peter Newbury, Tom Nicol, George Phillips, and Craig Stevens. There are many friends from both the Department of Mathematics and the Environ mental Fluid Mechanics group in the Department of Civil Engineering who made my stay at UBC enjoyable and the rain tolerable. I will miss you all. Finally, a special thanks to my husband James for putting up with my strange schedule as I worked towards completing this thesis and for making sure that I came home to a hot meal and a cold beer at the end of a long day. I received support from the Natural Sciences and Research Council of Canada and a Izaak Walton Killam Memorial Pre-Doctoral Fellowship from the University of British Columbia. xi Chapter 1 Introduction Mixing of stably stratified shear flows is an important aspect of many problems in oceanography, meteorology, and several branches of engineering. Properly quantifying this mixing, in order to accurately model turbulence, remains a problem of great im portance (see, for example, Gregg [19] and Fernando [14]). Studying mixing in a stably stratified flow has an advantage over other types of transition; stable stratification ap pears to sufficiently suppress the onset of turbulence to enable distinct events in the transition to turbulence to be classified (Sherman et al. [50]). Wave-like instabilities are the first to occur during transition. Thus understanding the types of wave-like in stabilities that can occur, and how they develop, is the first step in gaining complete understanding of mixing. The study of such instabilities remains one of the fundamental problems in the field of hydrodynamic stability. The traditional method used to determine the stability of a given flow is the normal mode approach to linear stability analysis. This method adds a perturbation (usually, but not necessarily, two-dimensional) to a background flow. The governing equations are then linearized about this state. The normal mode approach allows one to determine the wave number, growth rate, and spatial structure of the most unstable modes. A complete description of this method is given in section 2.2. The underlying assumption in linear stability analysis is that, if the perturbation from the mean flow initially consists of random noise, then it is the most unstable modes that grow fastest and dominate. All investigations into the stability of a flow start with determining its linear stability. 1 Chapter 1. Introduction 2 On its own, however, linear stability analysis does not provide an adequate portrait of the evolution of the flow, since it is only valid when the perturbation amplitudes are very small. It is well known that nonlinear effects rapidly become important and should not be ignored. With today’s computers it is possible to use numerical simulations to examine these nonlinear effects. One drawback, however, is that only specific cases can be examined. The results of linear stability analysis can be used as a guide, indicating which flow parameters may produce instabilities of interest. Also, since linear stability analysis tell us which modes are most likely to grow, the nonlinear simulation can be given a head start by using the results from linear stability analysis to initialize the perturbations. In this thesis, we restrict our attention to the stability of a statically stable stratified shear flow whose density interface is much thinner than the shear layer thickness, as depicted in figure 1.1. Such background flows occur in many physical situations, for ex ample, salt wedge flows (Yonemitsu [70]) and exchange flows (Zhu [72]). Yonemitsu [70] and Zhu [72] have conducted laboratory experiments to model such flows and have ob served wave-like instabilities at the density interface of the two fluids. Since the density gradient acts as a stabilizing force, the instabilities caused by the presence of the velocity shear do not always develop into fully turbulent flow. In this case, studying the onset and development of instabilities deals directly with many practical issues: for example, determining the interface location and the interfacial friction are of interest in many engineering problems. In addition, by determining the effect on the mean flow of small scale structures, a study of the development of interfacial instabilities can lead to a better understanding of large scale dynamics of two layer systems. We assume that the only non-zero component of the background velocity vector is in the x-direction and that the backgrollnd flow does not vary in this direction. In addition, we allow the background flow to develop in time only, since temporally developing flows Chapter 1. Introduction 3 are easier to compute than spatially developing flows. These assumptions result in the background flow having parallel streamlines, and hence the term parallel flow is used. A good first approximation to the background flow shown in figure 1.1 is the piecewise linear inviscid flow studied by flolmboe [23] (figure 1.2). This flow is characterized by a single parameter, the bulk Richardson number, J — 2(p _ 1 , 2 ) +-pi)(U pi)gh/(p U which is a measure of the relative strength of the stratification to the strength of the velocity shear. The flow depicted in figure 1.2 is linearly unstable to two types of instabilities. When J < 0.07, linear stability analysis predicts a Kelvin-Helmholtz instability. This instability is characterized by zero phase speed and is easily identified by the rolling up of the density interface. Since the original work of Helmholtz [22] and Kelvin [25], much has been accomplished towards understanding this instability, including the first nonlinear numerical simulations by Tanaka [61] and Patnaik et al. [44]. Recent research in this area includes the study of secondary instabilities (see, for example, Klaassen & Peltier [31] and Staquet [58]), and understanding how the Kelvin-llelmholtz instability develops into fully three-dimensional turbulence (Caulfield & Peltier [8]). The second type of instability for the flow depicted in figure 1.2 is a Holmboe insta bility. This is defined as the superposition of two unstable modes of equal growth rate and equal, but opposite, phase speeds, and is named after Holmboe [23]. A Holmboe instability is characterized by a right-moving wave whose energy is concentrated above the center of the shear layer and a left-moving wave whose energy is concentrated below the center of the shear layer. Although Holmboe instabilities occur for J > 0, they are the only possible instabilities when J > 0.07. These instabilities are of interest as they are the only mechanism through which the flow may become unstable when stratification is strong compared to the strength of the velocity shear (J > 0.07). flolmboe [23] studied the case where the flow is symmetric about the center of the shear layer. In nature, the background flow is not necessarily symmetric, e.g. the salt Chapter 1. Introduction 4 wedge flows described by Yonemitsu [70]. There are two ways that the flow may lose its symmetry: either by having horizontal boundaries placed at different distances from the center of the shear layer (Yonemitsu [70]), or by displacing the density interface with respect to the center of the shear layer (Lawrence et al. [38]). A more general case would allow both of these to occur. When the symmetry of the background flow is broken, one of the modes of the flolmboe instability dominates. Non-symmetric background flows have been used to explain why “one-sided” instabilities often occur in laboratory experiments, instead of the Holmboe instabilities predicted by the symmetric model (see, for example, Maxworthy & Browand [40], and Koop & Browand [33]). Since relatively little research has been done to understand these “non-symmetric” Holmboe instabilities, the main purpose of the present study is to understand how the displacement of the density interface affects the development of instabilities in a flow which, for the symmetric case, supports pure Holmboe instabilities. In chapter 2, we present the background needed for the present study, including an introduction to the governing equations, a description of linear stability analysis, and a review of the relevant literature. Chapter 3 focuses on how displacing the position of the density interface with respect to the center of the shear layer affects the linear stability of the flow. This is separated into two distinct parts: inviscid theory for piecewise linear background profiles, and viscous theory for smooth background profiles. The inviscid theory is an extension of the work of Lawrence et al. [38] who modified flolmboe’s original model by allowing a displacement of the density interface. In section 3.2.1 we examine the effect of imposing horizontal boundaries on this model. We also examine the possibility that the perturbations are initially three-dimensional (section 3.2.2). Since inviscid flows with piecewise linear profiles do not take into account such physical parameters as the Reynolds number and the Prandtl number, section 3.3 is devoted to the linear stability of viscous flows with smooth background profiles. In particular, we Chapter 1. Introduction 5 examine how varying the thickness of the density interface affects the stability of the flow when the displacement of the density interface is also varied. As with the inviscid case, we examine the effect of horizontal boundaries and the possibility of three-dimensional instabilities. We also examine, in detail, the shapes of the normal modes and how they govern the growth mechanism of the instability. In chapter 4, numerical simulations are used to examine the nonlinear evolution of non-symmetric Holmboe instabilities. Since a Holmboe instability consists of two waves propagating in opposite directions, all previous studies of the nonlinear evolution of Holmboe instabilities have initialized the perturbations with both modes predicted from linear stability analysis. As a starting point in our nonlinear investigations, we initialize the flow with one mode only. We start by perturbing the flow with only the right moving instability, because linear theory predicts that this mode dominates the flow as the offset increases. We then examine the response of the flow to left moving perturbations. Finally, we impose both a left moving and a right moving perturbation. By comparing the results of these different cases, we are able to determine how initial conditions affect the development of an instability and the interaction between the two modes. As a final test of the validity of the numerical simulations, chapter 5 concentrates on comparing the results of the nonlinear simulations with those of laboratory experiments. We start by comparing results of chapter 4 with the recent laboratory observations of Holmboe instabilities by Zhu [72]. Next we run simulations which match the flow con ditions of Lawrence et al. [38] and Guez & Lawrence [20]. These latter flows are clear examples of “one-sided” instabilities. Comparing results of the numerical simulations with both symmetric and non-symmetric instabilities observed in the laboratory helps us identify strengths and weaknesses of the present model. Chapter 1. Introduction z x Figure 1.1: Schematic of a two layer flow with a thin density interface 6 Chapter 1. Introduction 7 z 1 U 1 p 2 p Figure 1.2: Two-layer flow studied by Holmboe [23] x Chapter 2 Background 2.1 Governing Equations A two-dimensional, unsteady, stratified flow can be described completely in terms of the following variables: u and w (velocity components in the horizontal, x, and vertical, z, directions, respectively), p (density), 0 (concentration of stratifying agent, e.g. heat or salinity), and p (pressure). Assuming that the fluid is Newtonian, incompressible, and the equation of state linear, the governing equations are given by 1. conservation of momentum: PUt PWt + + PUUx UW + PWUz — (u + (u + —p + PWWz = Pz Pt + PUv + P’Wz — pg + wz)x), (w + (u + w)) (2.1) (2.2) 2. conservation of mass: + UPx + W,Oz 0, (2.3) 3. diffusion equation for stratifying agent: O + uO + wO = ,cA, (2.4) 4. the equation of state: p = po[l 8 — 7(0 — 0w)]. (2.5) Chapter 2. Background 9 In the above set of equations, g, t, i, and y are the gravitational acceleration, the molecular viscosity, the diffusivity of the stratifying agent, and the coefficient of expansion for the stratifying agent, respectively, and are all assumed to be constant. In addition, A represents the two-dimensional Laplacian operator, i.e. A = - + -. The equation of state (2.5) and the heat equation (2.4) are often good approximations when the fluid is a pure liquid, where the effects of compressibility are minor (see Pedlosky [45]). In (2.5), 0 is the density at a standard state &o. The validity of the two-dimensional assumption p will be examined in chapter 3. If p and p are expanded about a state of hydrostatic equilibrium Pa and Pa, as dis cussed, for example, in Turner [65], the equations of conservation of momentum can be written as PUt PWt + + PUUx PUWx + + PWUz = PWWz = where p = pa(Z) 2.1.1 Small Density Variations + p(x,z,t), P —p + —p — (Au + (u + p’g + i = Pa(Z) +p(x,z,t), wz)), (Aw+(u + w)), and -a-— (2.6) (2.7) —gp. For the flows we consider, the maximum density variation Ap, is small compared to the standard density p, i.e. Ap/po = 0(10—2), or less. Using this property, the above equations can be simplified considerably. For a fluid in which the effects of compressibility are minor, changes in density have little effect on the conservation of mass, and (2.3) can be approximated by the continuity equation for an incompressible and nondiffusive flow. u, + w 0. (2.8) In order to derive the above approximation, we use a scaling argument. For the mean background flow depicted in figure 1.1, motions in the x-direction are much larger than Chapter 2. Background 10 those in the z-direction. We can write u = 0(U) and w we let x 0(W), where W << U. Also, 0(L) and z = 0(D). We expect the time scale in the vertical direction to be of the same order of magnitude as the time scale in the horizontal direction. Hence 0(L/U) = 0(D/W), which implies that D <<L. Finally, since variations in density are small, the equation of state (2.5) suggests that we may write p = po(i + r(x, z, where p is constant and r = 0(e) < 0(10—2) << 1. Using the equation of state, we can write (2.4) as Pt 0 (PoU) + ‘UPx 0 (PoU) + WPz 0 + = (PoW) 0 () ‘Pzz 0 () (2.9) We notice that all terms on the left hand side of (2.9) have the same order of magnitude. The dominant term on the right hand side is lcpzz, since D << L. As diffusive terms on the right hand side cannot dominate, we have 0 (P). Subtracting (2.9) from the conservation of mass (2.3), we obtain pu o() + P’Wz oQ°°) = Pxx o() — Pzz o() (2.10) Now, (2.11) The relations in (2.11) indicate that we can neglect the terms on the right hand side of (2.10) and obtain (2.8). We can now simplify the momentum equations (2.6) and (2.7) using the Boussinesq approximation. This approximation assumes that variations in density are only important in the gravitational term of (2.7). This assumption is valid providing the maximum density variation Lp is small compared to po, as is the case here (see Turner [65] for a more detailed discussion of the Boussinesq approximation and the types of flows for which Chapter 2. Background 11 it is valid). Using the Bonssinesq approximation and the incompressibility condition (2.8) derived above, we can rewrite (2.6) and (2.7) as t 2 + uu + wu = — + vAu, (2.12) P0 Wt where v = + — + = — yAw, (2.13) ,u/po is the kinematic viscosity. Equations (2.12) and (2.13) along with (2.8) and (2.9) form or new set of governing equations. These equations are valid when the fluid studied is a pure liquid for which density variations are small and the effects of compressibility are negligible. Stream Function Representation 2.1.2 We can use the continuity equation for incompressible flow (2.8) to define a stream such that function thb u=— Taking 8(2.12)/8z — L’ 1 th ; w=———. (2.14) 8(2.13)/8x and using (2.14) one obtains the following equation for the stream function. (A + u(A) + w(A) = + (2.15) We have reduced the original system of five equations, (2.1) to (2.5) with five unknowns, ‘u, to, p, p, and 6, to a system of two equations (2.15) and (2.9), with two unknowns, ‘(x,z,t) 2.1.3 and p(x,z,t) = pa(z) + p’(x,z,t). Nondimerisional Equations It is convenient to examine the nondimensional equations. This has the advantage of classifying various flows by the means of certain nondimensional parameters. A schematic Chapter 2. Background 12 of the background flow that we will be examining is depicted in figure 1.1. Precise definitions of the background profiles are discussed in chapter 3. Defining the mean velocity as U )/2 and using the length scale L, the velocity scale 6U, and 2 1 + U (U the density scale 6 p, where L, 6U, and Sp will be defined as needed, the nondimensional variables, denoted by ‘, are defined as = 6Uu* + U b = 6ULb* x = Lx* = 6Uw* = pp + 6 po Z = Lz* p (2.16) L* — su Substituting the above definitions into the governing equations (2.15) and (2.9), we obtain (A) + u(A) + w() Jp + (2.17) RPAP. (2.18) = and Pt + UPx + WPz = In the above equation the * notation for nondimensional quantities has been dropped. The three nondimensional parameters, namely, the bulk Richardson number J, the Reynolds number Re, and the Prandtl number Pr, are defined as 6 gL p po(6U) The Prandtl number, Pr = v/ii, Re , = SU L v , Pr = . (2.19) is the ratio of kinematic viscosity to thermal diffusiv ity, whereas the Schmidt number, Sc = v/i, is the ratio of kinematic viscosity to the coefficient of mass diffusivity (we are interested in the diffusion of salt). Throughout the remainder of this study, however, we shall refer to the ratio regardless of the stratifying agent (heat or salt). v/it as the Prandtl number Chapter 2. Background 2.2 13 Linear Stability Analysis The normal mode approach to predicting linear stability is the most common analytic tool used to investigate the stability of a flow. For a parallel shear flow, we start with a steady flow in the x-direction which satisfies the governing equations. Perturbations are then superimposed onto the steady flow and linearized equations for the perturbations are derived. The normal mode approach assumes that the perturbations are proportional to exp[ic(x — ct)]. For example, for the governing equations (2.15) and (2.9) we write = (z) + ‘(x, z, t), (2.20) = Pa(Z) + p’(x, z, t), (2.21) {(z)ei_ct)}, (2.22) p’(x,z,t) = {(z)ei_ct)}. (2.23) z, t) p(x, Z, where z, t) = Here the nondimensional wave number, c, is real and positive, and c = Cr + ic is the complex nondimensional wave speed. It is the amplification factor cc which determines the linear stability of the problem. If c > 0, the flow is linearly unstable since the perturbations will grow exponentially. Conversely for c < 0 it is stable. If c = 0 then the flow is said to be neutrally stable. Choosing c real and c complex means that the solution will grow in time (known as the temporal modes). This is in contrast to the physical case where solutions tend to grow in space (known as the spatial modes). Spatial modes can be found by choosing c real and o complex (see Drazin Reid [12] for a discussion on the use of spatial modes). We will restrict our study to that of temporal modes, however, since the nonlinear evolution of these is more easily computed. Gaster [17] provides the following approximate transformation between the parameters of the temporal problem Chapter 2. Background 14 and those of the corresponding spatial problem. If the temporal solution has parameters a = a = ar and ac ar — = iarcj/cg a,cr + and iarcj, then the corresponding spatial solution has parameters = arcr, where Cg = ã(arcr)/ãar is the group velocity. This approximation is valid providing amplification rates are small. When linear stability analysis was first used, many investigators hoped that transition to turbulence would follow directly from linear instability. It was quickly shown, however, that this is not necessarily the case. Linear analysis is only valid for a very short time before nonlinear effects become important. Nevertheless, linear instability does correctly describe the onset and early evolution of infinitesimal perturbations. It also seems to give a qualitatively correct indication of the overall stability of the flow (Maslowe [39]). For this reason, much of the literature has been devoted to linear stability theory including books by Chandrasekhar [9] and more recently by Drazin & Reid [12], as well as review articles by Maslowe [39] for viscous shear flows, Pellacani [46] for stratified shear flows, and Gage [16] for viscous stratified shear flows. The normal-mode approach of linear stability analysis reduces the system of nonlin ear partial differential equations into a system of linear ordinary differential equations where the boundary conditions depend on the flow configuration. These new differential equations describe an eigen-problem where c is the eigenvalue. The coefficients of the resulting differential equations are, in general, not constant. For this reason these equa tions can only be solved exactly for a handful of special cases. These are mostly limited to steady flows approximated by either piecewise constant or piecewise linear velocity and density profiles. Although these flows are not physically realistic, they often admit simple solutions. For more complicated flow profiles it is usually necessary to use either asymptotic approximations, or to solve the equations numerically. Having described the linear stability problem, we will now discuss how it is applied. For a given flow configuration, the wave speed c depends on the value of the wave number Chapter 2. Background 15 a, which can take on any real value. Depending on the problem, there is, for each value of a, either a continuous spectrum of wave speeds c, a collntable infinite number of discrete values of c, a finite number of discrete values of c, or some combination of these cases (see Drazin & Reid [12]). The complete solution for the stream function perturbation b’ is of the form z, t) where c is the flth J A(z)ei_t)da (2.24) wave speed associated with the wavenumber a and A are coefficients which are determined by the boundary conditions. If there is a continuous spectrum of wave speeds, then the summation in (2.24) is replaced by an integration. Obviously for the solution to be stable, we need E{c} < 0 for all a and n whereas for instability we only need {c} > 0 for one value of a and n. In view of the above form of the general solution, it is easy to see that solving the complete linear problem is not a trivial task. Fortullately, one does not usually need to determine “ completely to gain some understanding of the stability mechanism. In fact, for the most part, the use of linear analysis is restricted to determining the dependence of the wave speed c on the wavenumber a. This dependence is usually described in the form of a dispersion relation. Once this relation has been obtained, either exactly, with asymptotics, or numerically, one of two problems is usually addressed. The first deals with transition from laminar to turbulent flow and is concerned with determining at what value of a given parameter, say the Reynolds number, does the flow first become linearly unstable. This value is often called a critical value. The second problem deals with a flow which is known to be unstable. In this case, given a fixed value of a parameter, for example, the Reynolds number, we wish to know at what wave number a does the most unstable mode occur (i.e., the largest value of act). Both approaches are used in the present study. Chapter 2. Background 2.3 16 Literature Review In this section we review the relevant literature. We concentrate on important results that hold for stratified shear flows in general and those that are specific to Holmboe instabilities. Many of the major results obtained in the study of the stability of stratified shear flows have been for inviscid flows. An advantage of performing linear stability analysis on inviscid flows is that exact solutions to the dispersion relation can often be found. The simplest model of a stably stratified shear flow is depicted in figure 2.1. Instability in such a flow was first remarked upon by Helmholtz [22] and calculated by Kelvin [25]. Details of the calculations are given in section 232 of Lamb [36] and section 4 of Drazin & Reid [12]. g(p 1 U — Results of linear stability analysis indicate that a mode is unstable if (Ui 2 p) < kpip — , where k is the wave number of the normal mode. When ) 2 U 2 we can always find a value of k which satisfies this condition, and hence the U flow is unstable. If, in addition, P1 = P2, then the flow is unstable at all wave numbers. This type of instability is called a Kelvin-Helmholtz instability and is characterized by a ‘rolling up’ of the shear layer into periodic train of symmetric vortices which travel with the mean flow (Smyth [51]). In 1962, Holmboe [23] examined the stability of a discontinuous density interface contained within a piecewise linear velocity shear of finite thickness (figure 1.2). The stability of such a flow is governed by the bulk Richardson number J. When J < 0.07, the velocity shear dominates and Kelvin-Helmholtz instabilities are predicted. He discovered, however, that when J > 0.07, the flow is unstable to two modes having equal growth rates and equal, but opposite, phase speeds. Using the method of symmetric waves he was able to determine the time dependence of the phase speeds for the two modes. The superposition of these two modes is usually referred to as a Holmboe instability, a term Chapter 2. Background 17 first used by Browand & Wang [4]. Smyth [52] offers the following physical explanation for the difference between a Holm boe instability and a Kelvin-Helmholtz instability. If one thinks of the density interface as a flexible boundary between the two fluids, then the boundary is almost rigid when stratification is strong (i.e., large J). In this case, the two modes of the Holmboe instabil ity, i.e., the right moving mode located in the upper layer and the left moving mode in the lower layer, propagate through the flow with very little interaction. When stratification is weak (small J), however, the interface is very flexible and disturbances in the upper layer interact with those in the lower layer as they pass through each other. Below some critical value of J (J < 0.07 for the case studied by Holmboe [23]), the phase speeds of the right and left moving modes become zero and there is a phase locking of the two modes. The two modes then wrap around each other, forming the Kelvin-Helmholtz instability. Hazel [21] developed a numerical method for solving the linear stability problem for continuous profiles. He examined the stability of flows where the density interface thick ness has a finite value and is characterized by R, the ratio of the shear layer thickness to the density interface thickness. When R Helmholtz instabilities. When R = = 1, the only possible instabilities are Kelvin 5, unstable modes with non-zero phase speed, i.e., Holmboe instabilities, were found. This indicates that there is some critical value of R below which Holmboe instabilities do not exist. By numerically computing stability diagrams for various values of R, Smyth [51] found the critical value of R to be approxi mately 2.4. Smyth [51] and Smyth et al. [53] used numerical simulations to examine the nonlinear evolution of Holmboe instabilities and were able to confirm Holmboe’s predic tion of varying phase speed as the right and left modes interact. In addition, Smyth and Peltier have conducted two studies on the transition from Kelvin-Helmholtz to Holmboe instabilities. The first concentrated on linear theory [54], while the second examined the results of nonlinear numerical simulations [56]. Chapter 2. Background 18 Observations of both Kelvin-Helmholtz and Holmboe instabilities have been docu mented. Perhaps the best known examples of Kelvin-Helmholtz instabilities are the billow clouds near Denver, Colorado (Colson [11]) and the tilting tube experiments by Thorpe [63]. Holmboe instabilities are, in general, less well-known than Kelvin Helmholtz instabilities. There are numerous examples, however, of laboratory obser vations of Holmboe instabilities. These illclude the experimental results of Browand & Winant [5], Koop [32] (reproduced in part in Tritton & Davies [64] figure 8.8), and, more recently, Poiiliquen et al. [48], and Zhu [72]. Browand & Wang [4] measured the position of the stability boundary experimentally and found good agreement with the predictions of Holmboe [23]. The calculated growth rates for their low Reynolds num ber flows (Re = 1 (U — )h/v 2 U = 40—300), however, were an order of magnitude lower. There are two possible explanations for this discrepancy. The first lies in the piecewise linear profiles used by Holmboe [23]. Smyth [51] has computed the stability of an in viscid flow with smooth background profiles. He observed that as the thickness of the density interface increases, the maximum growth rate in the Holmboe regime decreases. Secondly, Nishida & Yoshida [43] have computed the stability boundaries for a two-layer stratified shear flow at various Reynolds numbers. They show that the region of insta bility decreases significantly with decreasing Reynolds number, indicating a decrease in the maximum growth rates. These results suggest that comparison with a model that includes the effects of viscosity and has smooth background profiles might yield better agreement between experiments and theory. In addition to Kelvin-Helmholtz and Holmboe instabilities, a third type of instability has been observed in stratified shear flows with thin density interfaces. This instability occurs for larger Richardson numbers when Holmboe waves are expected. Instead, the instability exhibits a “one-sidedness” with the wave protruding into one layer only and is believed to be either, the result of a shift of the density interface with respect to the Chapter 2. Background 19 center of the shear layer, or, the result of horizontal boundaries being non-symmetrically placed with respect to the centers of the shear layer and density interface. This has been observed by Keiilegan [26], Yoshida [71], Koop & Browand [33] (see also Maxworthy & Browand [40]), Lawrence et al. [38], and Yonemitsu [70]. In order to better understand this “one-sideness”, Lawrence et al. [38] modified Holmboe’s model by allowing the den sity interface to be displaced with respect to the center of the shear layer. They found that, as the density interface is displaced, one of the two modes of a Holmboe insta bility dominates (we will refer to this as a non-symmetric Holmboe instability). For a viscous flow, Yonemitsu [70] examined the effect of non-symmetrically placed horizontal boundaries. Similar results were obtained as discussed below. Smyth [51] made some preliminary computations on the nonlinear evolution of non-symmetric Holmboe waves. He found that these waves evolve, after many oscillation periods, into weakly nonlinear interfacial waves which are qualitatively similar to waves that have been observed in the lower atmosphere by Gossard & Richter [18] and Emmanuel et al. [13]. Inviscid theory has been used to determine several general criteria for instability. In 1880, Rayleigh [49] proved his famous inflection-point theorem: a necessary condition for instability is that the background velocity profile have an inflection point. Howard [24] proved the semi-circle theorem which states that, for any unstable mode, the complex wave speed, c, must lie inside a semi-circle in the upper half-plane whose diameter is given by the maximum amplitude of the background velocity. Miles [41] and Howard [24] showed that a sufficient condition for stability in a stratified parallel shear flow is that the gradient Richardson number Ri —(gpa(z)/(paUz(z) ) 2 , where U(z) is the horizontal background velocity profile, be greater than 1/4 everywhere in the fluid. It is possible, however, to have unstable flows for Ri diffusion are included (see Gage [16]). > 1/4 when the effects of viscosity and thermal Also, for viscous, thermally dissipative flows, unstable waves whose phase speed lie outside the semi-circle predicted by Howard [24], Chapter 2. Background 20 for inviscid flows, have been found (Gage [16]). Since boundaries are always present in the physical world, it is of interest to know how they affect the stability of a flow. Hazel [21] studied the effects of moving the boundaries in from infinity on an inviscid flow with hyperbolic tangent background velocity and den sity profiles. Two effects were observed. When the total domain height is greater than or equal to ten times the thickness of the shear layer, longer wavelengths are destabilized and are the only wavelengths affected by the presence of the horizontal boundaries. This is reasonable as we only expect waves of lengths comparable to the domain height to be affected. In addition, as the domain height is decreases further, shorter wavelengths become more stable. If the boundaries are moved in closer than some critical distance apart, the latter effect dominates and the flow becomes stable for all wavenumbers and Richardson numbers. Smyth [51] studied the effect of imposing boundaries on the profile studied by Holmboe (figure 1.2). Transition from Kelvin-Helmholtz to Holmboe instabil ities occurs at a lower value of J when boundaries are present. When the boundaries are not symmetrically placed, pure Holmboe instabilities are not found. Instead one mode of the Holmboe instability has a larger growth rate. For a viscous flow, Yonemitsu [70] examined how the presence of a lower boundary affects Holmboe instabilities. The wave that protrudes into the upper layer is little affected by the presence of the boundary unless the boundary is very close to the shear layer. The wave that protrudes into the lower layer, however, is significantly altered by the presence of the boundary indicating that perhaps the two modes are independent. In section 2.2 linear theory was developed under the assumption that perturbations are two-dimensional. At some point in the transition from laminar to turbulent flow the perturbations do indeed become three-dimensional, but one must consider whether they are initially two-dimensional during the time when linear theory is valid. In 1933, Squire [57] showed that, for an unstratified flow, a three-dimensional disturbance at a Chapter 2. Background 21 given Reynolds number is equivalent to the problem of a two-dimensional disturbance at a lower Reynolds number. Since a flow becomes more stable with decreasing values of the Reynolds number, this implies that all instabilities are initially two-dimensional for a homogeneous flow. Koppel [35] derived a generalized version of Squire’s result, showing that a three-dimensional disturbance at given Prandtl, Reynolds and Richardson num bers is equivalent to a two-dimensional disturbance at the same Prandtl number, a lower Reynolds number and a higher Richardson number. Thus any three-dimensional prob lem can be transformed into a two-dimensional problem. Since, for Holmboe’s instability, there is a region where the flow is destabilized as the Richardson number increases, it may be possible that the most unstable modes are three-dimensional. This was first suggested by Browand & Wang [4]. Two approaches have been used to examine this pos sibility. Smyth & Peltier [55] solved the three-dimensional problem and showed that, for a class of dissipative, stratified, parallel shear flows, the perturbations evolve directly into three-dimensional flows without going through the two-dimensional stage. Caulfield [7] calculated maximum growth rates to show that primary three-dimensional instabilities are possible in a generalized version of the inviscid three-layer flow studied by Taylor [62]. In chapter 3 we use the approach used by Caulfield [7] to examine the possibility of the initial perturbations being three-dimensional for the inviscid piecewise linear background flow studied by Lawrence et al. [38]. Chapter 2. Background 22 z 1 U 1 p n 2 r 2 U Figure 2.1: Two-layer flow considered by Helmholtz [22] and Kelvin [25] Chapter 3 Linear Stability Analysis 3.1 Derivation of Linear Equations We wish to examine the linear stability of the flow governed by (2.17) and (2.18). As described in Section 2.2, we split the flow field into a parallel mean component and a perturbation field, i.e., z, t) p(x, z, t) where ‘ = (z, t) + b’(x, z, t), (3.1) (z, t) + p’(x, z, t), (3.2) and p’ are O() <<1. Substituting these expressions into the governing equations (2.17) and (2.18), and collecting terms of equal magnitude, we obtain equations for both the background flow and the perturbations. 1 Ut = (3.3) Pt RP er Pzz, and (‘) + u(’) + w’u = Jp + (3.4) p + Üp + where u(x, z, t) = WZ RePr’ ü(z, t) + u’(x, z, t). Equations (3.3) indicate that the mean velocity and density profiles diffuse over time due to viscosity and the diffusivity of the stratifying agent. For the purpose of solving 23 Chapter 3. Linear Stability Analysis 24 the linear problem (3.4), however, we will use the initial mean profiles U(z) and Pa(Z) = = ü(z, 0) ö(z, 0). Hence st” and p’ are the perturbations from the initial mean state. Using the method of normal modes, we write ‘(x,z,t) = {(z)ci_ct)}, (3.5) {(z)ei_ct)}. (3.6) z, t) Substituting these into the linear equations (3.4) we obtain the eigen-problem: (U — — ) 2 a — Uzz (U_c)7pa = J — aRe — —z (a2) aRePr 2a + (37 Notice that the above set of equations reduce to the Orr-Sommerfeld equation when no stratification is present. Furthermore, the Taylor-Goldstein equation is obtained if there is no diffusion. Similar equations have been derived by Koppel [35], and Baldwin & Roberts [1]. 3.2 Piecewise Linear Profiles The starting point of our investigation is to examine the stability of an inviscid flow with two layers of constant density and a layer of constant horizontal shear (figure 3.1). This background flow has the advantage of admitting an exact dispersion relation to the normal mode analysis. Also, other than the density interface displacement d/h, the only additional parameter is the bulk Richardson number J. Finally, results of inviscid theory with piecewise linear profiles often compare favourably with experimental results. For the reasons given above, an inviscid flow with piecewise linear background profiles is an obvious first approximation to our problem. Chapter 3. Linear Stability Analysis 25 We define the background velocity, U= (u, ) = (U, 0), and density profiles as 1 U h/2<z<H/2, U1hU2z+UU2 —h/2<z<h/2, —H/2 < z < —h/2, (12 Pa = 1 P1 z>—d, 1P2 z<—d, (3.8) where h is the shear layer thickness, H is the total depth of the flow, and d is the displacement of the density interface from the center of the shear layer (figure 3.1). If we nondimensionalize using L and po = = h/2, SU 1 (U — U2)/2, U = U ) /2, p 1+2 (U 5 (p2 — p1)!2, (p2 + P1)! 2 then the nondimensional background profiles are U 1 1<z<H/2, z —1<z<1, —1 —H/2<z<—1, (3.9) 1—i Pa = z<—e, where H is now the nondimensional depth of the flow, and e = 2d/h. Since the flow is inviscid, (3.3) indicates that there is no diffusion of the mean flow. Furthermore, (3.7) reduces to the nondimensional Taylor-Goldstein equation: (U — c)( — a2) — + 2J6(z + E) = U 0, (3.10) where 6(z) is the delta function. Assuming the rigid wall boundary condition techniques outlined in Drazin 0 at z = +H/2 and using the Reid [12] (see appendix A for details of the calculations), Chapter 3. Linear Stability Analysis 26 we obtain the following dispersion relation. 4 1 + c 2 c 3 +a a where a = 0, = a(e, c, J, H) and are given in appendix A. As H (3.11) —* oc we retrieve the results of Lawrence et al. [38]. Also, this is simply a special case of the dispersion relation derived by Smyth [51] for piecewise linear profiles with non-symmetric horizontal boundaries. 3.2.1 Results: Piecewise Linear Profiles Symmetric Case (s For e 0 and H = 0) oo, we reproduce the results of Holmboe [23] (figure 3.2). When J > 0.07, Rolmboe instabilities are the only instabilities that occur. In the unstable region, there are two modes of equal growth rates and equal but opposite phase speeds (figures 3.3 and 3.4). The growth rate of the Holmboe instability increases with increasing J until it reaches a maximum value (figure 3.3) after which it decreases as J increases. When J < 0.07, the flow is unstable to both Kelvin-Helmholtz and Holmboe instabilities. An enlargement of the stability diagram in this region is shown in figure 3.5. When s 0, (3.11) reduces to 4 + c 2 = 0 a Solutions to (3.12) are oscillatory, i.e., Cr 0, when a (3.12) — four solutions, two stable and two unstable, of the form c 4 < 0. In this case there are 4a = +Cr of the two unstable modes is Holmboe’s instability. When a there are two stationary (Cr = + ic. The superposition — 4 > 0 and a 4a 2 > 0, 0) unstable solutions to (3.12) with different growth rates. These are Kelvin-Helmholtz instabilities. The contour a — 4 4a = 0 gives the stability boundary for the Holmboe instabilities. Results in figures 3.2 and 3.3 show the growth rates for the dominant Kelvin-Helmholtz Chapter 3. Linear Stability Analysis 27 instability only. It is clear from figure 3.5 that the Kelvin-Helmholtz instability domi nates for J close to 0 whereas the Holmboe instability dominates when J is near 0.07. It is not obvious at which value of the bulk Richardson number, JT, the flow changes from having a Kelvin-Helmholtz instability as the most unstable mode to a Holmboe instability. To determine the value of JT, the maximum growth rate for a fixed value of J must be found for both Kelvin-Helmholtz and Holmboe regimes. In appendix B the ci—J curves along which the maximum growth rates occur are calculated analytically for each region. These results are used to determine the maximum growth rates in each regime (figure 3.6). It is evident that transition from Kelvin-Helmholtz to Holmboe instabilities occurs at JT 0.046. By varying the value of H, we examine the effect of horizontal boundaries on the sta bility of the flow (figure 3.2). Decreasing H reduces the growth rate of the most amplified mode (figure 3.3). Moving the horizontal boundaries closer together, however, also causes the flow to become unstable at small values of c. This mechanism was initially observed by Hazel [21]. The presence of horizontal boundaries does not significantly affect the sta bility of the flow, however, until H 10. Since Kelvin-Helmholtz instabilities occur at smaller wave numbers than Holmboe instabilities, they are stabilized more rapidly as H decreases. This results in the transition from Kelvin-Helmholtz to Holmboe instabilities being shifted to a smaller value of J as the distance between boundaries decreases. When H 4, the flow is unstable to Holmboe instabilities only. As the boundaries continue to move together, the flow becomes stable at small J (see figure 3.2, H = 3). These results are consistent with those of Smyth [51]. Non-Symmetric Case ( 0) When e > 0 (figures 3.3, and 3.2), the flow has no pure Kelvin-Helmholtz or Holmboe instabilities. When H = oc, there are two unstable modes, one moving to the right and Chapter 3. Linear Stability Analysis 28 one moving to the left. As e increases, the growth rate of the right mode increases while its phase speed decreases. The left mode exhibits the opposite behaviour. Lawrence et al. [38] showed that, for fixed a and J, 4i 1) where (r) and (1) correspond to right and left moving waves, respectively. For a fixed value of J, the range of a for which the flow is unstable decreases for the right moving wave and increases for the left moving wave as increases. On the stability diagram, the limb corresponding to the right modes tilts towards the J-axis and thus the most unstable right mode for a fixed value of J occurs at smaller values of a with increasing s. The limb corresponding to the left mode tilts towards the a-axis and its most unstable mode occurs at larger a with increasing (figure 3.2). This is the behaviour described by Lawrence et al. [38]. We now examine the effect of horizontal boundaries (i.e., a finite value of H) on the stability of non-symmetric flows. To date, this effect has not yet been studied in any detail (Smyth [51] studied the case corresponding to e and H = 8). As with the = = 0.05 and a = 0.2 with H = 00 0 case, the presence of horizontal boundaries causes the flow to become unstable at small wave numbers. The boundaries have little effect on the fastest growing mode until H drops below 10, at which point growth rates of the fastest growing right and left modes decrease with decreasing H. Since right moving instabilities occur at smaller wave numbers than left moving instabilities, they are stabilized more quickly by the presence of boundaries. The result is that growth rates of right moving waves decrease more rapidly than those of left moving waves. This effect becomes even more pronounced as E increases and causes the left mode to dominate when H 5 (figure 3.3). 3.2.2 Three-Dimensional Perturbations: Piecewise Linear Profiles Until now, we have assumed that the perturbations are two-dimensional. Here, we discuss the relationship between two-dimensional and three-dimensional linear instabilities. Chapter 3. Linear Stability Analysis 29 The most general form of a three-dimensional disturbance (u’,p’,p’), where u’ (it’, = v’, w’), is that of a traveling wave whose amplitude varies with z and which propagates at an angle 9 t with respect to the x-axis (see White [68]), i.e., (u’, p’, p’) {(ü(z), (z), (z))ei Starting with a background flow of the form u Pa(Z) c0s+ sin-ct) } (3.13) (U(z), 0, 0) with density distribution in hydrostatic eqililibrium with the pressure pa(Z), we superimpose perturbations of the form (3.13). The equations of motions for an inviscid flow with the Boussinesq approximation under the linear assumption reduce to (u — c cost ) — ) 2 a — Paz1 — cos2(U_ cos ) = 0, (3.14) which, for the background profile described by (3.9), becomes (u — C cos ) — 2) — + 2J6(z + e) cos29(U_ C cos ) = (3.15) The above is simply the Taylor-Goldstein equation for two-dimensional perturbations, (3.10), with bulk Richardson number J/ cos 2 and complex phase speed c/ cos 9. Finding the growth rate for the three-dimensional problem is equivalent to finding the growth rate of a related two-dimensional problem and multiplying it by cos , i.e., ac(U, 9, a, J) = cos acj(U, 0, a, J/ cos 2 9). (3.16) This is simply a special case of Koppel’s [35] result for the more general case of a stratified shear flow with viscosity and dissipation. Koppel’s result is given in a form similar to the above by Smyth & Peltier [55]. The assumption that perturbations are initially two-dimensional is valid providing the fastest growing two-dimensional mode is more unstable than the fastest growing Chapter 3. Linear Stability Analysis 30 three-dimensional mode, i.e., a*c(U0*J/ t 9) < a*ci(U0c*J) cos9 2 where (3.17) is the wave number at which the maximum growth rate occurs for a given value of the bulk Richardson number. For J > 0, this is equivalent to having cosi9 a*c.(U , 0, J/ cos 2) ac(U, 0, a, J). (3.18) In order to validate the two-dimensional assumption, we must show that a*cj/ de creases with increasing J along the curve of maximum growth rate for a given J. In appendix B we derive analytic expressions for c—J curves for the case originally studied by Holmboe (H oc and E 0). These are used to find the maximum growth rates in each region (figure 3.6). In the Kelvin-Helmholtz regime, the maximum growth rate decreases with increasing J and thus (3.18) is automatically satisfied. In the Holmboe regime, however, there is a region where the maximum growth rate increases for increas ing J. Thus it may be possible that (3.18) is violated in this region. Figure 3.6 shows, however, that a’c//T does indeed decrease with increasing J. This implies that the initial instabilities are two-dimensional. For the case studied by Lawrence et al. [38] (H = oc, e > 0), we were not able to find exact expressions for the a*_J curves along which the growth rate is maximized. We can, nevertheless, use the results used to produce the stability diagrams (figure 3.2) to determine these curves for right and left moving modes. These are shown in figure 3.7 as are the curves a*cj/../J for e 0, 0.25, and 0.5. For larger values of , the maximum growth rate of the right mode decreases with increasing J and thus (3.18) is automatically satisfied. For smaller values of e, however, and for the left mode, there are regions where the maximum growth rate increases with J. Nevertheless, c*c//J decreases with increasing J and thus we expect the initial instabilities to be two-dimensional. Chapter 3. Linear Stability Analysis 31 Although we have shown that, for the inviscid case, the two-dimensional assumption is valid, these results should be used with caution. Smyth & Peltier [55] have shown that it is possible for the most unstable initial instabilities to be three-dimensional when viscosity and thermal dissipation are included. We will discuss this possibility further in section 3.3.2. 3.3 Smooth Profiles Although results for the piecewise linear profile case often agree with experimental results, two parameters are neglected: the Reynolds number and the Prandtl number. Further more, initially smooth background profiles are required in the nonlinear numerical sim ulations used in chapter 4. We now consider background flows with smooth profiles and include the effects of viscosity and diffusion. The initial mean velocity, (a, a) = (U, 0), and density profiles are depicted in figure 3.8 and are (i(z) = 2 2z 2 + 1 U U tanh—+ h 2 (3.19) pa(z) = _P2_P1tanh(z+d)+P1+P2, where h is the shear layer thickness and is the density interface thickness. Wang [66] and Yonemitsu [70] have found that the hyperbolic tangent function is a good approximation to the mean velocity fields measured in their laboratory experiments of salt stratified flows with thin density interfaces. Also, for large values of R = h/i’, the above density profile is a good representation of thin density interfaces. If we nondimensionalize (3.19) using L = h/2, SU = 1 (U — )/2, U 2 U = )/2, 2 1 + U (U Chapter 3. Linear Stability Analysis = (p2 — , and po 2 p’)/ = (p2 32 + pi)/ , as described in section 2.1.3, then the nondimen 2 sional background profiles are U(z) = tanhz, Pa(Z) = —tanhR(z+e), (3.20) where, as before, J is the nondimensional distance of the density interface displacement. With these background profiles, (3.7) become (U — c)( — (U ) 2 a Uzz — Paz — — = = J ), 4 +a — aRe — p(zz where, as defined for the piecewise linear case, J = — (3 21) ) 2 a (p2 2 — 2 + pi)(Ui pi)gh/(p — ) is 2 U the bulk Richardson number. We solve (3.21) numerically using the method described in the next section. 3.3.1 Numerical Method In this section we describe the method used to solve (3.21) for a viscous flow with diffusion and initial background profiles given by (3.20). For the inviscid problem with piecewise linear profiles, it was possible to have an unbounded domain (i.e., —oc < z < oc). Since we solve the present problem numerically, we must restrict ourselves to a finite domain. The boundary conditions used are = = 0 at z = +H/2, where H is the total height of the domain (figure 3.8). These are equivalent to setting u = w’ = p’ = 0 at z = +H/2. Following the method of Klaassen & Peltier [30] we use a truncated Fourier sine series to approximate q 5 and and seek solutions of the form I \ — — pZj — N Zn=l N . fl ThZZ1PTh S111 n7r(z+H/2) H n7r(z+ /2) H ‘ ( ) Chapter 3. Linear Stability Analysis 33 Substituting (3.22) into (3.21) we obtain the following. N —z {—(U—c)DnqnSn—UzzqnSn—JpnSn} N = n=1 N {(U—c)pnSn—pabnSn} = n=1 where Dn +a and Sn (nK/H) 2 N z DnpnSn, aPrRe n=1 I (3.23) I n7r(z + H/2) In order to solve for the coefficients H g and Pn we use Galerkin’s method. We multiply the above equations by sin mK(z + H/2)/H for rn = 1 .. . = sin N and integrate from z = —H/2 to z = H/2. We obtain 2N equations that can be written as N n1 D 2 ti(m,n)+t ( m,n) Dm J + Dmcbm+pm = aRe CEbm N n=1 for m = (m, fl)n + ti(m, fl)Pn} + aReprDmPm 3 {—t = CPm Ii (3.24) J 1,... , N, and where ti(m,n) (m, n) 2 t (m, n) 3 t 2 = H/2 mTr(z+H/2) I Usin sin nlr(z+H/2)d z, H Hi-H/2 H — 2 sin j-H/2 = = pH/2 2 j I H/2 J-H/2 Paz Sill m7r(z + H/2) sin n7r(z + H/2) dz, H H mr(z + H/2) n7r(z + H/2) sin —dz H H (3.25) (3.26) (3.27) We have reduced the system of ordinary differential equations (3.21) to a simple matrix eigen-problem (3.24) which is of the form AV = cV, (3.28) where A is a 2Nx2N matrix defined by (m,n))+ 2 (Dnti(m,n)+t An,m = An,m+N = An+N,m = (m,n), 3 —t An+N,m+N = ti(m,n) + Dm Dm e (3.29) (3.30) 8nm (3.31) jDm aRePr 6nm, (3.32) Chapter 3. Linear Stability Analysis mn 6 34 is the Kronecker delta function, and V= (3.33) N 0 , It should be noted that although the above method finds 2N values of the complex phase speed c and the corresponding eigenvalues, it is only the two most unstable modes that are of interest. These correspond to the values of c with the two largest values of ac. It was determined that N = 100 is usually sufficiently large to resolve eigenvalues and eigenvectors and is the value used, unless otherwise stated. 3.3.2 Results: Smooth Profiles We are interested in how displacing the density interface affects the stability of the flow. This corresponds to varying the parameter e. Lawrence et al. [38] have studied this ef fect for the inviscid case originally examined by Holmboe [23] where the only additional parameter is the bulk Richardson number J. Their results have been reproduced in sec tion 3.2.1 where the effects of boundaries were also studied. For smooth profiles we have the following additional parameters: the ratio R between the shear layer thickness and the density interface thickness, the Reynolds number Re, and the Prandtl number Pr. Yonemitsu [70] studied the viscous case without diffusion, but only looked for the most unstable mode (as opposed to two unstable modes, one moving to the right and another moving to the left, which corresponds with the results of Lawrence et al. [38]). Chapter 3. Linear Stability Analysis 35 Smyth et al. [53] have shown that an initial mean flow with arbitrary initial values of R and Pr will eventually reach a state where R Pr = = For this reason, we will set , linking the effects of Pr and R on the stability of the flow. Below, we discuss 2 R the effect of specific parameters on the linear stability of the flow. Varying R The most significant difference between the piecewise linear profiles (3.9) examined in section 3.2 and the smooth profiles (3.20) examined here is that we now allow the density interface to have a finite thickness, which we change by adjusting the value of the param eter R. Since the expected value of R can vary greatly depending on the stratifying agent and the temperature of the fluid (table 3.1), it is of interest to see how this parameter affects the stability of the flow. In our study of the effect of R on the stability of the flow with varying , we use the values R 3, 4, 6, and 8 (with Pr = 9, 16, 36, and 64, respectively). Although salt stratified flows have much larger values of the Prandtl number than examined here, we are unable to numerically resolve such a flow as there is insufficient diffusion to compute with a feasible value of N (see section 3.3.1). Never theless, with the range of Prandtl numbers chosen, we are able to observe trends in how the stability of the flow changes from thermally stratified to salt stratified flows. Results of Hazel [21] and those of section 3.2.1 indicate that, when H 10, horizontal boundaries do not significantly affect the stability of an inviscid flow. A similar result for the viscous case is discussed below. Based on the above, we set H 10 to insure that the horizontal boundaries do not unduly influence the stability of the flow. We also use Re=300 which Smyth [51] found to be sufficiently large to give good agreement with his inviscid solutions, facilitating comparison with results of section 3.2.1. For e = 0 we have computed the maximum growth rate in the Holmboe regime for the above values of R (figure 3.9). As R decreases, the maximum growth rate decreases. Chapter 3. Linear Stability Analysis Temperature p 3 kg/rn °C kg/rn. s /s 2 m Molecular diffusion of heat in pure water 999.843 1.79 x i0 0 13.4 x 10_a 5 999.967 1.52 x 10 13.6 x 10_a 10 999.702 1.31 x 10 13.8 x 10_a a 15 999.102 1.14 x 10 14.0 x 10 20 998.206 1.01 x 10 14.2 x 10_a 25 997.048 0.89 x 10 14.4 x 10_8 995.651 0.80 x i0 30 14.6 x 10_a Molecular diffusion of NaC1 (0.01 molar) in water 1000.322 1.79 x i0 0.1415 x 10_8 0 1000.436 1.52 x 10 0.1441 x 10_8 5 10 1000.162 1.31 x 10 0.1467 x i0 15 999.554 1.14 x i0 0.1493 x 10_8 20 998.653 1.01 x i0 0.1519 x iD_a a 25 997.490 0.89 x 10 0.1545 x 10 30 996.089 0.80 x i0 0.1571 x 10—8 Molecular diffusion of salts in sea water 1028.106 1.89 x i0 0.1360 x 10_8 0 1027.676 1.61 x i0 0.1386 x 10_a 5 10 1026.952 1.39 x 10 0.1409 x 10_a 15. 1025.973 1.22 x 10 0.1434 x 10_a 20 1024.763 1.09 x i0 0.1474 x 10 25 1023.343 0.96 x 10 0.1484 x 10_a 1021.729 0.87 x i0 0.1509 x 10_8 30 36 Pr 13.4 11.2 95 8.2 7.1 6.2 55 1264 1054 893 764 666 576 511 1352 1130 961 829 722 632 564 Table 3.1: Values of Prandtl number for 1) diffusion of heat in pure water, 2) diffusion of NaC1 (0.01 molar) in water, and 3) diffusion of salts in sea water. In all cases p was calculated using the international equation of state of sea water (see appendix 3 of Pond & Pickard [47]) with salinities of 0, 0.5844, and 35 pss, respectively. 1u was determined from table 16 in Sverdrup et al. [60] using salinities of 0, 0.5844, and 35 ppt, respectively. Values for the coefficient of thermal diffusivity (i for case 1) are from Fischer et al. [15]. For cases 2 and 3, c corresponds to the diffusivity of NaC1 in water. Values were taken from CRC Handbook of Chemistry and Physics [67] and are for molarity of 0.01 and 0.6 for cases 2 and 3, respectively. Where necessary, linear interpolation was used to determine coefficients that are not listed in referenced tables. Chapter 3. Linear Stability Analysis 37 It is well known that Holmboe instabilities do not exist when R = 1 and so there is some critical valne of R below which Holmboe instabilities do not occur. Using the method of least squares, we fit a curve of the form (acj)max = a/(b — R) + c, where a, b, and c are constants, to the data in figure 3.9. This gives us a critical value of R agrees with the results of Smyth [51] who found a critical value of R 2.39. This 2.4 by computing the stability diagrams for several values of R. We start by examining the results for e = 0 (figure 3.10). We observe the same behaviour found by Smyth [51] for the inviscid case with smooth profiles. Increasing R increases the region of instability for both right and left travelling waves: at a given value of J, the range of a for which the flow is unstable increases as R increases; at a given value of a, the range of J for which the flow is unstable also increases with increasing R. In addition, the growth rates of the fastest growing Holmboe instabilities (figure 3.11) increase with R. Growth rates of Holmboe instabilities when R is finite are smaller than the piecewise linear case which corresponds to R —* cc. When J is small, stratification is weak and we expect the flow to behave like the homogeneous case for which the value of R has no effect. Consequently the stability of the flow for small J does not change significantly as R varies. For larger values of J, the wave number at which the most unstable mode occurs increases with R. Also, transition from Kelvin-Helmholtz instabilities to Holmboe instabilities occurs at smaller values of J as R increases. When 6 > 0, the region of instability for both left and right modes increase as R increases (figures 3.10). Also, the growth rate of the fastest growing mode increases with increasing R (figure 3.11). These results are the same as the e = 0 case. The wave number at which the most unstable mode occurs, however, changes with both If we fix 6, 6 and R. then the wave number at which the fastest growing mode occurs increases with increasing R; conversely, if we fix R and vary e the behaviour depends on the value of R, as described below. Chapter 3. Linear Stability Analysis 38 For a piecewise linear inviscid flow, Lawrence et al. [38] found that as e increases, the most unstable right travelling wave occurs at a smaller wave number, whereas the most unstable left travelling wave occurs at a larger wave number. On the stability diagrams, this corresponds to the tilting of the limb corresponding to the right moving wave (the right limb) towards the J-axis, and the limb corresponding to the left moving wave (the left limb) towards the a-axis (figure 3.2). This behaviour remains unchanged for large finite values of R (figure 3.10, R = 6 and 8). As R decreases below some critical value of approximately 4, however, this behaviour starts to change (figure 3.10, R R = = 4). For 4, there is very little tilting of the right limb and the wave number at which the most unstable right moving wave occurs remains almost constant as increases. When R = 3, however, the wave number of the most unstable right mode increases with increasing e (figure 3.10, R = 3). This change in behaviour with changing density interface thickness was documented by Caulfield [7] for three layer inviscid flows. It is seen that for all values of R, the left limb tilts towards the a-axis. Also, the left limb lifts off the a-axis as increases (figure 3.10). This effect becomes more noticeable as R decreases. This did not occur in the findings of Lawrence et al. [38] since their case corresponds to R —* c. The wave number of the most unstable left mode increases as e increases, except when R and 0 < < 0.25. In this case, a* decreases as E = 3 increases. This is a result of the large separation of the left mode from the a-axis. The separation of the left limb from the a-axis affects the range of J for which both left and right moving unstable modes exist. As E decreases, the range of J for which both modes exist decreases. Varying H Above, it was assumed that H = 10 is sufficiently large for the boundaries to have little affect on the stability of the flow. In section 3.2.1 we found this to be true for an inviscid flow with piecewise linear profiles. We now examine the validity of this assumption for Chapter 3. Linear Stability Analysis 39 smooth profiles. Having studied the effect of R in the previous section, we now set R = 3 since it is sufficiently large to allow Holmboe instabilities, yet small enough to converge without requiring an excessive number of Fourier modes. Also Pr = 2 R = 9 is a physically realistic value of the Prandtl number, corresponding approximately to the diffusion of heat in water at 10°C (see table 3.1). As above we set Re When e = = 300. 0 (figure 3.12), the effect of H on the stability of the flow differs markedly from the piecewise linear case. For H 10 the boundaries do not significantly affect the stability of the flow except at small values of both a and J. This is not surprising in light of the results from the piecewise linear case. The fastest growing mode is not affected by the presence of horizontal boundaries until H < 10, at which point its growth rate decreases as H decreases (figure 3.13). An interesting difference between the effect of H on the piecewise linear profiles and smooth profiles is how the transition between Kelvin Helmholtz instabilities and Holmboe instabilities changes as H varies. For the piecewise linear profiles, it was observed that as H decreases the transition occurs at smaller values of J until only Holmboe modes are present at H 4. This was explained by the Kelvin-Helmholtz instabilities being stabilized first as they occur at smaller wave numbers than the Holmboe instabilities. In contrast to the piecewise linear case where Holmboe instabilities exist for all J > 0, for smooth profiles there is a range of J between 0 and some critical value J (which depends on the value of R) for which Holmboe instabilities do not exist. In addition, for smooth profiles, instabilities near the transition between Kelvin-Helmholtz and Holmboe instabilities occur at smaller wave numbers and have the smallest growth rates. These are the first modes to be stabilized by the boundaries. This manifests itself by the separation of the two regions of instability with a stable region in between. As H decreases, the region separating the two modes increases while the growth rates of the two modes decrease. When e = 0.25 with H = 20 (figure 3.13), the growth rate of fastest growing left mode Chapter 3. Linear Stability Analysis 40 is much smaller than that of the fastest growing right mode. As H decreases, the growth rates of the two modes decrease. Overall, the growth rate of the right mode decreases more rapidly than that of the left mode as H decreases. Nevertheless, the right mode always remains the faster growing of the two modes. This is in contrast to the piecewise linear case where for small H the left mode dominates. For the piecewise linear case with e = 0.25 and large H (figure 3.3), the difference between growth rates of left and right modes is not great, and thus the rapid drop in growth rate of the right mode as H decreases allows it to fall below the growth rate of the left mode. This is not possible for the case examined here since the initial difference between the growth rates of right and left modes is so large. Varying Re In this section we examine the effect of varying the Reynolds number, Re, on the linear stability of a mean flow described by (3.20). We study the effect of Re by computing the stability diagram for Re = 300 and 100 with R 3. These are shown in figure 3.14. As Re decreases, viscosity becomes more important and we expect an overall stabi lization of the modes with decreasing Re. Since viscosity has the greatest effect on small wave lengths, it is primarily the instabilities at large wave numbers which are stabilized as Re decreases. This has the effect of decreasing the wave number at which the most unstable mode occurs. When e = 0, growth rates of the most unstable modes decrease with decreasing Re and are affected primarily in the transition between Kelvin-Helmholtz and Holmboe instabilities. When E > 0, the fastest growing right mode decreases as Re decreases. There is a slight increase in the growth rate of the most amplified left mode as Re decreases. Although this is surprising, the difference in growth rates is not significant. It is difficult to draw any conclusions about the left mode without further investigation. Chapter 3. Linear Stability Analysis 41 Varying e: summary We now discuss the overall effect of varying . It was seen that precise behaviour of the flow stability with increasing r depends on the thickness of the density interface (which is specified by the parameter R). Nevertheless, there are trends that were observed in all the cases examined. As with the piecewise linear case, there are two types of instabilities when 0. For small J, the flow is unstable to Kelvin-Helmholtz instabilities, whereas for larger J it is unstable to Holmboe instabilities. When e > 0, there are no pure Kelvin Helmholtz or Holmboe instabilities. The Kelvin-Helmholtz instability is replaced by an instability with positive phase speed. The growth rate of this instability increases with . For larger J, there are still two unstable modes, one moving to the right and another to the left. As e increases, however, the growth rate of the right moving wave increases whereas that of the left moving wave decreases. This is consistent with the results of Lawrence et al. [38]. On the stability diagram this corresponds to a splitting into two limbs. The limb corresponding to right travelling waves gets larger. Thus as e increases, the maximum value of J for which the flow is unstable increases. The limb corresponding to left travelling waves shrinks. Hence, the interval of J over which there are two unstable modes with opposite phase speeds decreases. It was also seen that as E increases, the phase speed of the right mode decreases whereas that of the left mode increases. The most significant difference between the smooth profile case and the piecewise linear case examined in section 3.2 is that when > 0, the difference in growth rates between right and left modes is much larger for the smooth profile case studied here. We expect this difference to be largely due to the different background profiles used in the two studies. This is supported by the similarity between our results for a smooth background flow with the effects of viscosity and diffusion included and with = 0 and those of Smyth [51] who studied the stability of an inviscid flow with the same smooth background profiles. Chapter 3. Linear Stability Analysis 3.3.3 42 Three-Dimensional Instabilities: Smooth Profiles We wish to examine the possibility of the instabilities initially being three-dimensional when viscosity and diffusion are included. In this case, the relationship between twodimensional and three-dimensional instabilities for an inviscid flow (3.16) becomes cc(U, ,a, Re, J, Pr) cos = ac(U, 0, c, Recos 9, J/ cos 2 9, Pr). This is Koppel’s result [35] which reduces to Squire’s theorem [57] when J = (3.34) 0. In order for the initial instabilities to be two-dimensional, we require cos 9. g*cj(U 0, c, Recos 9, J/ cos 2 , Pr) < c*cj(U, 0, a, Re, J, Pr), (3.35) where c is the wave number at which the growth rate is maximized for given values of J, Re, and Pr. We can conclude from the above that the initial instabilities may be threedimensional if the growth rate of the fastest growing mode increases with decreasing Re and/or increasing J. Above, we examined the effect on the stability of the flow of decreasing the Reynolds number. It was observed that this tends to decrease the growth rates of the right mode when e > 0, and both modes when 0. Therefore the only way that the three- dimensional instabilities can dominate is if the growth rate of the fastest growing mode increases with increasing J. When e = 0, there is a region where the most amplified growth rate increases with increasing J. This indicates that three-dimensional instabilities may grow faster than two-dimensional instabilities. In order to verify this, one would have to solve the three-dimensional linear stability problem to see if the initial perturbations are indeed three-dimensional. Smyth & Peltier [55] have done this for a range of Re with R = 3 and 0. They were able to show that, for smaller values of Re, the flow is initially unstable to three-dimensional perturbations. This includes the value of Re = 300 examined here. When > 0, the growth rate of the fastest growing right mode Chapter 3. Linear Stability Analysis 43 decreases with increasing J (figure 3.11) implying that the most unstable instabilities are two-dimensional. The behaviour of the left mode when e > 0 is different than that of the right mode. First, we observed a slight increase in the growth rate of the most amplified left mode when the Reynolds number decreased from 300 to 100. Also, there is a region where the growth rate of the left mode increases with increasing J (figure 3.11). Thus it may be possible that the left mode is most unstable to three-dimensional perturbations. Further investigation is required to determine whether or not this is true. Although the left mode may be unstable to three-dimensional instabilities, its small growth rate compared to that of the right mode when e > 0 make it unlikely that we will initially see threedimensional instabilities. Therefore, for the remainder of this study, we will concentrate on the development of two-dimensional instabilities. It may be possible, however, that the eventual growth of the left mode is one mechanism through which the perturbations become three-dimensional at later times. 3.3.4 Eigenfunctions So far, we have only examined the linear growth rates to determine the overall stability of the flow. In order to better understand the mechanism of the instability in the linear regime, we now consider the eigenfunctions found from linear stability analysis. Since, for this study, we are interested in flows that have two unstable modes, a right moving wave and a left moving wave, we must examine not only the eigenfunctions of each mode separately, but also the superposition of the two. We restrict our attention to the study of the modes found by linear analysis of viscous flows with hyperbolic tangent profiles as these are the starting point of the simulations described in chapter 4. Since linear stability analysis only determines the eigenfunctions up to a constant, we normalize the eigeufunctions so that they have the same amount of perturbation kinetic energy in each Chapter 3. Linear Stability Analysis 44 mode. In appendix C we derive the equation governing the evolution of the horizontally averaged perturbation kinetic energy, given by (C.14). In the linear regime, we can neglect higher order terms in the perturbations and obtain the following equation for the development of the perturbation kinetic energy. 1 (u12 + w12) = — (w’p’) — + _ (u’An’ + w’’). (3.36) The left hand side of (3.36) is the rate of change of perturbation kinetic energy which is equal to the sum of the following terms (see Bradshaw [3] for homogeneous viscous flows and Smyth & Peltier [54] for stratified inviscid flows). 1. 2. —ai7 represents the extraction of energy from the mean flow by the perturbation, (W’p’)z represents the change in perturbation kinetic energy associated with the transfer of perturbation kinetic energy by the pressure fluctuations, 3. —Jw’p’ represents the conversion of perturbation kinetic energy into potential en ergy, 4. u’Au’ + w’Aw’ can be divided into a viscous diffusion term and a viscous dissipation term. In order to gain a better understanding of the development of the instabilities, we examine terms on the right hand side of (3.36) by calculating the correlations u’Au’ + w’z.w’. In addition, we examine the eigenfunctions , ii’, I, and ‘, and determined from linear stability analysis. Although we do not graph the complex amplitude of the stream function perturbation related by = i’th/o. , it can be inferred from the graph of i’ since the two are Chapter 3. Linear Stability Analysis 45 Eigenfunctions of Individual Modes We start by examining the eigenfunctions for the Holmboe instability corresponding to R = 3, Pr = 9, Re = 300, J 0.3, and E = 0. This corresponds to the case studied by Smyth & Peltier [54] for an inviscid fluid and gives us a basis for comparison for the 0 case. The eigenfunctions and correlation terms in (3.36) for the right and left modes of the Holmboe instabilities are shown in figure 3.15. Since e 0 the growth rates of the two modes are equal and the phase speeds are equal and opposite. The horizontal velocity perturbation amplitude, ü(z), is concentrated in a narrow region centered near the critical level (the height z, such that U(z) = Cr) amplitude, ,(z), is concentrated near z of the mode. Also, the density perturbation 0, the center of the density interface, and is due to the strong gradient of Pa at this point. The vorticity perturbation amplitude & has two peaks, one just above (below) the density interface and another just above the critical level (below) for the right (left) mode. For both left and right modes the maximum value of the perturbation kinetic energy, 2 + w’ u’ , occurs near their critical levels. Since 7 2 < 0, energy is being extracted from the mean flow. This indicates that the modes are unstable. It can be shown (see Klaassen & Peltier [29]) that iii can be expressed in terms of the slope of the phase of z’, i.e., = where ?1’ = r(z) exp(i(z)). A positive slope of (3.37) indicates that the mode is growing, if it is negative then it is decaying. For the eigenfunctions shown in figure 3.15 the slope of the phase of ti’ is positive and thus the perturbations are unstable. This is in agreement with the positive growth rate &c = 0.030561 and the correlation The minimum of occurs near the critical level and so the perturbations extract energy from the mean flow most efficiently at this level. The correlation represents a vertical flux in the perturbation kinetic energy. If Chapter 3. Linear Stability Analysis 46 iF > 0 then the kinetic energy is directed upward. It is directed downward when < 0. For the right mode, 7j < 0 and hence perturbation kinetic energy is being moved downward from the critical level towards the centre of the density interface. For the left mode, the perturbation kinetic energy is being moved upward from the critical level to the density interface. The quantity corresponds to a vertical displacement of mass. When i77> 0 mass is moved against the gravitational restoring force, causing an increase ill the potential energy of the flow. Potential energy is converted into perturbation kinetic energy when mass is displaced in the direction of the gravitational restoring force (i < 0). For both the right and left modes > 0 indicating that perturbation kinetic energy is being converted to potential energy. From the above we can describe the growth of the perturbations as follows. The perturbations extract kinetic energy from the mean flow at the critical level. The perturbation kinetic energy is then transported towards the center of the density interface where a portion of it is converted to potential energy. When e = 0.25 (not shown), the general mechanism for growth and the shapes of the eigenvalues are similar to the e = 0 case. The horizontal velocity perturbation has a jet concentrated near the critical level. Also, the density perturbation amplitude, ,3, is concentrated near the centre of the density interface which now occurs at z —e. There are, however, some differences in the relative amplitudes associated with the two modes which did not occur for the = 0 case. The magnitude for the density perturbation of the left mode is greater than that of the right mode. The vorticity & has two peaks, one near the density interface and one near the critical level, but instead of being equal in amplitude, as in the = 0 case, the peak near the density interface is greater. Also, the left mode is starting to develop vorticity near the critical level of the right mode. As with the 0 case, the perturbation kinetic energy, u’ 2 + w’ , is largest near the 2 critical levels. The magnitude of iF is larger for the right mode, indicating that it is Chapter 3. Linear Stability Analysis 47 extracting energy from the mean flow more efficiently and hence will grow more quickly. Finally, and are larger for the right mode than the left. This indicates that the right mode is responsible for a larger portion of the vertical displacement of perturbation kinetic energy. Also, most of the conversion from perturbation kinetic energy to potential energy occurs through the right mode. When = 0.5 (see figure 3.16), the differences between the right and left modes are even more pronounced. The overall shapes of the eigenfunctions for the right mode are not significantly different from the e 0 case. The eigenfunctions of the left mode, however, have been altered dramatically by the displacement of the density interface. This is similar to results of Yonemitsu [70] who found that for non-symmetrically placed horizontal boundaries only one mode was affected. it is clear, from 7j and Z7 that the right mode is the main mechanism through which energy is exchanged between the various reservoirs (mean kinetic, perturbation kinetic, and potential energies). For the left mode, the most efficient extraction of energy from the mean flow no longer occurs at the critical level. Instead it has been shifted upward towards the density interface. Superposition of Right and Left Modes To better understand how right and left modes interact in the linear regime, we su perimpose the two modes and examine the evolution of the perturbation kinetic energy throughout one cycle of its period. In appendix D it is shown that the period of the horizontally averaged kinetic energy is 2K/(c) — 41)). We examine the evolution of the terms in (3.36) for one period of the kinetic energy of the superimposed modes in increments of Ax = r/4 where x = a(4r) — cO)t. For convenience, we suppress the linear growth components of the perturbations. The evolution of the terms in (3.36) throughout a cycle is shown in figure 3.17 for the = 0 case. In the linear regime, growth and decay of the perturbation kinetic energy is Chapter 3. Linear Stability Analysis controlled by 48 (which determines the transfer of kinetic energy from the mean flow to the perturbation). The perturbation extracts energy most efficiently at corresponds to a local maximum in the perturbation kinetic energy. x = 0. This ir corresponds to the time when the most amount of perturbation kinetic energy is being returned to the mean flow, indicating a local minimum in the perturbation kinetic energy. In between these local extrema there is a period of decay for 0 < for r < x < 2ir < r and a period of growth after which the cycle repeats itself. Although the perturbations spend equal amounts of time in the growth and decay portions of the cycle, the amplitudes of are largest when the perturbations are extracting kinetic energy from the mean flow. Thus if growth of the modes were included, we would observe a net gain in perturbation kinetic energy over one cycle. This oscillation in the growth pattern of Holmboe waves was first predicted by Holmboe [23] and later verified by Smyth et al. [53], and Smyth & Peltier [54]. In chapter 4, we also observe this result numerically, and examine how the exchange of energy changes in the nonlinear regime. Examination of the correlations w’p’, and w’p’, enable us to determine how en ergy is transferred between the three energy reservoirs: mean kinetic energy, perturbation kinetic energy, and potential energy. When energy is transferred from the mean flow to the perturbations (i.e., i7 < 0), the critical levels act as sources for the perturbation kinetic energy. Perturbation kinetic energy is then transported vertically away from the critical level towards the density interface. The density interface then acts as a sink, converting perturbation kinetic energy into potential energy. This is as described above from the results of the eigenfunctions for the individual modes. The opposite behaviour occurs when the perturbations return kinetic energy to the mean flow. The critical levels act as sinks and the density interface as a source. Potential energy at the density inter face is converted back into perturbation kinetic energy which is transported towards the critical levels where the perturbations return energy to the mean flow. There is a phase Chapter 3. Linear Stability Analysis difference, however, between i7 and 49 At = 37r/4 perturbation kinetic energy is still being transferred to potential energy even though it is also returning energy to the mean flow. Similarly at x = 37r/2, potential energy is being converted to perturbation kinetic energy as is mean kinetic energy. Figure 3.18 shows the evolution of the perturbation kinetic energy and the correlations for = 0.25. There is essentially the same oscillation between growth and decay as described for minimum when 0. The perturbation kinetic energy is maximum when x = 0 and ir with corresponding periods of decay and growth in between these local extrema. As with the e energy is in phase with flow x = 0 case, the growth and decay of the perturbation kinetic which determines the exchange of energy between the mean and the perturbations. Eigenfunctions shown in figure 3.16 indicate that we can, for the most part, use the line z = —e to separate contributions to left and right modes, particularly if we focus on the behaviour near their respective critical levels. In figure 3.18 we observe that the right mode extracts energy from the mean flow more efficiently than the left mode. During the time when the perturbations are returning energy to the mean flow, the right mode returns less than the left mode. These two features indicate that the right mode grows more rapidly than the left mode. This is indeed the case, as determined by the values of oc predicted from linear analysis. Upon examination of it is evident that more perturbation kinetic energy is being moved around in the upper layer than the lower layer. Thus the right mode is responsible for much of the vertical transport of perturbation kinetic energy. We believe that the right mode is probably responsible for most of the exchange between perturbation kinetic energy and potential energy, although this is difficult to determine from the evolution of 2 + w’ The evolution of u’ , and the corretations in equation 3.36 for e 2 = 0.5 are shown in figure 3.19. As with the two previous cases, e 0 and 0.25, the perturbation kinetic energy oscillates between a maximum value at —ir/4 and a minimum at x x = 3ir/4. Chapter 3. Linear Stability Analysis There is a phase difference, however, between n’w’ and 50 2 it’ , whereas for e 2 + w’ = 0 and 0.25 they was not. The most efficient transfer of energy from the mean flow to the pertnrbations occur at when x x= 0. The perturbations return the most energy to the mean flow Also, the perturbations spend more time throughout the cycle extracting energy from the mean flow than in the previous cases. If we continue to use z = —e to separate the behaviour associated with the right mode from that of the left mode, we observe that right and left modes no longer exhibit the same behaviour throughout the cycle. There is a slight phase shift between the two modes. The right mode start to release energy to the mean flow and resumes its extraction of energy from the mean flow at different times than the left mode. As with the e 0.25 case, the right mode extracts energy more efficiently from the mean flow and returns less than the left mode. Also, the right mode is the primary means through which perturbation kinetic energy is transported vertically towards the center of the shear layer. Chapter 3. Linear Stability Analysis Figure 3.1: Background flow for piecewise linear velocity and density profiles 51 Chapter 3. Linear Stability Analysis 52 e=o £O.25 H=ou H=ou 0.50 u.Du 0.40 0.40 H=ou 0.50 U YJ4’J 0.40 0.30 uhf 0.30 0.20 0.20 0.20 0.10 0.10 0.10 Ill wiiii iii — dill/f !A — — - 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 C 1 0.5 1.0 a 1.5 2 H=1O 0.50 0.40 0.30 0.20 0.10 0.0 0.5 1.0 a 1.5 2.0 0 H=5 H=5 0.50 0.50 0.50 0.40 0.40 0.40 0.30 0.30 0.20 0.20 0.10 vi, 00 0.5 1.0 1.5 0.10 ;E 0.00! S ) 00 1.0 0.5 1.5 .0 H =4 Il/Il I 0.40 — .Q,I/A / 0.40 / // ‘‘2,; / ‘2/I 0.30 0.30 / I IIaj 0.20 — 0.20 — — — 0.10 0.10 / — / — / / 1. — .c 5 ____. ‘%- 0.01 0.0 0.5 1.0 a 1.5 21 1 :r 0.5 0 1.0 1.5 l//// I 0.10 Q.OO 1.5 Figure 2.0 n on 0.0 ,7 — ‘ 0.20 1.0 0.5 1.0 1.5 .0 H=3 -1’ 0.40 •‘7 0.5 - U H=3 0.30 0.00 0.0 Out. —‘0.0 — 0.40 0.30 0.20 — — — — (/77/ 0.10 i’l’ __u. — — U= 0.5 3.2: Stability — 1.0 1.5 2.0 0.00 0.0 0.5 I 1.0 — —: - 1.5 2.0 diagrams for inviscid case stability diagrams for piecewise linear profiles with = 0, 0.25, and 0.5 for H = oo, 10, 5, 4, and 3. Contours are of constant growth rate For > 0, solid contours correspond to right travelling waves and dashed contours to left travelling waves. When 0, growth rates of left and right modes coincide in the Holmboe regime whereas in the Kelvin-Helmholtz regime only dominant modes are shown. Stability boundaries (ac = 0) are the outer most contour in each limb. Linear Chapter 3. Linear Stability Analysis 53 8=0.25 8=0 H=co H=oo 0.20 0.20 0.15 0.15 0.10 0.10 0.05 0.05 0.10 0.20 0.30 0.40 0.15 a / / 0.10 0.05 0.0c 0.00 0J 0 H=aa 0.20 S 0.00 8=0.5 0.10 0.20 H=1O 0.30 0.40 0J 0 0.00 0.10 H=1O 0.20 0.30 0.40 0.50 0.40 0.50 0.40 0. 0.40 0. H=1O u.zD 0.20 0.20 0.15 - 0.20 0.15 0.15 a a 0.10 0.10 0.05 0.05 0.00 0.10 0.20 0.30 0.40 0, 0 a / / 0.10 0.05 0.00 0.10 0.20 H=5 0.40 0.30 0.50 0.00 R=5 0.20 0.25 0.20 0.20 0.20 0.15 0.15 0.15 0.10 0.10 0.10 0.05 0.05 0.05 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 H=4 0.20 0.30 H=5 0.20 0.00 0.10 0.30 0.40 0.50 0.00 0.10 H=4 0.20 0.30 H=4 0.20 0.20 0.20 0.20 0.15 0.15 0.10 0.10 0.10 0.05 0.05 0.15 a a 0.00 0.10 0.20 0.30 0.40 0.50 /*— 0.00 0.10 0.20 H=3 0.30 0.40 0.05 0.1 0 0.50 0.10 H=3 0.25 0.25 0.20 0.20 0.20 0.15 0.15 a 0.10 0.10 — — 0.20 0.30 H=3 0.15 — a — — — — — 0.10 7 — / 0.05 0.05 0.00010 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.00 0.10 0.20 0.30 0.40 0.50 Figure 3.3: Growth rates of most unstable modes for inviscid case Growth rate of most unstable mode at a given J for stability diagrams shown in figure 3.2. Maximum growth rate in the flolmboe regime (e = 0) is indicated by the horizontal dotted line. Chapter 3. Linear Stability Analysis 54 8=0 8=0.25 H=co 0.00 0.00 0.10 0.10 0.20 0.30 8=0.5 H=co 0.40 0.50 0.00 0.10 0.20 0.30 H=co 0.40 0.50 ::__ 0.00 H=iO H=1O H=1O 0.20 0.30 0.40 0.50 0.10 0.00 H=5 . —0.5 0.20 ., 0.30 0.40 0.50 0.00 —is. 0. 0 .. 0.30 0.40 0.50 0.30 H=5 H=5 ç — 0.40 0.50 0.40 0.50 0.40 0.50 0.40 0.50 0.40 0.50 — .. 0.10 0.20 0.30 J J H=4 H=4 1.0 0.40 0.50 0.00 0.10 0.20 0.30 S H=4 1.0 —0.5 0.00 0.20 J —0.5 0.20 0.10 J — 0.00 0.30 J 1.0 0.10 0.20 J J —Is 0.00 0.10 J —0.5 0.10 0.10 0.20 0.30 0.40 0.50 0.00 —0.5 0.10 0.20 0.30 0.40 0.50 —1.0 0.00 0.10 0.20 0.30 J S S H=3 H=3 H=3 0.20 0.30 J 0.40 0.50 0.00 0.10 0.20 0.30 S 0.40 0.50 0.00 0.10 0.20 0.30 J Figure 3.4: Phase speeds of most unstable modes for inviscid case Phase speed of most unstable mode at a given J for stability diagrams shown in figure 3.2. Chapter 3. Linear Stability Analysis 0.0 0.2 55 0.4 0.6 0.8 1.0 0.6 0.8 1.0 a a 0.10 - 0.08 0.0 Dominant Kelvin—Helmholtz Mode 0.2 0.4 a Figure 3.5: Linear stability diagrams for small J when E = 0 Growth rates (solid contours) and phase speeds (dotted lines) for the case studied by Holmboe [23]. When J < 0.07 there are two regions of instability: Kelviri-Helmholtz instabilities (cr = 0) and Holmboe instabilities (the superposition of right and left moving waves). Chapter 3. Linear Stability Analysis 56 0.20 0.15- g 0.10 0.05- 0.00 L. 0.00 0.10 0.20 0.30 0.40 0.50 0.30 0.40 0.50 Li 4 -3 -7 0 0 0.00 0.10 0.20 J Figure 3.6: ac/Vi for Kelviu-Helmholtz and Holmboe waves Growth rate, ac,j and ac/i/J along a-J curve for which growth rate is maximized for H = oc, e = 0. Results for both Kelvin-Helmholtz (dotted line) and Holmboe (solid line) instabilities are shown. Chapter 3. Linear Stability Analysis 57 8=0 8=0.25 6=0.5 0.10 0.10 0.10 0.05 0.05 0.05 0.00 0.00 o.oc 0.10 0.20 0.30 0.40 0.E 0.00 0.10 0.20 J 0.30 0.40 0.50 0.00 2.G 2X 1.5 1.5 1.5 1.0 l.o 1.0 0.10 0.20 0.20 0.30 0.40 0.50 0.00 0.10 0.20 J 0.30 0.40 0.50 0.00 J Figure 3.7: ac//J for 0.30 0.40 0.50 0.30 0.40 0.50 J 2.0 0.00 0.10 J 0.10 0.20 J = 0, 0.25, and 0.5 Growth rate, cc and c/\/J along a-J curve for which growth rate is maximized for H = oo, with = 0, 0.25, and 0.5. Both right (solid line) and left (dashed line) modes are shown. Chapter 3. Linear Stability Analysis z Figure 3.8: Background flow for hyperbolic tangent velocity and density profiles 58 Chapter 3. Linear Stability Analysis 59 0.10 0.08 0.06 0 E -Th C) 0.04 0.02 0.00 0 2 4 6 8 10 R Figure 3.9: Maximum growth rate in Holmboe regime with varying R Maximum computed growth rate in flolmboe regime (+) for R 3, 4, 6, and 8 (Pr=9, 16, 36, and 64, respectively) with e = 0, Re=300, and H = 10. The solid line is the approximation (cj)max = 0.2234/(0.5544 R) + 0.1218 determined using the method of least squares. — Chapter 3. Linear Stability Analysis 60 8O.25 0.10 0.0 0.5 a 1.0 a 1.5 2.0 0.00 0.0 0.5 1.0 a 1.5 2.0 0.5 1.0 a 1.5 2.0 1.5 2.0 -3 0.10 0.0 0.5 1.0 a 1.5 2.0 1.0 a 2.0 0.00 0.0 R=6 0.50 0.40 0.30 0.20 0.10 fi 0.0 0.5 1.0 a 1.5 2.0 nfl 0.0 0.5 1.0 a 1.5 2.0 a R=8 R8 0.50 0.50 0.40 0.40 0.30 0.30 0.20 0.20 0.10 0.10 -3 n fin 0.0 0.5 1.0 a 1.5 2.0 0.0 fi 0.5 1.0 a 1.5 2.0 fin 0.0 0.5 1.0 a Figure 3.10: Linear stability diagrams for smooth profiles: varying R Linear stability diagrams for e = 0, 0.25, and 0.5 with Re=300, H = 10, and R 3, 4, 6, and 8 (Pr = 9, 16, 36, and 64, respectively). Contours are of constant growth rate act. For e > 0, the solid contours correspond to right travelling waves and the dashed contours to left travelling waves. When 0, growth rates of left and right modes coincide in the flolmboe regime whereas in the Kelvin-Helmholtz regime only dominant modes are shown. Stability boundaries (oc, 0) are the outer most contour in each limb. Chapter 3. Linear Stability Analysis 61 =O.25 R=3 R3 0.15 0.10 U 0.05 0.05 0.001 .-.,— 0.001 —— 0.00 0.10 0.20 0.30 0.40 0.50 J “.— :i.. — 0.00 0.10 0.20 0.30 0.40 0.50 J R=4 0.15 0.10 () : 0.10 0.05 0.05 0.05 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J 0.00 L 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J R=6 R=6 0.15 R= 6 0.15 0.10 0 0.15 0.10 0 0.10 0.05 0.05 0.05 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J / n fin 0.00 0.10 0.20 0.30 0.40 0.50 J 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J R=8 R=8 0.15 0 0.15 0.10 0 0.05 R=8 0.15 0.15 0.10 :i 0.10 0.05 0.05 / 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J fi fin 0.00 0.10 0.20 0.30 0.40 0.50 J . .. .N 0.00 0.00 0.10 0.20 0.30 0.40 0.50 J Figure 3.11: Maximum growth rates for smooth profiles: varying R Maximum growth rate at a given J for stability diagrams shown in figure 3.10. Maximum growth rate in the Holmboe regime (E = 0) is indicated by the horizontal dotted line. Chapter 3. Linear Stability Analysis 62 é=O.25 H =20 0.10 0.Oc 0.0 0.5 1.0 1.5 2.0 a H=10 0.0 0.5 1.0 1.5 2.0 2 a ,, 0.5 1.0 1.5 2 0.5 1.0 1.5 21 1.5 2.0 H=4 :: 0.10 0.00 0.0 0.5 1.0 1.5 21 0.00 0.0 H=3 0.51, E \ 0.20 0.0 ‘.,—‘——-.——_ 0.10 0.0 0.5 1.0 1.5 2.0 0.00 0.0 .2 — 0.5 1.0 Figure 3.12: Linear stability diagrams for smooth profiles: varying H Linear stability diagrams for R = 3, Pr=9, Re=300, = 0, and 0.25 with H = 20, 10, 5, 4, and 3, (N 150, 100, 100, 100, and 100, respectively). For > 0, the solid contours correspond to right travelling waves and the dashed contours to left travelling waves. When = 0, growth rates of left and right modes coincide in the Holmboe regime whereas in the Kelvin-Helmholtz regime only dominant modes are shown. Stability boundaries (ac = 0) are the outer most contour in each limb. Chapter 3. Linear Stability Analysis 63 S=O.25 H=20 H=20 o.01 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0. 0.30 0.40 0. 0.30 0.40 0. 0.30 0.40 0.0 -- H=1O 0.15 0.10 0.05 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 H5 H=5 0.15 0.15 0.10 a 0.05 0.05 10 0.10 0.20 0.30 0.40 0. 0.00 0.00 ..-. 0.10 0.20 H=4 H =4 0.15 0.15 0.10 0.10 0.05 0.05 a nn 0. 10 0.10 0.20 0.30 0.40 0. 0.00 0.00 0.10 0.20 H=3 H=3 0.15 0.15 0.10 0.10 0.05 0.05 a n on 0.00 0.l0 0.20 0.30 0.40 0.50 non 0.00 — 0.10 0.20 0.30 0.40 —. —. 0.50 Figure 3.13: Maximum growth rates for smooth profiles: varying H Maximum growth rate at a given J for stability diagrams shown in figure 3.12. Maximum growth rate in the flolmboe regime (e 0) is indicated by the horizontal dotted line. Chapter 3. Linear Stability Analysis 64 8=0 80.25 8=0.5 0.50 0.40 0.30 -3 0.20 0.10 0.0 0.5 1.0 a 1.5 2.0 0.0 0.5 1.0 a 1.5 2.0 0.00 0.0 0.5 Ro=300 0.50 0.40 0.30 0.20 0.10 (1 a 0.5 1.0 a 1.5 2.0 1.5 2.0 Re=300 “I Or 0.0 1.0 a 0 1.5 2.0 0.0 0.5 1.0 a R= 100 0.15 0.10 a 0.05 0.00 0.10 0.20 0.30 0.40 0.50 0. )O 0.10 0.20 0.30 0.40 0. 50 Re300 0.15 0.10 0.05 0.00 0.10 0.20 0.30 0.40 Figure 3.14: Smooth profiles: varying Re Linear stability diagrams for R = 3, Pr=9, H 10, Re=100, and 300 with E = 0, 0.25, and 0.5. Here, N = 50 and 100, for Re = 100 and 300, respectively. Also shown is the maximum growth rate at a given J. When > 0, solid lines correspond to right travelling waves and dashed lines to left travelling waves. When E = 0, growth rates of left and right modes coincide in the Holmboe regime whereas in the Kelvin-Helmholtz regime only dominant modes are shown. Stability boundaries (ac = 0) are the outer most contour in each limb. Maximum growth rate in the Holmboe regime (E = 0) is indicated by the horizontal dotted line. Chapter 3. Linear Stability Analysis 65 a) p phase/rr b) w phase/Tr c) u phase/7T 1.5 0.0 amplitude 0.00 0.15 amplitude 0.0 0.3 amplitude — d) c phase/IT — e) p phase/Tr N f) 2 + (u ) /2 w g) uw h) wp 0.00 0.65 amplitude 0.00 0.26 amplitude j) uu+ww 0.015—3x10 5 3x10 i) wp 5 ::.* NQ -5 0.00 0.02 —0.011 0.011 Right mode, ac a)p phase/IT = —0.003 0.030561, cr 0.003 —0.015 = 0.598220, and Zc = 0.69037. b)w phase/IT c)u phase/IT d)i phase/IT e)p phase/IT 0.00 0.15 amp1 it u d e 0.0 0.3 amplitude 0.00 0.65 amplitude 0.00 0.26 amplitude 5 N —0 —5 0.0 1.5 amplitude f) 2 + (u ) /2 w g) uw h) wp j) i) wp uu+wtw 5 N —0 -5 0.00 .. 0.02 —0.011 Left mode, ac 0.011 = —0.003 0.030561, Cr 0.003 —0.015 5 3x10 0.015—3x10 5 —0.598220, and z Figure 3.15: Eigenfunctions and correlations when —0.69037. E = 0 Eigenfunctions and correlations when e = 0, R = 3, Re = 300, Pr = 9, H = 10, J = 0.3, and a 0.54, a), bth, c) ü, d)’, e), f) (u’ )/2, g) iF7, h) 2 i) 2 + w’ i) u’L\u’ + w’/w’. z = —e is indicated by the dashed horizontal line and the critical levels are the horizontal dotted line. For a) through e), both amplitude (thick line) and phase (thin line) are shown. In a) the phase is only given when the amplitude is nonzero. Chapter 3. Linear Stability Analysis a)p phase/it 66 b)w phase/it c)u phase/it d)c phase/it e)p phase/it 0.00 0.15 amplitude 0.0 0.3 amplitude 0.00 0.65 amplitude 0.00 0.26 amplitude 5 N —0 —5 1.5 0.0 amplitude f) 2 + (u ) /2 w g) uw h) wp j) i) wp uu+ww NOz: 0.00 0.02 —0.011 0.011 Right mode, cc = —0.003 0.110458, cr 0.003 —0.015 = 5 3x10 0.015—3x10 5 0.296927, and z 0.30615. a)p phase/it b)w phase/it c)u phase/it d)c phase/it e)p phase/it 0.0 1.5 amplitude 0.00 0.15 amplitude 0.0 0.3 amplitude 0.00 0.65 amplitude 0.00 0.26 amplitude N f) 2 + (u ) /2 w g) uw h) wp j) i) wp uzu+wzw 5 N —0 —5 0.00 0.02 0.011 —0.011 Left mode, ac = —0.003 0.006855, Cr 0.003 —0.015 5 3x10 O.015—3x10 5 —0.943582, and z = Figure 3.16: Eigenfunctions and correlations when Eigenfunctions and correlations when J = 0.3, and o = 0.67. E = 0.5, R = 3, Re = —1.76975. = 0.5 300, Pr = 9, H = 10, Chapter 3. Linear Stability Analysis uw 67 wp wp’ uu+Ww’ 5 0 -5 0.000 S b 0 i 4 0.035 —0.024 0.024 —0.012 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 b 4x10 0 i 4 5 0.18 _ 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 5 0 —5p 0.012 —0.18 0.18 + + Figure 3.17: Evolution of perturbation kinetic energy and correlations with = 0. Evolution of perturbation kinetic energy and correlations when both modes are present for e = 0, R = 3, Re = 300, Pr 9, H 10, J 0.3, and c = 0.54. Chapter 3. Linear Stability Analysis + 2 (u ) /2 w uw 68 wp w’p 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 4x10 —4x10 0.18 5 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x10 LO 5 4x10 (N r) L Figure 3.17 continued. 5 4x10 Chapter 3. Linear Stability Analysis 69 uu+ww a 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x 5 0 —511 0.000 c’1 II - Figure 3.18: Evolution of perturbation kinetic energy and correlations with 5 4x10 0.25 Evolution of perturbation kinetic energy and correlations when both modes are present for e = 0.25, R = 3, Re 300, Pr 9, H = 10, J = 0.3, and a = 0.52. Chapter 3. Linear Stability Analysis uw 70 wp wp uu+wW - 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 0.18 —4x10 4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x 0.000 LO - Figure 3.18 continued. 4x1 Chapter 3. Linear Stability Analysis uw 71 wp uu+WW 0 II 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 -5 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.000 5 II 0 >< : 0.000 5 N) 0 II —5 0.000 I 5 4x10 5 0.18 —4x10 Figure 3.19: Evolution of perturbation kinetic energy and correlations with e = 0.5 Evolution of perturbation kinetic energy and correlations when both modes are present for e = 0.5, R = 3, Re = 300, Pr = 9, if 10, J = 0.3, and a = 0.67. Chapter 3. Linear Stability Analysis + (u ) 2 /2 w u,w 72 w ‘p wp uru+ww 5 0 —5 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 4x10 —4x10 0.18 5 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 4x10 5 0.18 —4x10 LI) r) 0HE 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18 —4x 0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 5 0.18 —4x10 4x10 N. Figure 3.19 continued. Chapter 4 Nonlinear Numerical Simulations In this chapter we wish to study the nonlinear evolution of perturbations in a flow with initial background profiles given by (3.20). When variations in density are small, (2.17) and (2.18), derived in chapter 2, govern the two dimensional evolution of the flow: (A) + u() + w() = Jp + (4.1) and Pt + Px + WPz (4.2) RePr Since we wish to calculate the perturbations from the initial mean state, we write b(x,z,t) = p(x,z,t) u(x, z, t) Pa(Z)+P’(X,Z,t), = w(x, z, t) where /“, p’, u’ aild (z,0) + ‘(x,z,t), U(z) + u’(x, z, w’(x, z, t), w’ are the perturbations from the initial mean state. Substituting (4.3) into (4.1) and (4.2) we obtain the following equations to be solved. (‘) + (U + u’)(A’) + w’(U + (‘)) = Jp + (AA’ + (4.4) and p + (U + u’)p + w’(pa + P’)z = RePr + pazz). (4.5) It is assumed that the flow is periodic in the direction of the flow (x-direction) with = = p’ = 0 at the horizontal boundaries z = +H/2. 73 Chapter 4. Nonlinear Numerical Simulations 74 The initial perturbations are of the form (sb’, p’) = J{(q, ) exp(i&x)}, where d is the real wave number which determines the period of the flow, and (q, ) are the complex amplitudes determined using the normal modes approach to linear stability analysis. As discussed in chapter 3, the most unstable modes for the right and left moving waves occur at different wave numbers when e Re 300, Pr = 9, and H = 0 (see figure 4.1 for the case J = 0.3, R = 3, 10). If these modes were used to initialize the flow, the computational domain would have to be very large in order to maintain periodicity in the x-direction. Since it is believed that it is the linear growth rate that determines the relative stability of the modes, we use the following criteria for selecting the wave number, &, of the flow. We choose & so that the ratio of the left to right growth rates at this wave number equals the ratio of the maximum left growth rate to the maximum right growth rate. The values of these are indicated in figure 4.1 for the above mentioned case. The growth rates at the selected wave numbers are close to the maximum values. Although both modes have the same period in the x-direction, the complex phase speed, c, and complex amplitudes, (, ), are different for left and right moving waves. The precise form of the initial perturbations were discussed in section 3.3.4. The initial amplitudes of the perturbations are determined by setting the initial amount of wave kinetic energy in a given mode, K’, sufficiently small (0.0005 for the simulations examined here) to guarantee that the nonlinear simulation spends some time in the linear regime. The wave kinetic energy K’ is defined as K’(t) — (u/2 + (4.6) w2). Since we are calculating the perturbations from the initial mean flow we write u(x, z, t) ü(z, t) + u’(x, z, 1) where ü(z, 1) Here () = U(z) + iZ(z, t), and u’(x, z, t) = u’(x, z, t) represents the change in the mean flow due to diffusion. (j = LH/2 = — = F(z, t). (j dx)/L and •dz indicate horizontal averaging and vertical integration respectively, where Chapter 4. Nonlinear Numerical Simulations L = 75 2K/& is the length of the domain in the x-direction. Equations (3.3) indicate that the mean velocity and density profiles diffuse over time due to the presence of viscosity and thermal diffusion. This means that the values of the bulk Richardson number J and the Reynolds number, Re, will increase over time (due to the increase in the shear layer thickness, Ii). Hence, when the simulation enters the nonlinear regime, the actual values of J and Re will have increased from their initial value. We will continue to describe, however, the flow by the initial values of these parameters. 4.1 Numerical Method In order to numerically solve (4.4) and (4.5) we use second order centered differencing in space and a fourth order Rnnge-Kutta scheme in time. We let Ax and Az denote the step sizes in the x- and z-directions respectively. Using the following centered differencing approximations — “-i ()i+,i — Dx 2Ax D D 2 8z — ()+i — — 2Az (.)j+lj — . — D = 2(.) 2 (Ax) tJJrhL 2(.) + 2 (Az) . (‘)i,j—l + we can write the semi-discretization of (4.4) as = —(Ui + (4 7) + D’J4(U, + DAhW,) + Chapter 4. Nonlinear Numerical Simulations + = F(’I”, ‘, 76 + U) (4.8) and (4.5) as = (Uj + + Dx,j(pazj + Dze,) + + RePr’ = G(1’, e’, U, p),j (4.9) ‘(iAx,jAz,t) and where p’(iAx,jLz,t) are the approximate values of the stream function and density perturbations at a given grid point, and () indicates differentiation with respect to time. Before discretizing in time, (4.8) must be written in the form = f(’, e’, (4.10) To this end, we use a method described by Strang [59] which makes use of periodicity of the function in the x-direction. The first step is to take the discrete Fourier transform of the semi-discretization (4.8) in the x-direction. Writing .F{.} = , we obtain the following. + { 2 cos(2kAx) -2 ()2 } + = , e’, U)k, (4.11) Now for each Fourier component, k, (4.11) describes a tridiagonal system. This is easily solved, giving = )(‘, ‘, U) (4.12) Taking the inverse discrete Fourier transform of the above yields an equation of the desired form (4.10) which along with (4.9) can be discretized in time using the fourth order Runge-Kutta scheme. Chapter 4. Nonlinear Numerical Simulations 4.1.1 77 Numerical Stability In this section we discuss the stability of the numerical scheme described above. Canuto et al. [6] give the stability region of the fourth order Runge-Kutta scheme used for the time discretization in the numerical solution of (4.4) and (4.5). A schematic of the region is shown in figure 4.2. In order to have stability, we require that the amplification factor of the semi-discretized scheme (4.8) and (4.9) lie inside the kidney-shaped region. For simplicity, we will require that the amplification factor, A, lie within the indicated rectangle, i.e., o ?{A} (4.13) o < {A} < where At is the size of the time step. To determine an approximate range for the amplification factor A, we first consider the original equations (4.4) and (4.5). Since these are nonlinear equations in ‘ and p’, we use the method of frozen coefficients to linearize the problem. We know that U 1. Also, since the velocity perturbations u’ and w’ are much smaller than the mean flow we can use u’ < 1 and w’ < 1 (note that these bounds are much larger than the physical bounds). Using the above bounds, we can reduce (4.4) and (4.5) to the following linear system. 2 (A’) + 2(A’) + U + (A’) Jp + (AA’ + U ), 2 (4.14) and p + 2p + Paz + P RePr + pazz). (4.15) Since the quantities U, U, Paz, and Pazz are independent of the grid size, we can restrict our attention to the stability of the following equations. (A’) = —2(A’) — (A’) + Jp + AA’, (4.16) Chapter 4. Nonlinear Numerical Simulations 78 and = —2p + (4.17) RePA Using the notation described by (4.7), the semi-discretization of (4.16) and (4.17) is = (4.18) + + — and = —2D, + — (4.19) RePr’ We use von Neumann analysis to determine the amplification factors of (4.18) and (4.19). Taking the two-dimensional discrete Fourier transform of (4.18) and (4.19) we obtain = [(Ah)2 + hk,l iAhD] ,l + iJ-k,l (4.20) ek,l where = 27rkx, 2 + 2(cos 1)/(Lx) = — = RePr 2irlz.z, 15 +D gk,l —2sin/Ax = — sin/Az, and /q, = 2(cos — , the discrete Fourier transform of the Laplacian operator. 2 1)/(L\z) We can write (4.20) in matrix form, -‘ . +zD sine / I 1hAX k,l (4.21) ek,l 0 RePr The amplification factors of the semi-discretization are given by the eigenvalues of the above matrix, ) = 2 Ah/Re + iD, and A = Lh/RePr + iD. The restrictions (4.13) give the following conditions for stability. , = 2 (cos—1 cosC—1” + Re 2 (x) 2 (z) ) —4 Re ( 1 (Ax)2 + 1__\\ )2} —2.79 2t 4.22 and {A}= —2sin sinC Lx Az Ax AzAt’ (4.23) Chapter 4. Nonlinear Numerical Simulations 79 providing Pr> 1 (and hence {A } > 1 From (4.22) and (4.23) we can deduce the following restriction on the time step, At, I At < mm < — I_L+J_’o tz x 1>. 2.79Re 1 1 (4.24) Although (4.24) only guarantees a stable numerical scheme for the linear equations (4.14) and (4.15) with constant coefficients, in general it is still used to determine an appropriate value of At for the numerical solution to the nonlinear equations (4.4) and (4.5). Results 4.2 4.2.1 Methodology For the background flow defined by (3.20), linear theory predicts that the right mode dominates when > 0. We start our investigation of the nonlinear evolution by examining the development of a. flow when the initia’ perturbation contains only the right mode. Next, we perturb the flow with the mode corresponding to the more slowly growing left moving wave. Finally, in order to examine the interaction between the two waves, we impose both modes. To examine the evolution of the perturbations, we examine the density and vorticity distributions. At different times throughout the simulation, we plot contours of constant density ranging from -0.72 in the upper layer to 0.72 in the lower layer. The contours are of equal increments of 0.18. We superimpose vorticity (w = — w) shading upon the density contours, where white corresponds to -1.8 (counter-clockwise vorticity) and black to 1.8 (clockwise vorticity). These plots enable us to visualize how the perturbations evolve. They also allow us to qualitatively compare results with those of laboratory experiments. It is desirable to compare results of the nonlinear simulation with those of linear Chapter 4. Nonlinear Numerical Simulations 80 stability analysis. Specifically, we wish to compute the growth rates and phase speeds of the perturbations. In appendix D, we showed that K’(t) oscillates in time when two modes are present. A lack of oscillations would indicate that the perturbation is composed of one mode only. To compute the growth rate of the simulations, any oscillations in K’(t) are removed using a low-pass filter. The linear growth rate ac can then be determined from the resulting perturbation kinetic energy (see Smyth [51]) K’. ac To determine the phase speed cr, = ld j-1nK’. (4.25) we track the position of maximum/minimum values of density at the critical Jevels determined from linear analysis of the right/left modes. The change in position with respect to time gives the phase speed. This method of tracking the right and left moving waves has been used by Smyth et al. [53] in their study of symmetric Holmboe waves. In appendix C, the equations governing the evolution of mean kinetic energy, pertur bation kinetic energy and potential energy were derived and are given by (C.22), (C.23), and (C.29), respectively. Adopting the notation of Klaassen & Peltier [29], we have = —C(K, K’) — D(K) C(k, K’) + C(P, K’) = —C(P, K’) — — D(K’) D(P) (4.26) (4.27) (4.28) where C(k,K’) —(u7), D(k) C(F,K’) -J(), D(K’) = D(P) = = ((u) + ()2 + (w)2 + (w)2), RePr (4.29) Chapter 4. Nonlinear Numerical Simulations 81 The notation C(R ,R 1 ) indicates a conversion of energy from reservoir R 2 1 to reservoir . D(R 2 R ) 1 > 0 represents a loss of energy from reservoir R 1 due to diffusion. We shall examine the evolution of the three energy reservoirs, K, K’, and P over time as well as the total energy T = K + K’ + P. In addition, we will examine the terms that govern the evolution of energy, i.e., the exchange of energy between reservoirs and the loss of energy due to dissipation. Equations (3.3) indicate that the mean velocity and density profiles diffuse over time due to the presence of viscosity and thermal diffusion. It is of interest to see how the presence of perturbations in the flow affect the development of the mean flow. To this end we compare the evolution of the computed mean flow ü(z, t), and ,ö(z, t) (where (z,t) pa(Z) +‘(z,t)) with the exact solutions to (3.3). Solutions to (3.3) with initial conditions given by (3.20) and boundary conditions H Tn H’ — sech2(2) — ‘H 1) U4— 2’ ‘ — — — — secni2l’H j 2 ,j — — (4.30) —tau1hR(—*+) pa(,t) pa(,t) Pa_ —tanhR(*+6) = Pa÷ = can be found using separation of variables (see, for example, Boyce DiPrima [2]). They are given by ü(z,t) = 0 A °° 7 I—n t 2 1 r frz 7r\ Uz+--+Aexp 1 1c0snff+), 2 ReH 00 (z, t) = a + PaZ B exp + { 2 :P:H2 } ( + ), (4.31) (4.32) where A 2 = H I Id/2 {U(z) H/2 — uz} cos n f { 1 Pa(Z)_(a irz + 7r +Apaz)}sinn dz, (+ ) dz, (4.33) (4.34) Chapter 4. Nonlinear Numerical Simulations with = (c-- + pa_)/2 and LPa = (f3a — 82 pa_)/H. Truncating the above series (4.31) and (4.32) at n = 15 is an excellent approximation for t > 100. 4.2.2 Parameter Values As already discussed in chapter 3, Smyth [51] found that R> 2.4 is required in order for Holmboe instabilities to exist. Since higher values of R mean a thinner density interface, more grid points are required in the vertical direction in order to numerically resolve the flow. We choose to study the case R = 3 which is sufficiently large to support Holmboe instabilities yet small enough to allow computations without requiring too fine a grid. This corresponds to Pr = 9 which is the approximate value for the diffusion of heat in water at 10°C (table 3.1). In the numerical simulations described we set Re = 300 which is sufficiently large to allow the instabilities to fully develop before they are damped by viscosity. Furthermore, it is of the same order of magnitude as that of many laboratory experiments. Smyth et al. [53] have examined in detail the nonlinear evolution of flolmboe waves when e = 0 with J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. We choose the same parameter values in our study of non-symmetric flolmboe instabilities. This allows us to verify our results for e = 0 against those of Smyth et al. It also gives us a good basis for comparison for the non-symmetric case. J = 0.3 corresponds approximately to the most unstable Holmboe instability (figure 3.11). Also, figure 3.12 indicates that by setting H = 10 we avoid the influence of the horizontal boundaries on the development of the instabilities. We examine the cases e = 0, 0.1, 0.25, and 0.5. To determine the appropriate number of grid points we compute the evolution of the perturbed flow at several grid levels and compare the results. As discussed in sec tion 4.1.1, the time steps selected satisfy (4.24). Perturbation kinetic energies are shown Chapter 4. Nonlinear Numerical Simulations 83 in figures 4.3, 4.4, and 4.5 for the right mode, left mode, and both mode cases, respec tively. Since the numerical scheme used is dispersive, it is not surprising that, for e = 0 and 0.1, there is a phase shift between the results computed on the nx=nz=64 grid, where nx and nz are the number of grid points in the x- and z-directions, respectively, and those from the nx=nz= 128 grid. Comparing density contours (not shown), we found that, although the phase speeds differed between the two grid levels, the flows looked qualitatively the same. Since there is little phase difference between the 128 by 128 and 256 by 256 grids used for e = 0.25 and 0.5, we believe that we can compute on a 128 by 128 grid with confidence up to t When = 400 for e = 0 and 0.1. 0.25 and 0.5 we encountered difficulties in numerically resolving the flow. This appeared to be associated with the extreme gradients in the density field and is demonstrated in figure 4.6. The perturbation in the density gets very steep. Denser fluid from the lower layer separates from the wave and is released into the less dense upper layer. For e 0.5, this results in a decrease in the perturbation kinetic energy (see, for example, figure 4.5 for E = 0.5) and occurs at different times for the nx=nz=128 and nx=nz=256 cases. Although this type of behaviour has been observed in experimental studies, we were not able to resolve the exact time or manner in which it occurs. It was observed that this ejection of denser fluid into the upper layer is most likely to occur when e is large. This is a result of the rapid growth rate of the right mode and explains why we have resolution problems when When e = = 0.5. 0.25, ejection of denser fluid into the upper layer only occurs on the coarser grid when the initial perturbations only contain one mode. We also notice that there is no sudden drop in the perturbation kinetic energy. Since the density contours remain qualitatively the same, except during the brief period of fluid breaking off from the perturbation, and the perturbation kinetic energy compares favourably between the two grid levels, we believe that we can compute with confidence on the finer grid up Chapter 4. Nonlinear Numerical Simulations to t = 84 400, when oniy one mode is initially present. When the flow is perturbed with both modes, however, fluid breaks off from the perturbation at both grid levels, but at different times. Although the perturbation kinetic energies do not differ significantly, the density contours quickly loose their similarities. Therefore we only have confidence in our results up until t = 170. A summary of the parameters used is given in table 4.1. e 0 0.1 0.25 0.5 c 0.54 0.53 0.52 0.67 lix At liZ 128 128 256 256 128 128 256 256 0.0625 0.0625 0.03125 0.03125 right mode 400 400 400 140 tmax left mode 400 400 400 240 both modes 400 400 170 140 Table 4.1: Parameters used for the numerical simulations when J 0.3, R = 3, Pr 9, Re = 300, and H = 10. The time step size is At. t rnax is the largest time for a given simulation for which the results are reliable. 4.2.3 Right Mode When the flow is initialized with the right mode only, we expect to see the evolution of a perturbation protruding into the upper layer and moving to the right. This is indeed the case as seen by the density contours shown in figures 4.7, 4.8, and 4.9 for 230 (e = 0), 190 < t 205 ( 0.25), and 125 t 140 ( t 245 0.5), respectively. This is further supported by the lack of oscillations in the initial evolution of the wave kinetic energy (figure 4.3). When e = 0, however, oscillations begin to develop at t indicates the presence of another mode. Figure 4.7, 385 200. This t < 400, shows that it is the left moving wave that is developing. This is not surprising as left and right modes have the same linear growth rate for this case and one expects the left wave to eventually obtain the same amplitude as the right wave. Another interesting observation is that, Chapter 4. Nonlinear Numerical Simulations 85 even though right and left waves have the same linear growth rates, once the left mode starts to grow it is in the nonlinear regime and grows much more slowly than if it had started its growth in the linear regime. As e increases, the right mode dominates for a longer amount of time before the development of the left moving wave. This is due to the reduced growth rate of the left mode as well as the increased growth rate of the right mode. An interesting phenomenon is the tilting of the wave in the opposite direction to which it is travelling (figures 4.7 and 4.8). This tilting is caused by a concentration of counter-clockwise vorticity that develops at the wave crest and has been observed in laboratory experiments by Lawrence [37]. The movement of the wave to the right is closely linked to the concentration of clockwise vorticity that proceeds the perturbation in the density interface and also moves to the right. Ripples in the vorticity to the left of the wave are a resolution problem caused by the dispersive nature of the numerical scheme and become less noticeable as the grid is refined. We calculate the effective linear growth rate of the perturbations. Since only the right mode is present in the initial development of the perturbations, no oscillations are observed in the perturbation kinetic energy. Thus we can compute the growth rate directly without filtering. In figure 4.10 these are compared to the growth rates predicted by linear stability analysis. When theory and the simulations. As = 0 we have excellent agreement between linear increases, the growth rate of the simulation is less than that of linear theory. The maximum amplitude achieved appears independent of the value of (figure 4.3). Consequently, the perturbations spend less time in the linear regime as e increases becailse of the increasing growth rate of the right mode. Holmboe [23] predicted from linear theory, and Smyth et al. [53] showed numerically, that the phase speeds of the two modes in a Holmboe wave are not constant but move fastest when they pass each other and slowest when they are widely separated. When Chapter 4. Nonlinear Numerical Simulations 86 the flow is perturbed with the right mode only, this oscillation in the phase speed is not initially observed since the left wave is not present to interact with the right wave (figure 4.11). In all cases, the perturbations at both critical levels are initially moving to the right. At later times, however, the perturbation moving at the lower critical level starts to exhibit oscillations in its phase speed. This coincides with the development of the left moving wave, as indicated by the oscillations in the perturbation kinetic energy (figure 4.3), and occurs at later times as e increases. When the lower critical level behaves erratically between 330 < t = 0, the perturbation at 370 after which it appears to be moving to the left. At this point, the phase speed of the right moving wave also has oscillations in its phase speed as the two waves pass through each other. In order to compare the initial phase speed with that predicted by linear theory, we compute the average phase speed over a cycle. Differences between predicted and computed phase speeds are shown in figure 4.12. In general, the initial phase speed is slightly faster than that predicted by linear theory. Also, the phase speed increases steadily with time. It appears that the phase speed may be asymptoting to a constant value, although further investigations are required to verify this conjecture. By studying the energy budgets given by (4.26), (4.27), and (4.28), we can gain a better understanding of the mechanism that controls the growth of the instability. We will start by discussing the trends that are present for all cases. First, we notice that the total energy decreases with time (figures 4.13 and 4.14) and that this is caused by the diffusion of the mean flow. In addition, diffusion of the mean flow causes the potential energy P to increase over time. During the initial linear growth of the instabilities, the decay in mean kinetic energy, K is primarily due to diffusion. As the perturbations enter the nonlinear regime, however, we see a decrease in the mean diffusion term D(K) and an increase in C(K, K’). At this point, the growth of the perturbation has a more significant effect on the total energy budget. Chapter 4. Nonlinear Numerical Simulations 87 Since the flow was perturbed with oniy the right mode, we initially see no oscillations in any of the quantities in the energy budget. During this period, C(K, K’) > 0 and C(P,K) <0 with C(f(,K’) > C(P,K’), indicating a transfer of mean kinetic energy to perturbation kinetic energy and a smaller transfer of perturbation kinetic energy to potential energy. Thus the perturbations are undergoing a net growth as indicated by the increase in K’. For = 0 oscillations occur at later times in the quantities K, K’, F, C(k,K’), C(P,K’), D(]?) and D(K’). As discussed above, this is due to the eventual growth of the left mode and indicates that the two waves are interacting. Notice that C(K, K’) remains almost always positive and hence the perturbation is constantly ex tracting energy from the mean flow. Conversely, C(F, K’) oscillates between positive and negative values. Therefore perturbation kinetic energy is converted to potential energy for some time and then the process is reversed. The cumulative transfer f C(F, K’)dt is always negative, however, and so there is a net transfer of perturbation kinetic energy to potential energy. When the left mode becomes sufficiently large, there is a drop in D(K) which is reflected by a flattening of the mean kinetic energy curve. A more detailed description of the energy budget when both waves are present will be discussed for the = 0 case initially perturbed with both modes as this permits easy tracking of the two waves. As e increases, the oscillations occur at later times. For = 0.5 (figure 4.14), there are no rapid oscillations in the short time for which the simulation was resolved. The behaviour is very different from the e = 0 case. We now see oscillations of much longer periods in C(]?, K’) and C(P, K’). This is reminiscent of the energy budgets observed by Klaassen &i Peltier [29] for the Kelvin-Helmholtz instability. One significant difference is that, in Klaassen & Peltier [29], C(]?, K’) and C(F, K’) oscillate between positive and negative values, whereas, in our case, C(K, K’) remains positive and C(P, K’) is negative for t < 130. Since the Kelvin-Helmholtz instability is characterized by the presence of Chapter 4. Nonlinear Numerical Simulations 88 interwinding fingers of heavier and lighter fluid, Klaassell & Peltier [29] attribute the transfer of potential energy to mean kinetic energy (via the perturbation kinetic energy) to the fingers of heavier fluid being pulled back down and the figures of lighter fluid being lifted upward. In the flow examined here, the perturbation is always displacing denser fluid upward and so there is no transfer of potential energy to mean kinetic energy. When 130 < t < 140, C(K, K’) is negative. During this time, there is a decrease in the wave kinetic energy. The wave amplitude, as seen by the density contours (figure 4.9), also decreases. The transfer quantities C(K, K’) and C(P, K’) are in phase and thus there is a continuous transfer of energy from mean kinetic energy to perturbation kinetic energy to potential energy. When characteristics of both the E = 0.1 and 0.25 (not shown) the energy budgets show 0 and e = 0.5 cases. Eventually rapid oscillations appear as the left wave grows, but there are more slowly varying oscillations superimposed upon these. One of the main concerns in the study of the evolution of instabilities is to determine how they influence the development of the mean flow. Mean velocity and density profiles at the final time in the numerical simulation and those given by (4.31) and (4.32) for the evolution of the mean flow with no perturbations are shown in figure 4.15. Thus, we are able to analyse how the presence of the perturbations affects the development of the mean flow. When e = 0, the mixing of momentum (as seen from the mean velocity profile) is primarily in the upper layer. As e increases, the amount of mixing in the upper layer does not change significantly, but it spreads into the lower layer. Mixing of momentum occurs throughout the upper layer, whereas changes in density due to mixing are concentrated near the crest of the wave with a smaller change occurring near the density interface. As e increases, the change in density due to mixing increases at both these locations. Chapter 4. Nonlinear Numerical Simulations 4.2.4 89 Left Mode Except for the 0 case, the flow behaviour is markedly different when it is initialized with the left mode only than when it is initialized with the right mode only. This is a direct consequence of the right mode having the largest linear growth rate when e > 0. As the difference between the growth rates of the left and right modes increases, i.e., as increases, the right mode influences the flow much more quickly and results in the onset of oscillations in the wave kinetic energy occurring at earlier times as (figure 4.4). In fact, for = increases 0.25 and 0.5, the right mode affects the flow while it is still in the linear regime. Furthermore, density contours when e = 0 and 0.1 (not shown), show that only the left mode is present in the early stage of the nonlinear regime, whereas for 0.25 and 0.5 the right mode, having a larger growth rate, has had sufficient time to grow and is clearly present (figures 4.16 and 4.17) even at early times in the nonlinear regime. When = 0 the flow behaves exactly the same way as the right mode only case, except that the roles of the right and left moving waves are reversed. Results for > 0 indicate that both left and right modes will eventually develop in the flow. We also see that both left and right moving waves have a concentration of counter-clockwise vorticity in their peaks once their amplitudes become sufficiently large. When > 0 the left mode only case spends longer in the linear regime than the right mode only case due to the smaller growth rates. The linear growth rates from the nonlinear simulation are given in figure 4.10 and are only slightly larger than those predicted by linear stability analysis. As is seen from the wave kinetic energy (figure 4.4), the = 0.5 case behaves differently from the smaller values of e. It only spends a very short amount of time growing at the rate of the left mode. After a short transition period, it continues to grow linearly but with a growth rate of 0.045. As this value lies midway between the growth rates for right and left modes predicted from linear analysis, Chapter 4. Nonlinear Numerical Simulations 90 we conclude that the change in growth rate is due to the presence of the rapidly growing right mode. As with the right mode only case, the phase speed of the left mode does not oscillate initially (figure 4.18). Also, perturbations at both critical levels move to the left. It is only when the right wave affects the flow that we observe oscillations in phase speed at the critical level of the right wave. When e 0.1, the motion of the perturbation at the critical level corresponding to the right moving wave exhibits interesting behaviour (figure 4.18). Initially the wave moves to the left, as expected. As the wave evolves, oscillations in the phase speed become more dramatic until a turn about occurs at t 200. The wave then travels to the right. Notice that the wave in the lower layer is unaffected by the turn about and continues to move to the left. Similar behaviour is observed for = 0.25 and 0.5 except that the turn about occurs sooner as increases. Although only the right mode is noticeable in the density contours in figure 4.17 ( = 0.5) the fact that the perturbation at the lower critical level is still moving to the left would suggest that the left moving mode is growing and would eventually become noticeable in the flow. This is further supported by the results shown in figure 4.16 where both modes are clearly present. Figure 4.12 indicates that the initial average phase speed of the left mode is slightly smaller than that predicted by linear theory. Also, for e = 0, 0.1, and 0.25, the phase speed increases steadily over time. This effect becomes less pronounced as and reverses when E = increases 0.5 for which the phase speed decreases with time. This supports the hypothesis that the phase speeds are asymptoting to a single value. When = 0, the energy budget for the nonlinear simulation initially perturbed with the left mode only is identical to the right only case shown in figure 4.13. This result is not surprising in view of the symmetry of the flow in this case. As e increases, oscillations begin sooner since, with increasing e, the right mode becomes important more quickly. Chapter 4. Nonlinear Numerical Simulations 91 Because the difference in linear growth rates between right and left modes increases with increasing , the interactions between the two waves is weaker for larger values of e. This is probably due to the difference in amplitudes of the two waves for the times considered and is evident in the reduced amplitude of oscillation in C(P, K’). In general, however, the energy budgets look similar to the When E = 0 case. 0, momentum mixing, as indicated by the difference between the evolution of the mean flow with and without perturbations, is the same as the right only case (figure 4.15) with the amounts of mixing in upper and lower layers reversed. As e increases, there is a decrease in the amount of mixing in the lower layer and and increase in the amount of mixing in the upper layer. To explain this phenomenon, we examine how the growth rates of right and left modes change with increasing e. Since the growth rate of the left mode decreases as increases, it grows less as increases which results in less mixing in the lower layer. A similar argument holds in the upper layer which is influenced by the right moving wave. The right wave causes more mixing with increasing e since it is the dominant mode when e> 0. The change in density due to mixing exhibits the same trends observed for the momentum mixing, although, as with the right mode only case, the magnitude of the change in density is much smaller than that of the change in the mean velocity field. 4.2.5 Both Modes For the symmetric case (e the wave kinetic energy = 0) when the flow is initialized with both modes, growth of (figure 4.5) is similar to those of the left and right mode only cases except that oscillations are present at all times. These oscillations are expected since both modes are present from the start and grow at the same rate. The flow looks similar to that computed by Smyth et al. [53] with a right moving wave protruding into the upper layer and a left moving wave protruding into the lower layer. One interesting Chapter 4. Nonlinear Numerical Simulations 92 difference is that, as seen in figure 4.19, our simulations have not lost the symmetry at later times whereas those of Smyth et al. [53] did. We expect the flow to remain symmetric even at large times since the background flow is initially symmetric, as are the horizontal boundary conditions. Results when the flow is initially perturbed with only one mode further supports the results observed here. It was seen that, when e = 0, the flow tends towards a state where both modes are present. As e increases, oscillations in the perturbation kinetic energy decrease (figure 4.5), es pecially in the linear regime. This observation indicates that, although the perturbations are initialized with the same amount of kinetic energy in both left and right modes, the right mode quickly dominates the flow for large values of e. Furthermore, the effective linear growth rate lies in between the values predicted by linear theory for the right and left modes, but is dominated by that of the right mode (figure 4.10). When e = 0 we have excellent agreement between the predicted linear growth rate and that of the linear simulation. Even for relatively small values of e, i.e. e = 0.1, linear theory predicts a difference in the growth rates for left and right modes (figure 4.1). This has a significant effect on the development of the perturbations. The right mode grows more quickly and initially has a larger amplitude (figure 4.20, 195 t 210). The left mode, however, continues to grow and is almost as large as the right mode at the end of the simulation. The numerical simulation is starting to have resolution problems at later times. In order to compute for longer times, a finer grid should be used. When 0.25, there is a rapid growth of the right mode which dominates the flow in the early development of the perturbations (figure 4.21, 100 t 115). Since the flow is initialized with both modes, the left moving wave is still present but, due to its smaller growth rate, takes longer to develop. At t = 170 (figure 4.21), we observe both right and left modes, although the right mode still has the larger amplitude. When e = 0.5, the Chapter 4. Nonlinear Numerical Simulations 93 flow initially evolves in a similar manner (figure 4.22) with only the right mode present. In fact, for the time in which we have confidence in our results, the flow looks similar to the right mode only case. This is probably a direct result of the reduced growth rate of the left mode at larger values of e. Although the left mode is still present, it has not yet had a chance to grow to an observable level for the results presented here. We discuss the behaviour of the flow for t > 140 when E = 0.5. This is beyond the range of time for which the flow was adequately resolved, but is still instructive from a qualitative point of view. The right mode continues to grow until t 185. At this time, the protrusion into the upper layer separates from the density interface and diffuses into the ambient fluid (figure 4.6), causing a drop in the perturbation kinetic energy (figure 4.5). After this sudden decrease in magnitude of the perturbation, the growth of the right mode appears to be stunted. The left mode, however, continues its slow growth and starts to dominate the flow (figure 4.22). Although we were not able to resolve this behaviour numerically, it occurred in several of the simulations. We believe that the numerical simulation helps us understand the overall mechanism of the instability for large values of 6 even though the precise details were not resolved. Furthermore, Koop et al. [34] have observed a similar phenomenon in laboratory experiments. Their flow is initially unstable to Kelvin-Helmholtz type instabilities. The billows formed mix quickly into the ambient fluid. Further downstream, weaker instabilities of a Holmboe nature appear at the interface. This leads us to believe that the phenomena observed in our simulations are physically possible. The positions of the perturbations over time in both top and bottom layers are shown in figure 4.23. For 6 = 0 our results are in agreement with those observed by Smyth et al. [53] and predicted by Holmboe [23] for the linear regime. During this stage, the waves speed up as they pass through each other, and slow down when they are far apart. As the perturbations enter the nonlinear regime, however, there is a shift in the behaviour of Chapter 4. Nonlinear Numerical Simulations 94 the phase speed. The phase speed still fluctuates during a cycle, but now the waves move fastest as they approach each other, slowing down as they pass through each other. The method for tracking the phase speed breaks down as the waves start to interact in the nonlinear regime. Just after the waves pass through each other, there is a jump in the position of maximum/minimum density (figure 4.23). This jump occurs when density contours within the wave have more than one local extremum (figure 4.19, t 245), and indicates an inaccurate measurement of the instantaneous phase speed. The average phase speed throughout the cycle appears to increase with time (figure 4.12). This observation is consistent with the results of the right mode only and left mode only cases, but the interaction between the two waves slows down their average phase speed. In light of the above discussion, further investigations are required to better resolve this behaviour. When = 0.1, the flow initially behaves the same as the e = 0 case. It has a right moving perturbation in the upper layer and a left moving perturbation in the lower layer. It also exhibits oscillations in the phase speed as the two waves interact. For 40 <t <100, the right moving wave appears to he unaffected by the presence of the left moving wave and does not have any oscillations in its phase speed throughout the cycle. On the other hand, the left moving wave has large changes in its phase speed throughout the cycle. This is most likely a result of the difference in amplitude of the two modes at this stage of their development. The right mode has had sufficient time to grow and is clearly present in the density contours, whereas the left mode has not yet grown enough to be seen (not shown). After t for e = = 100 the two waves behave as in the nonlinear regime 0, speeding up as they approach each other. For t> 330, the two waves appear to be moving in a non-smooth way. This behaviour is due to the resolution problem already discussed above. As for the = 0 case, the average phase speed of the right moving wave increases over time (figure 4.12) and is initially the same as for the right mode oniy case. Chapter 4. Nonlinear Numerical Simulations 95 As the left wave grows, however, its presence starts to affect the phase speed of the right wave which is slower than the right mode only case. The phase speed of the left wave initially increases but then decreases as the influence of the right wave strengthens. For larger values of (i.e., E = 0.25 and 0.5), perturbations in the upper and lower layers exhibit the same characteristics in the initial linear regime as the two previous cases. Oscillations in the phase speeds as the waves pass through each other are only present for a very short time, however, after which the right moving wave behaves very much like the case of the right wave only, with no oscillations in the phase speed. The perturbation at the left critical layer looses it distinctive motion to the left. It bounces between moving to the left and moving to the right in no particular fashion, indicatillg a complicated interaction between the slowly growing left moving wave and the rapidly growing right moving wave. For t > 140 in the e = 0.25 case, the left moving wave has grown large enough for the two waves to interact in the nonlinear regime as described above. The reason it takes longer to reach this state than the e = 0.1 case is the left mode has a much slower growth rate and hence takes longer to achieve a large enough amplitude to affect the motion of the right mode. For e = 0.5, the simulation does not last long enough for the left mode to grow sufficiently large to affect the evolution of the right mode, but the above results suggest that the flow would eventually exhibit the same behaviour as the other cases with smaller values of e. For both r 0.25 and 0.5, the average phase speed of the right mode increases steadily over time (figure 4.12) and appears unaffected by the left moving wave. The phase speed of the left mode, when = 0.25 appears to be slowed down by the presence of the right wave when compared to the left only results. The simulation for E — 0.5 is too short to draw any conclusions about how the phase speed of the left mode varies with time. The energy budget for e = 0 is shown in figure 4.24 and was discussed in detail by Smyth et al. [53]. We will review the results here and discuss how they tie into the Chapter 4. Nonlinear Numerical Simulations 96 previous cases. As before, we see that the total energy T decreases with time due to diffusion. Although there are rapid oscillations in K, K’, and F, there are no oscillations in the total energy indicating a balance in the exchange of energy between the reservoirs. The long term behaviour of K, K’, and P are as in the right mode and left mode only cases. Since j C(k, K’)dt > 0 for any value of t, there is a net transfer of energy from the mean flow to the perturbations. The greatest transfer of energy, however, is between the potential energy and the perturbation kinetic energy. Although C(P, K’) oscillates between positive and negative quantities, we have that J C(P, K’)dt < 0 at any instant. Thus there is a net conversion of perturbation kinetic energy to potential energy. Nevertheless, f C(K, K’)dt > — J’ C(P, K’)dt indicating that a portion of the energy transferred to the perturbation kinetic energy from the mean flow remains in the perturbation reservoir. Diffusion of the mean kinetic energy decreases as K decreases. Also, since both modes are present at all times, there is no sudden drop in D(K), as when the flow is initially perturbed with one mode only, but a steady decrease over time. Initially D(K) > G(K, K’) but as the flow enters the nonlinear regime, the two quantities have the same order of magnitude and both contribute to the overall decrease in K. Diffusion of the perturbation kinetic energy increases as K’ increases whereas D(P) remains constant throughout the simulation. From figure 4.23 we can determine when the right and left wave pass through each other and when they are maximally separated (for example, figure 4.19, at t and t 235 240, respectively). Studying the energy budgets we can determine how the interactions of the two waves affect the growth of the perturbations. We will examine this interaction during the linear stage and in the fully nonlinear regime (see figures 4.25 and 4.26, respectively, for enlargements of positions and energy budgets during these times). While the waves are growing linearly, the growth and decay of the perturbation kinetic energy is controlled by the transfer of mean kinetic energy to the perturbation kinetic Chapter 4. Nonlinear Numerical Simulations 97 energy, and vice versa. Extraction of energy from K is optimal just after the waves have passed each other. This corresponds to a maximum in the perturbation kinetic energy. As the waves approach each other, there is a brief period when C(K, K’) < 0, corresponding to a minimum in K’. In general, the transfer of energy is from mean kinetic energy to perturbation kinetic energy and then from perturbation kinetic energy to potential energy (or in the reverse direction). There is, however, a slight phase shift between C(K, K’) and C(P, K’). This was observed when studying the evolution of the eigenfunctions in chapter 3. Once the waves enter the nonlinear regime (figure 4.26), transfer of energy from K to K’ has two maxima. One occurs just before the waves pass through each other and the other just after. During the earlier stages of the nonlinear regime the most efficient extraction of energy from the mean flow occurs just after the waves have passed through each other. At later times, however, it occurs as the two waves are approaching each other. When the waves are maximally separated, there is a brief period when C(k, K’) < 0 and thus the perturbations are giving up energy to the mean flow. As left and right moving waves approach each other, C(P, K’) > 0. This corresponds to a growth period for K’ and a decay period for P. During the part of the cycle when the waves are moving away from each other, there is a reversal in the behaviour of C(P, K’) corresponding to a decrease in the perturbation kinetic energy and an increase in the potential energy. In contrast to the linear regime where oscillations in K’ are governed by C(K, K’), in the nonlinear regime they are linked to C(P, K’). Also, in the linear regime, maxima in K’ occur just after the waves have passed through each other, minima occur as the waves are approaching each other. In the nonlinear regime, however, maxima in K’ occur as the waves pass through each other whereas minima occur when the waves are maximally separated. When E = 0.1 and 0.25 (not shown), the energy budgets look similar to those of the Chapter 4. Nonlinear Numerical Simulations e 98 0 case. Therefore, the two waves need not have equal amplitudes for the nonlin ear interaction to exhibit the same behaviour described above. If e is sufficiently large, however, the energy budget can look very different even when the flow is initially per turbed with both right and left modes, as seen by examining the energy budget when = 0.5 case. The 0.5 (figure 4.27) looks similar to that of the right mode only (figure 4.14). This similarity is not surprising as the density contours (figures 4.9 and 4.22) are virtually identical for the two cases. One difference between the two energy budgets is that when the flow is initialized with both modes, there are small oscillations in the transfers between the reservoirs, C(]?, K’) and C(P, K’). These oscillations indicate that the left mode, however small it may be, is interacting with the right mode. Initially, these interactions are very weak and we see a continuous transfer of energy from K to K’ with a smaller transfer from K’ to P. With increasing time, C(P, K’) oscillates between positive and negative values. This behaviour in the energy budgets looks similar to that of the right mode only when = 0.1 and 0.25, suggesting that, given enough time, the left mode will grow and the energy budget will behave as those described above when the amplitudes of the two waves are similar. Mean velocity and density profiles are shown in figure 4.28. When E = 0, the changes due to mixing in the mean velocity and density profiles are the same in both upper and lower layers. This is as expected since the perturbations remain symmetric even near the end of the simulation (figure 4.19). As observed earlier, momentum mixing tends to occur throughout the entire domain whereas changes in density due to mixing are concentrated near the waves’ crests. When the perturbations initially contained the right or left mode only, there was a small change of density due to mixing near the center of the density interface. Changes were in opposite directions for right and left modes. When both modes are present, there is no change in density due to mixing at the density interface. Hence, any changes mean density at the interface are due to diffusion and changes Chapter 4. Nonlinear Numerical Simulations 99 in density due to mixing for left and right modes cancel at this level. As E increases, the amount of momentum mixing increases in the upper layer and decreases in the lower layer. This result is consistent with results of the right and left mode only cases. For the profiles shown it is difficult to determine the effect of varying E on the density changes due to mixing. Based on previous results, we expect the same type of behaviour as for the momentum. It also appears that as e increases, there is a larger change in density due to mixing at the interface. When is larger, the left mode is not sufficiently strong to produce an opposing change in density at the interface and thus there is no cancellation of the change in density due to the right mode at this level. Chapter 4. Nonlinear Numerical Simulations 100 0.1 0.12 0.12 0 0 0,00 0.0 0.00 0.0 1.5 1.5 a a 0.25 0.5 0.12 0.12 0 0 0.00 0.0 1.5 0.00 0.0 a 1.5 a Figure 4.1: Growth rates from linear stability analysis Growth rates from linear stability analysis for J = 0.3 with Re = 300, Pr = 9, R 3, and H 10. The solid line corresponds to the right moving wave and the dotted line to the left moving wave. When E = 0 the right and left growth rates coincide. The vertical dashed line indicates the value of ö 0.54, 0.53, 0.52, and 0.67 for E = 0, 0.1, 0.25, and 0.5, respectively, the wave number used in the simulations. Chapter 4. Nonlinear Numerical Simulations 101 2f2 2 -plane -2.79 Figure 4.2: Absolute stability region of fourth order Runge-Kutta scheme Schematic of absolute stability region of fourth order Runge-Kutta scheme. The scheme is stable inside the kidney shaped region. See Canuto et al. [6] for the precise shape of the stability boundary. Chapter 4. Nonlinear Numerical Simulations 102 8=0 :: nxnz=128, dtO.O625 nx=nz=64, dt=O.125 0 200 — 300 4C 300 40 300 40 0 300 400 8=0.1 U -2- —4 —6 nx=nz=128, dt=O.0625 —8 nx=nz=64, dt=O.125 100 0 200 8=0.25 0 -2- —4 —6 nx=nz=256, dt=O.03125 —8 nx=nz=128, dt=O.0625 I —10 100 0 200 8=0.5 U —2 —4 —6 —8 100 0 Figure 4.3: 200 Perturbation kinetic energy for right mode Perturbation kinetic energy when flow is initialized grid size for J = 0.3, R = 3, Pr = 9, Re = 300, with the right H = 10. and mode oniy with varying Chapter 4. Nonlinear Numerical Simulations 103 8=0.1 o ioo 200 300 400 =O.25 0 0 100 200 300 400 0 100 200 300 400 Figure 4.4: Perturbation kinetic energy for left mode Perturbation kinetic energy when flow is initialized with the left mode only with varying gridsizeforJ=0.3,R=3,Pr=9,Re=300, andH=10. Chapter 4. Nonlinear Numerical Simulations 104 C, i -6 :: nxnz=128, dt=O.0625 nx=nz=64, dt=o.125 100 0 200 300 40 8=0.1 C) :: —6— — : nx=nz=128, dt==O.0625 -8 nx=nz=64, dt=O.125 —Ic’ 100 0 200 300 4 C=025 0 -2- ,,.. —4 —6 nxnz=256, dt=O.03125 —8 — nx=nz=128, dt=O.0625 ci —10 0 100 200 300 4C 0 300 400 8=0.5 0 —2 —4 —6 —8 —In 0 100 200 Figure 4.5: Perturbation kinetic energy for both modes Perturbation kinetic energy when flow is initialized with both modes with varying grid sizeforJ=0.3, R=3, Pr=9, Re=300, andH=10. Chapter 4. Nonlinear Numerical Simulations 105 4 2 —2 —4 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 Figure 4.6: Effect of grid size when e = 0.2 0.4 0.6 0.8 1.0 0.5 Density contours for both modes for e = 0.5 with J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. x has been scaled with respect to its period 2ir/. Top and bottom rows are results for nx=nz= 128 and nx’nz=256, respectively. Breaking off of fluid parcel was not resolved numerically. Chapter 4. Nonlinear Numerical Simulations 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 106 0.2 0.4 0.6 0.8 1.0 0.0 Figure 4.7: Density and vorticity for right mode oniy with 0.2 0.4 = 0.6 0.8 1.0 0 Results of nonlinear simulations when the flow is initialized with the right mode only for J = 0.3, R = 3, Pr = 9, Re = 300, if = 10, and e = 0. Contours are of constant density p = p(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 2ir/&. Chapter 4. Nonlinear Numerical Simulations 0.0 0.2 0.4 0.6 0.6 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 107 0.2 0.4 0.6 0.8 1.0 0.0 Figure 48: Density and vorticity for right mode only with e 0.2 = 0.4 0.25 Results of nonlinear simulations when the flow is initialized with the right mode only for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and a = 0.25. Contours are of constant density p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 2K/&. Chapter 4. Nonlinear Numerical Simulations Figure 4.9: Density and vorticity for right mode only with e 108 = 0.5 Results of nonlinear simulations when the flow is initialized with the right mode oniy for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and = 0.5. Contours are of constant density p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 2ir/. Chapter 4. Nonlinear Numerical Simulations 109 I 0.15 I I I I I + 0.10 U + 0.05 0.00 I 0.0 I I 0.2 I I 0.4 I 0.6 Figure 4.10: Linear growth rates from nonlinear simulations Linear growth rates from: 1. linear stability analysis (A), 2. nonlinear simulations with left or right mode only (+), 3. nonlinear simulations with both modes () when J 0.3, R 3, Pr = 9, Re = 300, and H = 10. Upper values correspond to the right mode and lower values to the left mode. Chapter 4. Nonlinear Numerical Simulations 0 100 200 110 300 400 Figure 4.11: Positions of waves for right mode only Horizontal position of maximum density at critical level for right wave (solid line) and minimum density at critical level for left wave (dashed line) when the flow is initially perturbed with the right mode oniy for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. The horizontal position has been scaled with respect to its period 27r/a*. Chapter 4. Nonlinear Numerical Simulations 111 : 0.0— cycle 1•0: 0.5 c3 = 0 0.0— 0 o o - cycle 1.0 8=0.25 0.5 o’ Et 0.0- 00 + 30 cycle 1.0 0.5 ci- - 0-0 —0.5- 0 0 0 tt.t.t -1.0 0 10 20 30 cycle Figure 4.12: Average phase speeds from nonlinear simulations Average phase speed over a cycle for: 1. nonlinear simulations with right or left mode only (+), 2. nonlinear simulations with both modes (0) when J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. Phase speed predicted by linear theory is indicated by a dashed line. Positive phase speeds correspond to right moving waves and negative wave speeds to left moving waves. Chapter 4. Nonlinear Numerical Simulations 112 2.5 2.0 1.5 1.0 0.5 — 0.0 0 100 200 0 100 200 300 45 0 100 200 300 400 300 400 0.04 000 I —0.02 0.004 0.002 ......... I — — 0 0.000 Figure 4.13: Energy budgets for right mode only with e = 0 Energy budget for nonlinear simulation initially perturbed with the right mode only for J = 0.3, R = 3, Pr = 9, Re = 300, and H 10 with e = 0. Reservoirs: T + 6.1 (solid line), K—2.1 (dashed line), K’ (dotted line), and F+9.1 (dash-dot-dash line). Transfers: C(K, K’) (solid line), and C(F, K’) (dotted line). Diffusions: D(K) (dashed line), D(K’) (dotted line), and D(F) (dash-dot-dash line). Chapter 4. Nonlinear Numerical Simulations 2.5 2.0 ,, 113 - — - 1.5 - 0 70 0.02 140 I 0.00 70 0 0.004 140 — 0.002 — 0 - 0.000— 0 — 140 70 tine Figure 4.14: Energy budgets for right mode oniy with e = 0.5 Energy budget for nonlinear simulation initially perturbed with the right mode only for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10 with e = 0.5. See figure 4.13 for description. Chapter 4. Nonlinear Numerical Simulations U 114 U U 4 —1.0 p —0.5 0.0 p 0.5 1.0 Figure 4.15: Evolution of mean profiles for right mode only Mean velocity and density profiles at 1. t = 0 (dashed line), 2. from numerical simulations initialized with right mode only (solid line), and 3. diffusion of mean flow without perturbations (dotted line). Parameter values are J = 0.3, R = 3, Pr = 9, Re = 300, and H=l0. Chapter 4. Nonlinear Numerical Simulations 115 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 4.16: Density and vorticity for left mode only with e = 0.25 Results of nonlinear simulations when the flow is initialized with the left mode only for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and e = 0.25. Contours are of constant density p Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 27r/&. Chapter 4. Nonlinear Numerical Simulations Figure 4.17: Density and vorticity for left mode only with 116 = 0.5 Results of nonlinear simulations when the flow is initialized with the left mode only for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and e = 0.5. Contours are of constant density p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 2ir/&. Chapter 4. Nonlinear Numerical Simulations 0 100 200 117 300 400 Figure 4.18: Positions of waves for left mode oniy Horizontal position of maximum density at critical level for right wave (solid line) and minimum density at critical level for left wave (dashed line) when the flow is initially perturbed with the left mode only for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. The horizontal position has been scaled with respect to its period 21r/a*. Chapter 4. Nonlinear Numerical Simulations Figure 4.19: Density and vorticity for both modes with 118 E = 0 Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3, R = 3, Pr = 9, Re 300, H 10, and e = 0. Contours are of constant density p p(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 27r/&. Chapter 4. Nonlinear Numerical Simulations 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 119 1. Figure 4.20: Density and vorticity for both modes with e = 0.1 Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and = 0.1. Contours are of constant density p = Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 27r/. Chapter 4. Nonlinear Numerical Simulations 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 120 0.2 0.4 0.6 0.8 1.0 0.0 Figure 4.21: Density and vorticity for both modes with 0.2 0.4 0.6 0.8 1.0 0.25 Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and = 0.25. Contours are of constant density p = Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respect to its period 27r/&. Chapter 4. Nonlinear Numerical Simulations 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 121 1.0 Figure 4.22: Density and vorticity for both modes with e = 0.5 Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and E = 0.5. Contours are of constant density p Pa(Z) + p’ iii the x-z plane and shading is of vorticity. x has been scaled with respect to its period 27r/&. Chapter 4. Nonlinear Numerical Simulations 100 200 122 300 400 Figure 4.23: Positions of waves for both modes Horizontal position of maximum density at critical level for right wave (solid line) and minimum density at critical level for left wave (dashed line) when the flow is initially perturbed with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. The horizontal position has been scaled with respect to its period 27r/a*. Chapter 4. Nonlinear Numerical Simulations 0 100 123 200 300 400 300 400 0.02 000 —0.02 — —0.04 3 I I 100 200 0.004 J/J\/r 0.002 J\/ 1/J\1/\ — 0.000 — 0 — 100 200 400 300 Figure 4.24: Energy budgets for both modes with E = 0 Energy budget for nonlinear simulation initially perturbed with both modes for J R = 3, Pr = 9, Re = 300, and H = 10 with e = 0. See figure 4.13 for description. = 0.3, Chapter 4. Nonlinear Numerical Simulations 124 1.0 0.8 0.6 0.4 0.2 0.0 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 0.002 0 > a) (I) a) 0.000 0.0004 0.0002 a-) —0.0000 C) F— —0.0002 —0.0004 0.004 (I) 0 U) D 0.002 0.000 t Figure 4.25: Wave positions and energy budgets in linear regime for e = 0 Wave positions and energy budgets in linear regime for 0 as described in figure 4.24. Times of maximum wave separation and wave intersection are indicated by vertical lines. Chapter 4. Nonlinear Numerical Simulations 125 1.0 0.8 0.6 x 0.4 0.2 0.0 200 210 220 230 240 250 2.5 2.0 U) 0 > 0 1.5 (I) 0 1.0 - - - - 0.5E o.0_ 200 -. 210 220 230 240 250 230 240 250 230 240 250 0.04 0.02 U) G) U) C 0.00 -. 0 F- —0.02 —0.04 200 210 220 210 220 0.004 U) C 0 U) D 0.002 0.000 200 t Figure 4.26: Wave positions and energy budgets in nonlinear regime for e = 0 Wave positions and energy budgets in nonlinear regime for e = 0 as described in fig ure 4.24. Times of maximum wave separation and wave intersection are indicated by vertical lines. Chapter 4. Nonlinear Numerical Simulations 2.5 126 - 2.0 — 1.5 - 0.0• 0 - - 70 0.04 140 I 0.02 — 0.00 - —0.02 - - - - -. . - - - - - - - - - -.‘ - - - - - - - - - - - - - - - - - - ‘.- - - - - - - - - - - -- - - - -- - S.. -- - - - - — —0.04 I 70 0 0.004 0.002 — 0.000 — 0 — -- 140 70 Figure 4.27: Energy budgets for both modes with e = 0.5 Energy budget for nonlinear simulation initially perturbed with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10 with = 0.5. See figure 4.13 for description. Chapter 4. Nonlinear Numerical Simulations 127 8O.25 U 4 2 0 —2 —1.0 —0.5 0.0 p 0.5 1.0 —1.0 —0.5 0.0 p 0.5 1.0 —1.0 —0.5 0.0 p 0.5 1.0 -1.0-0.5 0.51.0 Figure 4.28: Evolution of mean density profiles for both modes Mean velocity and density profiles at 1. t = 0 (dashed line), 2. from numerical simu lations initialized with both modes (solid line), and 3. diffusion of mean flow without perturbations (dotted line). Parameter values are J = 0.3, R = 3, Pr = 9, Re = 300, and H=10. Chapter 5 Comparison with Lab oratory Experiments Although numerical simulations provide a vast amount of quantitative information about a flow, we are using it to study an isolated phenomenon. In nature, such a phenomenon is unlikely to occur on its own. Also, due to numerical restrictions, it is often difficult, if not impossible, to reproduce the high Reynolds number characteristic of natural systems. These issues lead to the question of whether or not the results of the numerical simulations are indeed physically realistic. There have been field observations of Kelvin-flelmholtz billows, for example, in cloud formations, Colson [11], and in the Mediterranean Sea, Woods [69]. There is, however, a paucity of quantitative data on Holmboe instabilities in the field and in the labora tory. Also, it is difficult to make field observations of Holmboe instabilities due to the non-stationary nature of the instabilities. Thus, although we would like to make de tailed quantitative comparisons between the numerical results and both laboratory and field observations, we restrict our attention to qualitative comparisons with laboratory observations. In order to compare our results with those of laboratory experiments, we draw on three sources: the results of Lawrence et al. [38], Zhu [72], and Guez &? Lawrence [20]. First, we examine the results of Zhu who provides one of the best sequences of photographs of a nearly symmetric Holmboe instability. This instability occurs at a relatively large value of the bulk Richardson number. Next, we compute the evolution of two flows with weak stratification and large density interface displacement and compare these with the results 128 Chapter 5. Comparison with Laboratory Experiments 129 of Lawrence et al. [38] and Guez & Lawrence [20]. Iii all three cases, the stratification is a result of having a layer of fresh water over a layer of salt water. Full descriptions of the experimental setups are given in the individual references. Symmetric Case 5.1 The results from Zhu [72] are shown in figure 5.1. It is evident that there is a wave protruding into the upper layer which is moving to the right and one into the lower layer moving to the left. Since the mean flow is moving to the right, it appears that the right wave is moving faster than the left wave. Flow conditions for this experiment are given in table 5.1. From the data, we calculate the relative speeds of right and left moving waves with respect to the mean flow. The relative phase speeds are 0.6 cm/s and -0.6 cm/s for right and left moving waves, respectively. The fact that the two waves are moving in opposite directions at the same speed relative to the mean flow is indicative that the density interface is at the same level as the center of the shear layer (i.e., = 0). Since the measured waves speeds are approximate, these wave speeds do not contradict the measured range of as it is still possible that the density interface displacement is small and negative. For the purpose of comparison, however, we set e = 0 in the numerical simulation. Due to the large values of both the Reynolds and Prandtl numbers (table 5.1), we cannot match the flow conditions in the numerical simulations since doing so would result in insufficient diffusion to adequately resolve the flow. Also, selecting 10 Re = R 15 with 460 requires a large number of grid points in the z-direction. For the purpose of comparison, we use the same parameters used in chapter 4 to examine the evolutioll of the symmetrie case. The values of these parameters are J R = 3, and = 0.3, Re 300, Pr = 9, 0.54. Parameters relating to the numerical scheme are given in table 4.1. Chapter 5. Comparison with Laboratory Experiments J Re Pr 0.45 0.06 0.03 460 25 25 700 700 700 g’ 2 cm/s 1.6 0.4 0.13 Ui cm/s 2.5 3.7 3.1 U 2 cm/s -1.7 1.7 1.5 h cm 5 0.58 0.6 A cm 14 3.9 4.1 130 R [-0.1,0] 0.5 — cm/s 1.0 0 0 [10,15] 10 10 c$!) cm/s -0.2 * * Table 5.1: Approximate flow conditions of laboratory experiments from Zhu [72] (J = 0.45), Lawrence et al. [38] (J 0.06), and Guez & Lawrence [20] (J = 0.03). A is the observed wavelength of the instability. g’ is the modified gravitational accelera tion. When J = 0.03 and 0.06 there are no observable left moving waves (indicated by * in the table). The lower value of the initial bulk Richardson number used in the numerical simulations is not an unreasonable choice since J increases as the flow evolves. Also, there is some uncertainty in determining the value of J in the laboratory experiments. In order to qualitatively compare our results with those of Zhii [72], we must output our results at the same frequency as the photographs were taken. During the 0.5 second between photographs, a wave has travelled (cr — U) 0.5/A 0.02 of its total wavelength. In non-dimensional time, we must output our results at intervals of t 0.02 . 2n/c 0.389. Here * = 0.02 A*/c = indicates the non-dimensional quantities used in the numerical simulation. Results of the numerical simulation for J = 0.3 with e = 0 are shown in figure 5.2. As in the experimental observations of Zhu [72], there is a wave moving to the right which protrudes into the upper layer and one moving to the left protruding into the lower layer. There are, however, several qualitative differences between the two results. The wave peaks are sharper in the laboratory results of Zhu [72] (reproduced in figure 5.1) than ill the results of the numerical simulations. A possible explanation for this difference is that the density interface is much thinner for the experimental flow than in the numerical Chapter 5. Comparison with Laboratory Experiments 131 simulation. Also, it is difficult to determine at what stage of the perturbation evolution the flow shown in the photographs occurs. As seen in figure 4.19, the shape of the waves changes throughout the evolution of the perturbation. Another observable difference is that in the experiments, the amplitude of the left mode appears larger than the amplitude of the right mode, whereas in the numerical simulations they have the same amplitude. There are two possible explanations for this difference. First, as shown in figure 4.20, the two waves can have significantly different amplitudes even when the density interface displacement is small. Thus it may be possible that the density interface shift is a small negative number. This possibility is supported by the range of e given in table 5.1 for the resilits of Zhu [72]. Second, in chapter 4 it was discovered that the relative amplitudes of right and left modes during early stages of their nonlinear evolution depends on the relative initial amplitude of the two modes. Since the Holmboe instability observed by Zhu [72] ocdllrred in a facility designed to study an exchange flow through a channel with an underwater sill, the Holmboe instability is just one of several mechanisms that are taking place (see figure 7.la of Zhii [72]). It is possible that one of the other mechanisms has an influence on the amplitude of the two modes. In spite of these differences, we feel that the comparison between the numerical and experimental results are favourable. 5.2 Non-Symmetric Case Experimental results from Guez & Lawrence [20] and Lawrence et al. [38] are shown in figures 5.3 and 5.4 respectively. Approximate flow conditions are given in table 5.1. These laboratory observations are clear examples of flows exhibiting “one-sidedness”. In both cases, perturbations take the form of fluid being dragged from the lower layer into the upper layer. This behaviour is characteristic of the non-symmetric case with > 0. Chapter 5. Comparison with Laboratory Experiments 132 Since, for both cases, the bulk Richardson number is small, billowing is observed and increases as J decreases. Although this set of experiments occurs at a lower Reynolds number, the Prandtl number and the ratio R are still too large to numerically model. The lower value of the Reynolds number does, however, allow us to easily compute the evolution of the flow for several values of R (and hence several values of Pr) to see if any trends are present. For this purpose, we use the flow parameters from the experimental results of Guez & Lawrence [20]. Since the value of e is not given for this experiment, we set e = 0.5. This value is given by Lawrence et al. [38]. Since experiments of Guez & Lawrence [20] were performed in the same facility as those of Lawrence et al. [38], it is not unreasonable to assume that the density interface displacements are approximately equal. When J = 0.03 and e = 0.5, the right moving instability is the only unstable mode and is the only mode contained in the initial perturbation. At low Reynolds numbers, the initial amplitude of the perturbation has a significant effect on how the flow develops. If the initial amplitude is too small, then diffusion acts to stabilize the flow before the perturbation has a chance to fully develop. This phenomenon was observed by Klaassen & Peltier [29] in their study of Kelvin-Helmholtz instabilities. In the experiments of Guez & Lawrence [20], the flow is forced at a frequency corresponding to a = 0.45, allowing us to initialize the perturbations with a relatively large amount of kinetic energy (K’ = 0.005). To compare our results with the photographs in figure 5.3, we wish to output the results at intervals equal to the time required for the wave to travel one wavelength. This time interval, t0 is given by = — A U+Cr 2irh/2 (U+c.6U) 1 where c is the non-dimensional phase speed obtained from linear stability analysis. In Chapter 5. Comparison with Laboratory Experiments 133 non-dimensional time, we obtain — — = irh &(U + c 6U) 27r6U &(U + c 6U) 4.75 Results from the numerical simulations with R 5’U h/2 (5.2) = 3, 6, and 10 (and hence Pr 36, and 100, respectively) are shown in figures 5.5, 5.6 and 5.7. In all cases, nx nz = 128, and At = = = 9, 64, 0.0625. Although the results are shown at intervals of 4.75, as calculated above, the billows nevertheless move across the frame as time increases. This movement is explained by recalling that, as discussed in Chapter 4, the phase speed of a right moving wave increases with time when e 0. For all values of R, the flow behaves qualitatively the same. A finger of heavier fluid is drawn up into the upper layer and forms a billow. The finger becomes thinner with a sharper tip as R (and hence Pr) increases. For R separated by t = = 6 and 10, the billow has 57. As the billow wraps around itself, it is stretched and becomes thinner. Eventually, the upper portion of the billow becomes so thin that it diffuses into the ambient fluid. At this point, the billow divides into two disconnected parts. When R = 3, the billow remains connected at all times since the initial thickness of the billow is much larger than that of the other two cases. It is observed that the diffusion rate increases as R decreases. This change in diffusion rate can be related to the decreasing value of the Prandtl number and is consistent with the results of Klaassen & Peltier [28] who studied the effect of the Prandtl number on Kelvin-Helmholtz billows. For R = we observe a resolution problem at the selected grid size. Since the results for R 6 are qualitatively the same as for the R = = 10 10 case and do not change when the grid is refined, we use this case for comparison with the experimental results of Guez & Lawrence [20]. Although the results from the numerical simulations depicted in figure 5.6 and those Chapter 5. Comparison with Laboratory Experiments 134 of the laboratory experiments of Guez & Lawrence [20] reproduced in figure 5.3 are qualitatively the same, there are several differences between them. In the numerical simulations the tips of the figures in the billows are not as sharp as those of the laboratory experiments. This difference between experimental and numerical observations is a result of the values of R and Pr used in the numerical simulations being much smaller than those of the experiments. This difference becomes less pronounced as R and Pr are increased. Another noticeable difference between the two results, is that there is not as much wrapping around of the billows in the numerical simulations as in the laboratory experi ments. There are two possible causes of this discrepancy. Firstly, Klaassen & Peltier [28] found that as the Prandtl number increases, there is more wrapping around of fingers in a Kelvin-Helmholtz instability before the displaced fluid diffuses into the ambient flow. Since we have used Pr = 36, which is much smaller than the actual value, to model the flow, perhaps the large amount of diffusion suppresses the wrapping around of the billow. Secondly, there is some uncertainty in determining the flow conditions for the laboratory experiments. It is possible that the actual value of the bulk Richardson number is slightly smaller than 0.03. Since the amount of billowing is sensitive the the value of the bulk Richardson number, we expect more curling up of the billows as J decreases. Breakdown of the billows occurs at approximately the same time in both the numerical simulations and the laboratory experiments. The cause of this breakdown is different, however, for the two cases. In the laboratory experiment the billows collapse when the flow becomes three-dimensional. Since the numerical simulations restrict the flow to two-dimensions, transition to three-dimensional flow cannot possibly be the cause of the breakdown of the billows. Instead, it is due to the top part of the billow diffusing into the ambient fluid. Because of this difference, it does not make physical sense to compare the two flows beyond this point. Nevertheless, it is of interest to note that the two flows remain qualitatively the same. Chapter 5. Comparison with Laboratory Experiments 135 Next we examine the experimental results of Lawrence et al. [38] when J ure 5.4). Since the results of the nonlinear simulations with R = = 0.06 (fig 6 compared favourably with the experimental results of Guez & Lawrence [20] and both experiments were con ducted in the same facility, we use the same value to compare with the results of Lawrence et al. [38]. The only difference is that no forcing is present in the laboratory study of Lawrence et al. [38]. Nevertheless, we found that a fairly large initial amplitude of the perturbation was still required in order for the perturbation to become fully developed. Having a large initial amplitude poses no problems since only the right moving mode is an unstable mode when J = 0.06 with = 0.5. When both right and left modes are un stable, more care must be taken in choosing the initial amplitudes of the perturbations, as is discussed below. Results of the numerical simulations with J = 4.95). As with the J = = 0.06 are shown in figure 5.8 (here 0.03 case, the fingers of the billows are not as thin as in the experimental results of Lawrence et al. [38] (figure 5.4). Also, the fingers do not curl up as much in the numerical simulations. The possible explanations for these differences are the same as the J = 0.03 case. Although the results of the numerical simulations and the laboratory experiments differ at later times, qualitative agreement as the perturbation initially develops is excellent. Lawrence et al. [38] examined the development of perturbations for a range of Richard son numbers (figure 5.4). Comparing numerical results with those of Lawrence et al. [38] is difficult when J is large and both right and left moving modes are unstable. The prob lem lies in selecting the correct initial amplitude of the perturbation. If the amplitude is not sufficiently large, then diffusion damps the perturbation before it has a chance to grow. If the initial amplitude is too large, the perturbation does not spend much time in the linear regime, resulting in the left mode being of approximately the same amplitude Chapter 5. Comparison with Laboratory Experiments 136 as the right mode. The experimental observations indicate, however, that the initial am plitude is such that the perturbation spends sufficient time in the linear regime for the right mode to dominate. Finding the correct initial amplitude is a matter of trial and error. We have not attempted to numerically reproduce these experiments. Chapter 5. Comparison with Laboratory Experiments 137 Figure 5.1: Experiments of Zhu [72] Sequence of photographs showing Holmboe waves (figure 7.lb of Zhu [72]). The upper layer is from left to right, the lower layer from right to left. The grid markings are 5 cm apart and the photographs were taken at 0.5 second intervals. See table 5.1 for flow conditions. — . C C C I! c_ CD ) II CjD CD Cl) C Ci) CD C Ci) -4 CD CD -4. b p p 0 p I. I I Chapter 5. Comparison with Laboratory Experiments 139 Figure 5.3: Experiments of Guez & Lawrence [20] Sequence of photographs showing development of Kelvin-Helmholtz type billows from Guez & Lawrence [20]. The flow is from left to right. The field view of each photograph is approximately 25cm with an overlap of approximately one wavelength between each successive photograph. See table 5.1 for flow conditions. Chapter 5. Comparison with Laboratory Experiments 140 I I In II II II II II Figure 5.4: Experiments of Lawrence et al. [38] Sequence of photographs showing development of shear instabilities at various Richardson numbers (figure 6 from Lawrence et al. [38]). The flow is from left to right. The field view of each photograph is approximately 25cm with an overlap of approximately 5cm between each successive photograph. As the bulk Richardson number increases, less billowing is observed. See table 5.1 for flow conditions when J = 0.06. II II CD C :J) C CD CD o C II CD - a-. C Ci) C. CD CD CD CD —1 o c± C)) 0 -> CD cn. C CD cJ 0 CD U) CD CD C Ci) Ci) CD C)1 NJ P 0 p b p p 0 0 p p p p NJ 0 p NJ NJ 0 0 NJ • 0 0 0 C N) 9) C C a p 0 0 p NJ 0 p NJ 0 0 0 C, N) C, 0 0 0 S p o S K....,. 1 S NJ 0 NJ .i[ C C N) P S 0 P 0 0 P P NJ p o p p NJ p 0 p p p 0 • I • NJ I NJ NJ 0 0 0 NJ NJ NJ • • C C C b ID 0 N) CI, 0 C), 0 0 C 0 0 C) Chapter 5. Comparison with Laboratory Experiments 9cqo•o t= 14.2500 t= 4 = 19.0000 .7500 . = 4 oo2o:4o:6 142 0:8 4 4 1.0 28.5000 t 33.2500 t 4 4 = 38.0000 4 t = 42.7500 4 : . = 4 04 .•. 08. •1.O 47.5000 t = 52.2500 t = 57.0000 t = 61.7500 66.5000 t = 71.2500 t = 76.0000 t = 80.7500 4 4 4 0.0.26081.0 o.o:20o.8l.o o.oo:24’o:6o81.o Figure 5.6: Results from numerical simulations for J = 0.03 with R = 6 Density contours of results from numerical simulations with J = 0.03, Re = 25, R = 6, Pr = 36, = 0.5, with & 0.45. Contours are from -0.72 to 0.72 in equal increments of 0.18. x has been scaled with respect to its period 27r/&. Chapter 5. t Comparison with Laboratory Experiments = 950000 t 4 4 2 2 = 14.2500 I = 143 190000 4 = 23.7500 = 42.7500 4 0 —2 —4 0 2 —4 0.2 0.4 = 0.6 0.8 1.0 00 0.2 28.5000 0.4 = 4 0.6 0.8 1.0 33.2500 I = 38.0000 4 4 2 0 —2 —4 0.0 —4 0.2 0.4 0.6 0.8 1.0 0.0 —4 0.2 0.4 0.6 2.8 1.0 0.0 0.2 =__47.5000 0.4 = 0.6 0.8 1.0 61.7500 57.0000 4 4 4 2 0 —2 —4 0.2 0.4 = 4 0.6 0.8 1.0 0.2 66.5000 = 0.4 0.6 0.8 1.0 C 3 0.2 71.2500 0.4 t 4 = 0.6 0.8 1.0 0.8 1.0 80.7500 4 2 2 0 0 —2 —4 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 Figure 5.7: Results from numerical simulations for J 1.0 = 0.0 0.2 0.4 0.03 with R 0.6 = 10 Density contours of results from numerical simulations with J = 0.03, Re = 25, R = 10, Pr = 100, e = 0.5, with & = 0.45. Contours are from -0.72 to 0.72 in equal increments of 0.18. x has been scaled with respect to its period 27r/&. Chapter 5. Comparison with Laboratory Experiments = 0.00000 t = 4.93750 t = 144 9.87500 t 4 4 4 4 2 2 2 2 —4 0.00.2 —4 0.4 0.60.8 1.0 0.0 —4 0.2 0.60.81.0 19.7500 4 0.4 0.6 1.0 0.00.20.40.60.81.0 29.6250 t= 34.5625 t = 54.3125 t = 74.0625 4 0!8 t = 44.4375 t = 4 o!4 t 59.2500 4 49.3750 4 0:6 = 4 0.2 4 4 O.OO:2 z 1.0 39.5000 4 0.0 14.8125 —4 87 4 0:6 ---- — - = 64.1875 t = 0:6 69.1250 4 4 1.0 Figure 5.8: Results from numerical simulations for J 0.06 Density contours of results from numerical simulations with J = 0.06, Re = 25, R = 6, Pr = 36, e = 0.5, with & = 0.45. Contours are from -0.72 to 0.72 in equal increments of 0.18. x has been scaled with respect to its period 27r/&. Chapter 6 Summary and Conclusions We have studied the effect of displacing a thin density interface with respect to the center of the shear layer on the stability of a parallel flow. Linear stability analysis was used to determine how increasing the density interface displacement changed the overall stability of the flow. In particular, we concentrated on how this displacement affects Holmboe instabilities. Results of both inviscid theory and viscous theory were considered. Numerical simulations were then performed to examine the nonlinear evolution of these instabilities. Finally, results of the nonlinear simulations were compared with those of laboratory experiments. We only considered displacements of the density interface below the center of the shear layer. 6.1 Linear Stability Analysis When the density interface is displaced with respect to the center of the shear layer, results of linear stability analysis indicate that right and left moving modes, which have equal growth rates when the density interface displacement is zero, split into a faster growrng mode in one direction and a more slowly growing mode moving in the other direction. When the effects of horizontal boundaries are negligible and the density inter face is placed below the center of the shear layer, it is the right moving wave that has the larger growth rate. 145 Chapter 6. Summary and Conclusions 6.1.1 146 Inviscid Results Inviscid theory was used to study the linear stability of piecewise linear background profiles. This has the advantage of admitting an exact dispersion relation. Also, other than the density interface displacement, the only additional parameter of this flow is the bulk Richardson number. For the non-symmetric case, decreasing the height of the vertical domain damps the growth rate of the right mode more quickly than the growth rate of the left mode. This difference in rate of decrease for the two modes results in the left mode having the larger growth rate when the total height of the domain is less than or equal to five times the shear layer thickness. When the density interface corresponds to the center of the shear layer (the symmetric case), the results of Smyth [51] were reproduced with transition from Kelvin-Helmholtz to Holmboe instabilities occurring at smaller values of the bulk Richardson number as the domain height decreases. For pure flolmboe waves, there is a region where the maximum growth rate increases with increasing Richardson number, indicating that three-dimensional instabilities are possibly more unstable than two-dimensional instabilities. It was found, however, that it is indeed the two-dimensional instabilities that dominate. For non-symmetric Holmboe instabilities, the maximum growth rate of the right mode decreases monotonically as the Richardson number increases and so there is no possibility of the three-dimensional instabilities dominating. The left mode still has a region for which the growth rate increases with the Richardson number, but, as with the symmetric case, two-dimensional instabilities remain more unstable. 6.1.2 Viscous Results Since the nonlinear numerical simulation requires that there are no discontinuities in the initial data or its derivatives, results of linear stability analysis for a viscous flow with Chapter 6. Summary and Conclusions 147 hyperbolic tangent background profiles were examined. Three additional parameters were introduced: the Reynolds number, the Prandtl number, and the ratio of shear layer thickness to density interface thickness, R. As the density interface displacement was increased, the precise behaviour of right and left modes depended on the value of R. For large finite values of R (R > 4), the behaviour was similar to the discontinuous density interface case which corresponds to R = cc. The most unstable mode of the right moving wave occurs at a smaller wave number whereas the maximum growth rate of the left mode occurs at a larger wave number. When R < 4, however, there is a change in this behaviour. The most unstable right wave now occurs at a larger value of the wave number. The limb on the stability diagram corresponding to the left moving wave continues to tilt towards the wave number axis, but also separates itself from the axis. This behaviour becomes more pronounced as R decreases. As with the inviscid case, the effect of decreasing the domain height was examined. For the symmetric case, the behaviour was significantly different than the inviscid case. Instead of the transition from Kelvin-Helmholtz to Holmboe instabilities being shifted to a smaller Richardson number as the domain height decreases, the unstable modes near this transition are the first to be stabilized. This is a direct result of the modes in this region being the least unstable. When the total domain height is less than or equal to five times the height of the shear layer thickness, there is a stable region separating Kelvin-Helmholtz instabilities from Holmboe instabilities. For the non-symmetric case, the growth rate of the right mode decreases more rapidly than that of the left mode as the boundaries get closer to the shear layer. Unlike the inviscid case, however, the right mode always remains more unstable than the left mode. For R = 3 with Pr = 9 and Re = 300, Smyth & Peltier [55] have shown that three dimensional instabilities are more unstable than two-dimensional instabilities when the density interface displacement is zero. As the density interface displacement increases, Chapter 6. Summary and Conclusions 148 two-dimensional instabilities dominate for the faster growing right mode. The left mode, however, still exhibits characteristics necessary for three-dimensional instabilities to be dominant and thus may be a mechanism through which perturbations become threedimensional. Further investigation is required to determine if three-dimensional instabil ities are indeed dominant for the left moving wave. As R increases, three-dimensional instabilities are less likely to be primary instabilities as the results tend to those of the piecewise linear case. Studying the shapes of the eigenfunctions predicted from linear stability analysis pro vided insight into the growth mechanism of the instabilities. For the symmetric case, Smyth & Peltier [54] showed that when both right and left modes are present, the insta bilities oscillate between growth and decay phases of equal duration. During the growth phase, the critical level acts as a source, extracting energy from the mean flow and the density interface as a sink, converting perturbation kinetic energy to potential energy. The same behaviour is observed for the non-symmetric case. One important difference, however, is that, as the density inlerface displacement increases, the right mode is re sponsible for a larger portion of the energy exchange between the mean kinetic energy, the perturbation kinetic energy, and the potential energy reservoirs. When the density interface displacement is small (0 0.25), the growth and decay of the perturbations are in phase with and governed by the transfer of energy between the mean flow and the perturbations. The maximum in the perturbation kinetic energy occurs just after the waves have passed through each other. The minimum occurs just after the waves have reached the stage of maximal separation and are starting to approach each other. When the density interface displacement is large (E = 0.5), the growth and decay of the perturbations are still governed by the exchange of energy between the mean flow and the perturbations, but the two quantities are no longer in phase. Chapter 6. Summary and Conc1u.ions 6.2 149 Nonlinear Simulations Numerical simulations were used to study the nonlinear evolution of the perturbations. In order to examine the effect of initial conditions on the development of the perturbations and to study nonlinear interactions between the right and left moving waves, we ran three sets of simulations. We first perturbed the flow with the more unstable right mode. Next, we tracked the evolution of a flow initially perturbed with the more slowly growing left mode. Finally, we perturbed the flow with both right and left modes. For the symmetric case, the same behaviour is observed when the flow is initially perturbed with either the right or left mode. The mode that is initially present dominates the flow at early stages of the evolution. Eventually, the other mode grows and the flow tends to a state where both modes are present. When the perturbation initially contains both modes, they grow at the same rate and the flow remains symmetric throughout its evolution. Thus, although the early behaviour depends on which modes are used to initialize the flow, long term behaviour appears to be independent of initial conditions. When the flow is initially perturbed with the right mode only, the amount of time spent in the linear regime decreases as the density interface displacement increases. This trend is due to the increased growth rate of the right mode, as predicted by linear stability analysis. The left mode does eventually develop, but takes longer to become noticeable as the density offset is increased. Once again, the behaviour can be attributed to the results of linear stability analysis which predict a smaller growth rate for the left mode as the density offset is increased. When the perturbations initially contain only the left mode, the flow spends more time in the linear regime as the density offset increases. Also, the right mode affects the flow more quickly. For very large density offsets, the right mode affects the flow while it is still in the linear regime and dominates the flow at an early stage in its evolution. Chapter 6. Summary and Conclusions 150 The behaviour of the flow when both modes are initially present is not surprising in view of the above results. When the density offset is large, the right mode quickly dominates the flow. The left mode, however, continues its slow growth and becomes noticeable at later times in the simulation. One problem encountered was our inability to resolve the separation of a fluid parcel from the perturbation. Resolution problems occurred when the density offset was large and restricted the time for which the simulation was valid. Thus little can be said about the long term behaviour of such flows. Growth rates in the linear stage of the simulations were compared to the results of linear theory. Excellent agreement was found for the symmetric case. For the non- symmetric case, the right mode grows more slowly and the left more quickly than the predicted values. When the flow has both modes present, the overall linear growth rate lies between those of the right only and left only cases, but is dominated by the right mode. Initial phase speeds agree with those predicted from linear stability analysis. As time increases, the phase speed of the right moving wave increases. That of the left moving wave increases for small density offsets and decreases for large density offsets. This behaviour indicates that perhaps all waves are tending to the same phase speed in the nonlinear regime. When oniy one mode is present, there are no oscillations in the phase speed. When both modes are present, however, the phase speeds of the waves change as they pass through each other. In the nonlinear regime, the waves speed up as as they approach each other and slow down as they pass each other. Also, there is a decrease in the average phase speed of the waves when both waves are present. In contrast to the linear regime where oscillations in the perturbation kinetic energy are governed by the exchange between the mean flow and the perturbations, in the nonlinear regime, oscillations in the perturbation kinetic energy are closely linked to the exchange of energy between the potential energy reservoir and the perturbation kinetic Chapter 6. Summary and Conclusions 151 energy reservoir. Maxima in the perturbation kinetic energy occur when the waves are passing through each other, whereas the minima occur when the waves are maximally separated. 6.3 Comparison with Laboratory Experiments Results of the numerical simulations were compared with the experimental observations 0.4, e of Zhu [72] (J al. [38] (J 0.06, e 0), Guez & Lawrence [20] (J 0.03, e 0.5), and Lawrence et 0.5). Due to numerical considerations, we were unable to match the high Prandtl number and thin density interface thickness of the salt stratified laboratory experiments. To determine the effect of the Prandtl number on the development of instabilities, we ran a series of numerical simulations with several values of the Prandtl number. The density interface thickness used in the simulations are related the Prandtl number by R = /N, where 1/R is a measure of the density interface thickness. As the Prandtl number increased and the density interface thickness decreased, the waves or billows formed by the instabilities became more sharply defined. Good qualitative agreement was obtained between the results of the numerical simu lations and those of the laboratory experiments. Main differences were due to the thicker density interface and the lower Prandtl number used in the numerical simulations. Since obtaining quantitative information about a flow from experimental observations is often difficult, simulations may be of use to quantify the instabilities studied here. 6.4 Future Work Although the above study has allowed us to gain better insight into the nonlinear devel opment of non-symmetric Holmboe instabilities, many issues have been brought up and remain unresolved. Some of the more important ones are mentioned here. Chapter 6. Summary and Conclusions 152 Due to numerical restrictions, we were unable to resolve the long term behaviour of the flow. This behaviour is of interest as it appears that the eventual state of the flow is independent of the initial conditions. Also, computing the long term evolution of the flow would allow one to verify the conjecture that the phase speed of the waves asymptotes to a single value. When the density offset was large, denser fluid from the lower layer separated from the wave and was released into the less dense upper layer. Since we were unable to consistently compute the evolution of this phenomenon, it is desirable to find out the time and manner in which it occurs. Resolving this mechanism would also permit the study of long term behaviour of these instabilities. The maximum amplitude of the perturbation kinetic energy appears to be indepen dent of the initial conditions. Since the perturbation kinetic energy contains the total contribution of both modes, determining the individual contributions of right and left modes would allow one to better study the interaction between the two waves. The nonlinear simulations examined here forced right and left modes to have the same wave length. For the non-symmetric case, linear stability analysis indicates that the most unstable right and left modes occur at different wave numbers. It would be of interest to see how this would affect the evolution of the instabilities. Finally, we restricted our study to the two-dimensional development of instabilities. Although this is often valid for the initial development of the perturbations, it is well documented that the flow eventually becomes three-dimensional. Hence determining the manner in which it becomes three-dimensional is of importance. List of Symbols (j horizontally averaged (mean) qnantity, unless otherwise noted 74 (.) vertically integrated quantity 74 (“) discrete Fourier transform of a quantity 76 () differentiation with respect to time 76 a nondimensional wave number 13 real part ofa 14 wave number used in numerical simulation 74 wave number at which the maximum growth rate occurs 30 & 7 coefficient of expansion for stratifying agent 8 two-dimensional Laplacian operator 8 tIh discrete two-dimensional Laplacian operator Lip maximum density variation Lit time step size in numerical simulation 77 size of horizontal space step in numerical simulation 75 Liz size of vertical space step in numerical simulation 75 6 delta function 25 density scale used for non-dimensionalizing 12 velocity scale used for non-dimensionalizing 12 amplitude of perturbations 23 e nondimensional vertical density interface displacement 25 ‘7 density interface thickness 31 8 concentration of stratifying agent SU 75 9 153 8 List of Symbols 154 concentration of stratifying agent at standard state angle with respect to x-axis of three-dimensional perturbations 8 . . .29 . it diffusivity of stratifying agent A dimensional wave length A amplification factor of semi-discretized scheme p molecular viscosity 8 ii kinematic viscosity 11 p density of fluid 8 130 77 8 complex vertical amplitude of p’ from normal mode analysis 13 p’ deviation of density from hydrostatic equilibrium 9 Pa vertical density distribution of fluid in hydrostatic equilibrium 9 Fourier coefficient for 3 32 Po density at the standard state 8 P1 mean density in upper layer 3 P2 mean density in lower layer 3 approximate value of p’(iAx, j/Xz, t) from numerical simulation complex vertical amplitude of x i/” from normal mode analysis w 76 13 Fourier coefficient for 32 nth stream function amplitude function for wave number a 15 phase angle of perturbation kinetic energy 47 initial mean stream function distribution 13 approximate value of “(iAx, jz, t) from numerical simulation 1’ ... .. . 76 stream function 11 perturbation of stream function from steady state 13 vorticity 79 List of Symbols complex vertical amplitude of vorticity 155 ( — iaii) 45 A matrix obtained from Galerkin’s method to solve linear problem A coefficients to complete solution to linear problem a th coefficient of dispersion relation for piecewise linear profiles ,R 1 G(R ) 2 conversion of energy from reservoir R 1 to reservoir R 2 81 c complex nondimensional wave speed 13 th wave speed associated with the wave number a 15 Cg group velocity 14 c, imaginary part of c 13 c. real part of c 13 complex wave speed associate with right moving mode 28 complex wave speed associate with left moving mode 28 loss of energy from reservoir R 1 due to diffusion 81 discrete first derivative in x-direction 75 discrete second derivative in x-direction 75 discrete first derivative in z-direction 75 discrete second derivative in z-direction 75 d vertical density interface displacement 25 g gravitational acceleration g’ modified gravitational acceleration h shear layer thickness 25 H total depth of flow 25 imaginary part of a complex value 15 ) 1 D(R J bulk Richardson number J value of J below which Holmboe instabilities do not occur .. 33 15 .. .26 8 130 3 39 List of Symbols 156 JT transition from Kelvin-Helmholtz to Holmboe waves occurs at JT K kinetic energy of mean flow 80 K’ perturbation (or wave) kinetic energy 74 K’ filtered perturbation kinetic energy 80 L length scale used for non-dimensionaliziug 12 N number of Fourier modes used to solve linear problem 32 nx number of grid points in x-direction for numerical simulation 83 nz number of grid points in z-direction for numerical simulation 83 F potential energy 80 Pr Prandtl number 12 p pressure 8 deviation of pressure from hydrostatic equilibrium 9 27 . complex vertical amplitude of p’ from normal mode analysis 29 Pa vertical pressure distribution of fluid in hydrostatic equilibrium R ratio of shear layer thickness to density interface thickness 31 real part of a complex value 13 Re Reynolds number 12 Sc Schmidt number 12 Ri gradient Richardson number 19 T total energy 81 t time 9 .... 8 time interval between outputs for nonlinear numerical simulation 130 . U initial mean horizontal velocity distribution 24 (I average horizontal velocity between upper and lower layers 12 1 U meau velocity in upper layer 3 List of Symbols 157 2 U mean velocity in lower layer u velocity vector 29 ü complex vertical amplitude of u’ from normal mode analysis 29 u’ perturbation velocity vector 29 u horizontal component of fluid velocity ü complex vertical amplitude of u’ from normal mode analysis 45 streamwise component of velocity perturbation 23 3 8 instantaneous perturbation of horizontal velocity from mean flow .74 V w vector containing Fourier coefficients q. and pn 33 spanwise component of velocity perturbation 29 vertical component of fluid velocity 8 complex vertical amplitude of w’ from normal mode analysis 29 vertical component of velocity perturbation 23 x horizontal coordinate in streamwise direction y horizontal coordinate in spanwise direction z vertical coordinate 8 29 8 Bibliography [1] P. BALDwIN AND P. H. ROBERTS, The critical layer in stratified shear flow, Math ematika, 17 (1970), pp. 102—119. [2] W. E. BOYCE AND R. C. DIPRIMA, Elementary Differential Equations and Bound ary Value Problems, John Wiley & Sons, New York, 3rd ed., 1977. [3] P. BRADSHAw, An Introduction to Turbulence and its Measurement, Pergamon Press, Oxford, 1971. [4] F. K. BR0WAND AND Y. H. WANG, An experiment on the growth of small distur bances at the interface between two streams of different densities and velocities, in Proceedings of the International Symposium on Stratified Flows, Novosibirsk, Soviet Union, August 29 311972. - [5] F. K. BR0wAND AND C. D. WINANT, Laboratory observations of shear-layer in stability in a stratified fluid, Boundary-Layer Meteorology, 5 (1973), pp. 67—77. [6] C. CANUTO, M. Y. HuS5AINI, A. QuATER0NI, AND T. A. ods in Fluid Mechanics, Springer-Verlag, Berlin, 1988. ZANG, Spectral Meth [7] C. P. CAuLFIELD, Stratification and Buoyancy in Geophysical Flows, PhD thesis, Cambridge University, Cambridge, England, 1991. [8] C. P. CAuLFIELD AND W. R. PELTIER, Three dimensionalization of the stratified mixing layer, Physics of Fluids, 6 (1994), pp. 3803—3805. [9] S. CHANDRASEKHAR, Hydrodynamic Press, Oxford, 1961. [10] B. W. CHAR, AND S. M. York, 1991. [11] D. CoLsoN, B. Hydromagnetic Stability, Oxford University K. 0. GEDDEs, M. B. MONAGAN, G. H. GONNET, Maple V Language Reference Manual, Springer-Verlag, New LE0NG, WATT, Wave-cloud formation at Denver, Weatherwise, 7 (1954), pp. 34—35. [12] P. G. DRAzIN AND W. H. Press, Cambridge, 1981. REID, Hydrodynamic Stability, Cambridge University 158 Bibliography 159 [13] C. B. EMMANUEL, B. R. BEAN, L. G. MCALLISTER, AND J. R. POLLARD, Observations of Helmholtz waves in the lower atmosphere with acoustic sounder, Journal of the Atmospheric Sciences, 29 (1972), pp. 886—892. [14] H. J. S. FERNANDO, Turbulent mixing in stratified fluids, Annual Review of Fluid Mechanics, 23 (1991), pp. 455—493. E. J. LIsT, R. C. Y. KoR, J. IMBERGER, AND N. H. BROOKS, Mixing in Inland and Coastal Waters, Academic Press Inc., San Diego, 1979. [15] H. G. FIsHER, [16] K. S. GAGE, Linear viscous stability theory for stably stratified shear flow: A review, Boundary-Layer Meteorology, 5 (1973), pp. 3—17. [17] M. GAsTER, A note on the relation between temporally-increasing and spatiallyincreasing disturbances in hydrodynamic stability, Journal of Fluid Mechanics, 14 (1962), pp. 222—224. [18] E. E. GOSSARD AND J. H. RICHTER, The shape of internal waves of finite am plitude from high-resolution radar sounding of the lower atmosphere, Journal of the Atmospheric Sciences, 27 (1970), pp. 971—973. [19] M. C. GREGG, Diapycnal mixing in the thermocline: A review, Journal of Geo physical Research, 92 (1987), pp. 5249—5286. [20] L. G. GuEz AND G. A. LAwRENCE, Hydrogen bubble tracking and the determi nation of velocity and vorticity fields in a stratified mixing layer. Submitted to Experiments in Fluids. [21] P. HAzEL, Numerical studies of the stability of inviscid stratified shear flows, Journal of Fluid Mechanics, 51(1972), pp. 39—61. [22] H. VON HELMHOLTz, On discontinuous movements of fluids, The London, Ed inburgh, and Dublin Philosophical Magazine and Journal of Science, 36 (1868), pp. 337—346. [23] J. H0LMB0E, On the behavior of symmetric waves in stratified shear layers, Geofys. Publ., 24 (1962), pp. 67—113. [24] L. N. HOWARD, Note on a paper of John W. Miles, Journal of Fluid Mechanics, 10 (1961), pp. 509—512. [25] W. KELVIN, Hydrokinetic solutions and observations, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42 (1871), pp. 362—377. Bibliography 160 [26] G. H. KEuLEGAN, Interfacial instability and mixing in stratified flows, Journal of Research of the National Bureau of Standards, 43 (1949), pp. 487-500. E27] G. P. KLAAssEN, The Transition to Turbulence in Stably Stratified Parallel Flows, PhD thesis, University of Toronto, Toronto, Canada, 1982. [28] G. P. KLAAssEN AND W. R. PELTIER, The effect of Prandtl number on the evolu tion and stability of Kelvin-Helmholtz billows, Geophysical and Astrophysical Fluid Dynamics, 32 (1985), pp. 23—60. [29] , Evolution of finite amplitude Kelvin-Helmholtz billows in two spatial dimen sions, Journal of Atmospheric Sciences, 42 (1985), pp. 1321—1339. [30] The onset of turbulence in finite-amplitude Kelvin-Helmholtz billows, Journal of Fluid Mechanics, 155 (1985), pp. 1—35. [31] The role of transverse secondary instabilities in the evolution of free shear layers, Journal of Fluid Mechanics, 202 (1989), pp. 367—402. , , [32] C. G. Koop, Instability and turbulence in a stratified shear layer, Tech. Report Rep. USCAE 134, Univ. Southern Calif., Dept. Aerospace Ellg., 1976. [33] C. G. KooP AND P. K. BROWAND, Instability and turbulence in a stratified layer with shear, Journal of Fluid Mechanics, 93 (1979), pp. 135—159. [34] C. G. Koop, C. D. cation. WINANT, AND F. K. BR0wAND, 1995. Personal communi [35] D. K0PPEL, On the stability of flow of a thermally stratified fluid under the action of gravity, Journal of Mathematical Physics, 5 (1964), pp. 963—982. [36] H. LAMB, [37] G. A. Hydrodynamics, Cambridge University Press, Cambridge, 6th ed., 1932. LAwRENcE, 1995. Personal communication. F. K. BR0WAND, AND L. G. REDEK0PP, The stability of a sheared density interface, Physics of Fluids A, 3 (1991), pp. 2360—2370. [38] G. A. LAWRENCE, [39] S. A. MAsL0wE, Shear flow instabilities and transition, in Hydrodynamic Insta bilities and the Transition to Turbulence, H. L. Swinney and J. P. Gollub, eds., Springer-Verlag, Berlin, 2nd ed., 1985, pp. 181—228. [40] T. MAxw0RTHY AND F. K. BR0WAND, Experiments in rotating and stratified flows: oceanographic application, Annual Review of Fluid Mechanics, 7 (1975), pp. 273—305. Bibliography 161 [41] J. W. MILES, On the stability of heterogeneous shear layers, Journal of Fluid Me chanics, 10 (1961), PP. 496—508. [42] A. S. M0NIN AND A. M. YAGL0M, Statistical Fluid Mechanics, The MIT Press, Cambridge, Massachusetts, 1971. [43] S. NIsHIDA AND S. YOSHIDA, Stability and eigenfunctions of disturbances in strat ified two-layer shear flows, in Proceedings of the Third International Symposium on Stratified Flows, Pasadena, California, 3—5 February 1987 1987. F. S. SHERMAN, AND G. M. CoRcos, A numerical simulation of Kelvin-Helmholtz waves of finite amplitude, Journal of Fluid Mechanics, 73 (1976), Pp. 215—240. [44] P. C. PATNAIK, [45] J. PEDL0SKY, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1979. [46] C. PELLAcANI, Shear instabilities in stratified fluids; linear theory, La Rivista Del Nuovo Cimento, 6 (1983), pp. 1—26. [47] S. POND AND G. L. PICKARD, Introductory Dynamical Oceanography, Pergamon Press, New York, 2nd ed., 1983. J. M. CH0MAz, AND P. HuERRE, Propagating Holmboe waves at the interface between two immiscible fluids, Journal of Fluid Mechanics, 266 (1994), pp. 277—302. [48] 0. P0uLIQuEN, [49] On the stability, or instability, of certain fluid motions, Proc. London Math. Soc., 11(1880), pp. 57—70. LORD RAYLEIGH, J. IMBERGER, AND G. M. CoRcos, Turbulence and mixing in stably stratified waters, Annual Review of Fluid Mechanics, 10 (1978), PP. 267—288. [50] F. S. SHERMAN, [51] W. D. SMYTH, Hb1mboe waves, Master’s thesis, University of Toronto, Toronto, Canada, 1986. [52] Instability and Turbulence in Density-Stratified Mixing Layers, PhD thesis, University of Toronto, Toronto, Canada, 1990. , G. P. KLAASSEN, AND W. R. PELTIER, Finite amplitude Holmboe waves, Geophysical and Astrophysical Fluid Dynamics, 43 (1988), Pp. 181—222. [53] W. D. SMYTH, [54] W. D. SMYTH AND W. R. PELTIER, The transition between Kelviri-Helmholtz and Holmboe instability: An investigation of the overreflection hypothesis, Journal of Atmospheric Sciences, 46 (1989), pp. 3698—3720. Bibliography 162 [55] Three-dimensional primary instabilities of a stratified, dissipative, parallel flow, Geophysical and Astrophysical Fluid Dynamics, 52 (1990), pp. 249—261. [56] , Instability and transition in finite-amplitude Kelvin-He lmholtz and Holmboe waves, Journal of Fluid Mechanics, 228 (1991), Pp. 387—415. , [57] H. B. SQUIRE, On the stability of three-dimensional disturbances of viscous flow between parallel walls, Proceedings of the Royal Society A, 142 (1933), pp. 621—628. [58] C. STAQUET, Two-dimensional secondary instabilities in a stably-stratified shear layer, in Proceedings of the 4th International Symposium on Stratified Flows, Greno ble, France, 29 June 2 July 1994. — [59] G. STRANG, Introduction to Applied Mathematics, Wellesley-Cambridge Press, Cambridge, MA, 1986. pp. 453—456. M. W. JOHNsON, AND R. H. FLEMING, The Oceans, Their Physics, Chemistry, and General Biology, Prentic-Hall, Inc., Englewoods Cliffs, N.J., 1942. [60] H. U. SvERDRUP, [61] H. TANAKA, Quasilinear and nonlinear interactions offinite amplitude perturbations in a stably stratified fluid with hyperbolic tangent shear, Journal of the Meteorological Society of Japan, 53 (1975), pp. 1—31. [62] G. I. TAYLOR, Effect of variation in density on the stability of superposed streams of fluid, Proceedings of the Royal Society A, 132 (1931), pp. 499—523. [63] S. A. THORPE, A method of producing a shear flow in a stratified fluid, Journal of Fluid Mechanics, 32 (1968), pp. 25—48. [64] D. J. TRITT0N AND P. A. DAvIEs, Instabilities in geophysical fluid dynamics, in Hydrodynamic Instabilities and the Transition to Turbulence, H. L. Swinney and J. P. Gollub, eds., Springer-Verlag, Berlin, 2nd ed., 1985. [65] J. S. 1973. TURNER, Buoyancy Effects in Fluids, Cambridge University Press, Cambridge, [66] Y. H. WANG, An experimental study of the instability of a stably stratifi ed free shear layer, Journal of Fluid Mechanics, 71(1975), pp. 563—575. [67] R. C. WEAST, ed., CRC Handbook of Chemistry and Physics, CRC Press, Inc., Boca Ratoii, Florida, 1st student ed., 1988. [68] F. M. WHITE, Viscous Fluid Flow, McGraw Hill, New York, 2nd ed., 1991. Bibliography 163 [69] J. D. WOODS, Wave-induced shear instability in the summer thermocline, Journal of Fluid Mechanics, 32 (1968), pp. 791—800. [70] N. Y0NEMITsu, The Stability and Interfacial Wave Phenomena of a Salt Wedge Flow, PhD thesis, University of Alberta, Edmonton, Canada, 1991. [71] S. YO5HIDA, On a mechanism for breaking of interfacial waves, Coastal Engineering in Japan, 20 (1973), pp. 7—15. [72] Z. Zilu, Exchange flow through a channel with an underwater sill, PhD thesis, University of British Columbia, Vancouver, Canada, 1995. Appendix A Dispersion Relation for Piecewise Linear Profiles In section 3.2 we considered an inviscid flow with the following nondimensional back ground velocity aild density profiles. U= 1 l<z<H/2 z —l<z<l —l Pa = —H/2<z<—1 1—1 (A.1) z>—e (1 It was shown that the normal modes approach to linear stability produces the nondimen sional Taylor-Goldstein equation. (U — c)( — 2) — U + 2J(z + = 0, (A.2) For the background flow described by (A.1), (A.2) can be solved exactly. Except at z = —1, —e, and 1, the Taylor-Goldstein equation reduces to = 0. — Thus we can write the solution — (A.3) to (A.2) as e+B 1 A e’ 1 1<z<H/2 e+B 2 A e 2 —E e+B 3 A e 3 —l <z < e+B 4 A e’ 4 —H/2<z<—1 164 < Z < 1 —E (A.4) Appendix A. Dispersion Relation for Piecewise Linear Profiles 165 In order to determine the coefficients A 1 through B , we impose the following bound 4 ary conditions and matching conditions. We require q = 0 at z +H/2 giving 2 = 0, +Bie” 72 Aie° (A.5) e” + 2 4 A 2 e’ = 0. 4 B Next, we require that q be continuous, giving Aie+Bie = e e+B 2 A , (A.6) 6+B e 2 A e = A 2 e+B 3 e, 3 ea = 4 3 ea. 4 e+B 3 A A +B e To derive the matching conditions, we rewrite the Taylor-Goldstein equation (A.2) as [(U — c) — Uz] — (U 2 — c) + 2J6(z + e) On integrating (A.7) across a discontinuity z from z — to z + (A.7) 0 u and letting —* 0, we obtain the following jump conditions. Atz=+1 Z\ [(U — c)qz — Uzq] = 0, (A.8) and at z = —e [(U— c) — U] + 2JU = 0. (A.9) Imposing the jump conditions (A.8) and (A.9) we obtain the following equations for A 1 through B . 4 (—1 + c) [a(Ai 3 (e + c) [a(A (—1 — — c) [a(A 3 — )e 2 A — a(Bi — )e + a(B 2 A 2 — )e 4 A — 3 a(B )e] 2 B = 1 A+B e e, 1 )e] 3 B = — — )ej 4 B e+B 2 (A e), 2 ea. 3 = 3 A +B e (A.10) Appendix A. Dispersion Relation for Piecewise Linear Profiles 166 Equations (A.5), (A.6), and (A.10) form a complete set of equations for A 1 through 4 and can be written in the form B 1 A 2 A 3 A M 4 A = 0 (A.11) 1 B 2 B 3 B 4 B The dispersion relation is obtained by setting det M = 0. Using Maple [10] (a sym bolic programming language) to simplify calculations, the following dispersion relation is obtained. c+a 1 c+a 2 c+a 3 4 +a 4 C = 0 (A.12) where 1 a 2 a 28, = (A.13) J 82 2 c(e + 1 (e 4a 2 — — 1) +4e(sinh2a 3 a = (1 + e 2 1) {(i + — —J sinh(2) e 2 1) ( 2 2a) — 2e° cosh 2ae) 4 e 2 [(i e — — — 2 2a) — e4j 2ccosh2c)}, H_fl 2 +e — (A.14) 2e) — (e 2 — 1) {_(l + 2a) 2+e 2 [(1 4+e +4e(— sinh 2c + 2a cosh 2)}, and — 2)2 — e4j (A.15) Appendix A. Dispersion Relation for Piecewise Linear Profiles 167 _2 4 a = 42(e2H — 1) {_(1 + 2 2a) + 4 e+ e 2 [(1 — 2 2a) — e4j +4e(—sinh2c + 2ccosh2a)} —J 2 ( 3 4 e — 1) +2 [2e(2a2 {_(i + 2)2 — +4e(cosh 2c — 4+e e 2 _1)(1 2 1) + e — 2c sinh 2a)}. — — 2 2a) — e4j 2a) + 2 e ( 2 + 1)] cosh2e (AJ6) Appendix B Maximum Growth Rates for Holmboe’s Equation In chapter 3, we determine the value of the bulk Richardson number, J, at which tran sition from Kelvin-Helmholtz to Holmboe instabilities occurs for the two-layered flow originally considered by flolmboe [23]. In addition, we examine the possibility that in stabilities are initially three-dimensional. To this end, we must determine the maximum growth rate ac for a fixed value of the bulk Richardson number. Note that Smyth et al. [53] found the maximum growth rate for a fixed wave-number. For the case studied by Holmboe we set 0 and let 11 —* cx in the dispersion relation obtained in appendix A for the piecewise linear profile described by (3.9). We obtain the following dispersion relation, 4 + c 2 = 0, a (B.1) where 2 a 4 a —J a+a_ 2 = (B.3) = 1 a = (B.2) — 2o + e 2 2 (B.4) Since we are interested in finding the most amplified growth rate crc, for a fixed value of J, we define o- = —icrc to be the complex growth rate. The above dispersion relation becomes oA 4 o+b 2 +b 168 = 0, (B.5) Appendix B. Maximum Growth Rates for Holmboe’s Equation 169 with = aJ+aa_, (B.6) . 2 aJa (B.7) Equation (B.5) yields two types of instabilities when J when b — 4 4b < > 0, and Kelvin Helmholtz instabilities when b 0: Holmboe instabilities — 4 4b > 2 0 and b < 0. We wish to address the following question: given a bulk Richardson number, J, at what wave number a does the most unstable mode occur (i.e., the largest value of ar)? Due to the nature of the different instabilities, it is necessary to consider the two types separately. B.1 Holmboe Instabilities When the flow is unstable to Holmboe instabilities b four roots of the form a = — 4 4b < 0 and equation (B.5) has +a,. + ia. Solving for a 2 in equation (B.5), one can write (a + ia) 2 = 2+4 —b i/4b 2 — (B.8) Equation B.8 can be used to derive the following equation for ar: a+ + -(b — ) 4 4b = (B.9) 0 In order to determine the maximum growth rate and the wave number at which this growth rate occurs, we differentiate (B.9) with respect to a and set 9 ar/ 3 a resulting equation defines the curve in the a-J plane along which — — r 0. The is maximized: (B.10) 0 Rewriting the above in terms of a and J we obtain J +(1 2 —2a” 1 A — 6a — e+ 2 A 4ae_2) Ju/2 +ah/2(2 —4a 3 A — a) 0 4 2e (B.11) Appendix B. Maximum Growth Rates for Holmboe’s Equation Since it is a quadratic in J1/2, 170 we find that, solving for J, there are two possible solutions. The solution determining the curve along which the growth rate is maximized for a fixed J is given by — 2 (—A — — 3 A 1 4A (B.12) ) — 2A It should be noted that this solution is only valid in the region of Holmboe instabilities. To determine the minimum value of (cr, J) for which (B.12) is valid, one must solve — 4 4b = 0 with J given by (B.12). The minimum value (, J) = (0.488, 0.007) is obtained. Kelvin-Helmholtz Instabilities B.2 For Kelvin-Helmholtz instabilities b equation (B.5) are of the form u = — 4 > 0 and b 4b 2 < 0. In this case, solutions to ar. Thus we can use the above methods on (B.5) directly. The curve along which the growth rate is maximized is defined by: —B + 3 4BB — (B.13) 1 2B where 1 B 16a ( 2 7 2 B 8a — (—4 3 B —1 (2 2 32a + 32a + 8o — 1), — — (B.14) + (2 + — 88c2 2 24a — )7 + (2 3 32a +2 + 32€ — = e a 2 . — — 48)7 + 8a + 3 )7 2 32a , 4 + (1 16a 10cr + 8a 2 + 83)72 + (—2 + 6c (—1 + 2 + 4 )y + -y 2 8a , 5 and 32a + 96o2 — — 6a + 20a 2 2 12a — (B.15) )’y + 4 40€ + 32a )-y + 16a +3 (B.16) Appendix B. Maximum Growth Rates for Holmboe’s Equation 171 By examining the sign of 2 ur/öa it is possible to determine the range of a for which ô the growth rate along the curve given by equation (B.13) is indeed a maximum. The curves defined by (B.12) and (B.13) are shown in figure B.1. 0.50 0.40 0.30 -) 0.20 0.10 0.00 0.0 0.5 1.0 1.5 2.0 cc Figure B.1: Holmboe’s [23] stability diagram Contours are of constant growth rate, act. Dashed line indicate values of a for which the growth rate is maximum for a given J. Appendix C Energy Equations Here we derive the energy equations used for diagnostics of the results from linear stability analysis and from the nonlinear simulations. The equations derived here are only valid for flows with the same boundary conditions as the flow under present consideration. These are essentially the same equations as those derived by Klaassen [27] in his study of Kelvin-Helmholtz instabilities. For a more general description of energy equations for turbulent flows see IVlonin & Yaglom [42]. C.1 Kinetic Energy We start with the nondimensional form of the momentum equations with the Boussinesq approximation (2.12) and (2.13). uu Wt Assuming that p + wu = —p + nw + ww (C.1) —p + — Jp’ + Aw. constant (a valid approximation when 6 p << (C.2) p0), the equation gov erning the evolution of the kinetic energy is found by multiplying (C.1) by u and (C.2) by w and then adding the two together. We obtain UUt = + WWt + u(uu + ww) + w(uu + ww) —(up + wp) — Jwp’ + (uu + wAw). (C.3) Equation (C.3) is the general form for the evolution of kinetic energy in a flow. The 172 Appendix C. Energy Equations 173 last term in (C.3) is usually divided into a viscous diffusion term and viscous dissipation term: nu + w’w + = _[(u)2 ) 2 w 2 + (w) 2+ + (u) (w)2] (C4) viscous dissipation viscous diffusion For simplicity in the following derivations, however, we leave these lumped together and deal with this issue as needed. Since we wish to separate the kinetic energy of the mean flow from that of the per turbations, we write u(x, z, t) w(x,z,t) where, as defined in chapter 4, (j = = ü(z, t) + u’(x, z, t), (C.5) w’(x,z,t), (C.6) (j dx)/L, where L is the length of the domain in the x-direction, indicates that the quantity has been horizontally averaged. Substituting the above expressions for u and w into (C.3), and taking the horizontal average, we obtain the following. 1 (2) + + () + + ‘‘u + () + u (2 + w’w + + + ui (C7) )/2 to be the perturbation kinetic energy and 2 If we define (u’ 2 + w’ /2 2 ü the mean kinetic energy, then (C.7) describes the evolution of the sum of the two kinetic energies. In order to separate the mean and perturbation kinetic energies, we consider the equa tion governing the conservation of momentum in the x-direction. Taking the horizontal average of (C.1) we get (C.8) Appendix C. Energy Equations 174 which when multiplied by ü gives the evolution equation for the mean kinetic energy: (u) + u+ uw’ -uu. — (C.9) To determine the equation for the horizontally averaged perturbation kinetic energy, we subtract (C.9) from (C.7), which gives (u!2 + + a u’ + u’w’w + u’w’u + w’ u w+ 2 +2 + = — — — 77 + JZ -- (u’Au’ + w’.’). (C.1O) An excellent description of the terms in (C.1O) for the homogeneous case can be found in Bradshaw [3]. Equations (C.9) and (C.1O) can be simplified by using the periodicity in the x direction. It is easily shown that = w’w = 0. (C.11) Also, u’p + w’p = u’p + (u + w’jp’ + w’p, by incompressibility, = (u’p’) H- (w’p’), = (w’p’), (C.12) by periodicity in x. Using the above simplifications, the equation for the evolution of the mean kinetic energy can be written as (c.13) and the equation for the evolution of the horizontally averaged perturbation kinetic energy becomes 1 (I2 = + — !2) + (w’p’) + — J+ + w 2 w’ (u’u’ + + w’w’). (C.14) Appendix C. Energy Equations 175 Equation (C.14) is used in chapter 3 to help understand the linear evolution of the flow. In order to understand how the mean flow and perturbations interact in the similla tions described in chapter 4, we study the energy budgets for the horizontally averaged and vertically integrated mean kinetic energy, perturbation kinetic energy, and potential energy. Vertically integrating (C.13) and (C.14) we obtain K + (ü7) = —(uu), (C.15) and K + (u’w’w) + (u’w’u) + (w’w’w) + = — where K = ), K’ = 2 (u ((w’p’)) — 2 + w’ (u’ ), and 2 (a) J() + (u’u’ + w’w’), (.) = f—H/2 (C.16) •dz. Equations (C.15) and (C.16) can be simplified using the boundary conditions = p’ = 0 at z = +H/2 (or equivalently w’ 0 at z = = ‘ = +H/2) as well as the periodicity in the x-direction. It is evident that (w’w’wØ = ((w’p’)) = 0. (C.17) Also, K u’w’w + u’w’u) = u + u’ 2 —(w’ w), integration by parts, 2 = w + u’ 2 (w’ L’) 2 by incompressibility, = w), 2 (w’ by periodicity in x, = 0 w’=Oatz=+H/2. (C 18) Appendix C. Energy Equations 176 and 1 (üw’u) = L 71 1 L H/2 üw’udzdx, I LL H/2 integration by parts, H/2Z, 0 at z 1 1 = L L +H/2 H/2 IH/2 L LI (u’üw’ — u’üu) dzdx, by incompressibility, (C.19) H/2 by periodicity in x, üu’w’dxdz H/2 Finally, Ku’u’ + w’w’) = = _((u)2 + w!2)) + (u)2 — + 2+2 ((u) (u) + 2 (w) + (w) ) 2 (w)2 + (w)2) (C.20) (()2) (C.21) and (aa) - = (()2) Thus (C.15) and (C.16) become (C.22) — and K C.2 = -(u) - J() - + (u)2 + (w)2 + (w)2). (C.23) Potential Energy At a given point in the fluid domain, the potential energy per unit volume is given by gpz, or, in non-dimensional form, Jpz. In order to compute the potential energy of the Appendix C. Energy Equations 177 flow, we consider the nondimensional form of the heat equation (2.9), Pt + UPx + WPz (C.24) RePr Separating the density into a mean and perturbation, p(x, z, t) — (z, t) + p’(x, z, t) (C.25) and using similar expansions for u and w, (C.24) becomes Pt + Pt + UP + + WPZ + Wp = RePr + + pj. (C.26) Multiplying by Jz to obtain the non-dimensional potential energy, taking the horizontal average, and integrating vertically gives (Jzt) + J(z(u’p + w’p)) = RePr (C.27) We can simplify (C.27) using the following: (z(u’p + w’p)) = (zu’p + (u + w)p’ + w’p), by incompressibility, = (z(u’p’) + (w’p’)), = (z(w’p’)), by periodicity in x, = —(7), integration by parts, (C.28) w’ = 0 at z = ±H/2. Equation (C.27) simplifies to J(wp) + RePr’ (C.29) where P = J(z) is the total potential energy of the flow. Unlike the kinetic energy, the potential energy cannot be separated into mean and perturbation parts. Instead, the potential energy combines the effects of both the mean flow and the perturbations. Equations (C.22), (C.23), and (C.29) are used in chapter 4 to compute the energy budget of the flow as the perturbations evolve. Appendix D Period of Perturbation Kinetic Energy 0, linear stability analysis predicts that the phase speed of right and left modes When are different. Using the normal mode approach we can write the velocity perturbation as the sum of the right and left modes, i.e., z, t) where (‘1 and tively and (r) = u)(x, z, t) + u’(x, z, t), = G(r)(z)cosx(r) + S(T)(z)sinx + C(z)cosx + 1 )(z)sinx( S( ) , (D.2) (l) (D.l) are the components corresponding to the right and left modes, respec = — c)t), = — )t). For simplicity we consider only the 1 c horizontal velocity component u’, but the results are valid when the vertical velocity component w’ is also present. The periods of the right and left modes are given by T(r) = 27r/o4r) and T 1 27r/oc!). When the two modes are superimposed, the period of the resulting wave is given by T = lcm(T,T), (D.3) the lowest common multiple of the two periods. In chapter 3, we silperimpose the two modes and examine how the horizontally aver aged perturbation kinetic energy changes with time, giving a better understanding of the linear evolution of the flow. Here, we compute the period of the horizontally averaged perturbation kinetic energy. 178 Appendix D. Period of Perturbation Kinetic Energy = JLU,2dx 1 = 1 = 1 = 2o 179 dx, 2 ) 1 + u’( jL ) + s(r) 2 + [c( + S(r) +2 (C(l)G(r) + ) + c( + x]dx, + cos a(c + 2 (C(r)S(l) + C(1)S(r)) sin (e — — c1))t], (D.4) where L = 2ir/a is the wavelength of the flow (i.e., the period in the x-direction). From (D.4), the period of the horizontally averaged perturbation kinetic energy is 2ir/a(cT) —
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Non-symmetric Holmboe waves
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Non-symmetric Holmboe waves Haigh, Susan Patricia 1994
pdf
Page Metadata
Item Metadata
Title | Non-symmetric Holmboe waves |
Creator |
Haigh, Susan Patricia |
Date Issued | 1994 |
Description | When two flows of different velocity and density meet, a shear layer with a density gradient is formed. Under certain conditions this flow can be unstable. A statically stable stratified shear flow in which the density interface is much thinner than the shear layer thickness can be linearly unstable to two modes with equal growth rates and equal and opposite phase speeds. The superposition of these two modes is called a Holmboe instability. This instability is only possible when the flow is symmetric about the center of the shear layer. We examine the effect of breaking this symmetry by allowing the center of the density interface to be displaced with respect to the center of the shear layer. There are three major components to this study: linear stability analysis, nonlinear numerical simulations, and comparison with laboratory experiments. Linear stability analysis is used to examine the effect of the density interface offset on the overall stability of the flow. Both inviscid theory with piecewise linear background velocity and density profiles, and viscous theory with smooth background profiles are used. As in previous studies, it was found that the growth rate of one mode increases and that of the other mode decreases as the density interface displacement is increased. The precise behaviour depends on the relative thickness of the density interface with respect to the shear layer thickness. For inviscid theory with piecewise linear background profiles, we show that the initial perturbations must be two-dimensional. When the effects of viscosity and diffusion are included, however, it may be possible that the weaker mode is initially three-dimensional. Detailed analysis of the energy transfer in the linear regime indicates that when the background flow loses its symmetry, it is the mode with the larger growth rate that is primarily responsible for the extraction of energy from the mean flow. Two-dimensional numerical simulations are used to examine the nonlinear development of “non-symmetric” Holmboe instabilities. We start by perturbing the flow with the stronger mode predicted by linear theory. We then examine the response of the flow to weaker mode. Finally, we impose both modes. By comparing the development of these three flows we are able to study the interaction between the two modes and the effect of initial conditions on the development of instabilities. Although the initial development of instabilities depends on the initial conditions, this dependence weakens as the density interface offset is increased. Also, preliminary results indicate that long term behaviour of the perturbations are independent of initial conditions. Results of the numerical simulations are compared to both symmetric and nonsymmetric Holmboe instabilities that have been observed in laboratory experiments. Since we are unable to compute the flow at the high Prandtl number of the salt stratified experimental flows, we run the simulation for a range of increasing Prandtl numbers to determine how instabilities of thermally stratified flows (with low Prandtl number) and salt stratified flows differ. Results of these simulations indicate that differences between the experimental and numerical flows can be attributed to the thicker density interface and lower Prandtl number used in the numerical simulations. When the density interface is thicker and the Prandtl number lower, the waves or billows formed by the instabilities are not as sharply defined as in the experiments. |
Extent | 3261191 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080013 |
URI | http://hdl.handle.net/2429/7225 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1995-059707.pdf [ 3.11MB ]
- Metadata
- JSON: 831-1.0080013.json
- JSON-LD: 831-1.0080013-ld.json
- RDF/XML (Pretty): 831-1.0080013-rdf.xml
- RDF/JSON: 831-1.0080013-rdf.json
- Turtle: 831-1.0080013-turtle.txt
- N-Triples: 831-1.0080013-rdf-ntriples.txt
- Original Record: 831-1.0080013-source.json
- Full Text
- 831-1.0080013-fulltext.txt
- Citation
- 831-1.0080013.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080013/manifest