Kelvin-Helmhlotz instabilities in sheared density stratified flows by Mona Rahmani B.Sc. Civil Engineering, AmirKabir University of Technology, 2002 M.Sc. Civil Engineering, Sharif University of Technology, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2011 c© Mona Rahmani 2011 Abstract Kelvin-Helmholtz instabilities are the most commonly studied type of instabil- ity in sheared density stratified flows. Turbulence caused by these instabilities is an important mechanism for mixing in geophysical flows. The primary ob- jectives of this study are the evolution of these instabilities and quantifying the mixing they generate using direct numerical simulations. The results are presented in three chapters. First, the evolution of primary Kelvin-Helmhlotz instabilities in two dimes- nions is studied for a wide range of Reynolds and Prandtl numbers, represent- ing real oceanic and atmospheric flows. The results suggest that some prop- erties of KH billows are predictable by a semi-analytical model. It is shown that a new Corcos-Sherman scale is a useful guide when simulating turbulent KH flow fields. The details of the mixing process generated by the evolution of Kelvin- Helmholtz instabilities as it goes through different stages, is analyzed. As the Reynolds number increases a transition in the overall amount of mixing is found, which is in agreement with previous experimental studies. This transition is explained quantitatively by the entrainment and mixing caused by three-dimensional motions, in addition to those resulted from the two- dimensional growth of the instability. The effect of Prandtl number on mixing is studied to understand the char- acteristics of high Prandtl number mixing events in the ocean; these cases have usually been approximated by low Prandtl number simulations. The increase in the Pradtl number has some significant implications for the evolution of the billow, the time variation of mixing properties, and the overall mixing. ii Preface The three chapters describing the research in this thesis have been submitted or are in preparation for journal publication. We outline the contributions of the co-authors in the following. Chapter 2 has been submitted for publication. The authors of this chap- ter are myself, B. R. Seymour and G. A. Lawrence. I was responsible for performing the numerical simulations, analyzing the data and preparing the manuscript. B. R. Seymour and G. A. Lawrence editted and revised the manuscript. Chapter 3 has been submitted for publication. The authors of this chapter are myself, G. A. Lawrence and B. R. Seymour. The ideas for this research were developed by myself, G. A. Lawrence and B. R. Seymour, and the research was supervised by G. A. Lawrence and B. R. Seymour. I was responsible for conducting the research and preparing the manuscript, with revisions from G. A. Lawrence and B. R. Seymour. Chapter 4 is in preparation for journal publication. I contributed to devel- oping the research with the guidance of G. A. Lawrence and B. R. Seymour. The manuscript was prepared by myself and revised by G. A. Lawrence and B. R. Seymour. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Kelvin-Helmholtz instabilities . . . . . . . . . . . . . . . . . . 2 1.1.1 Theory of shear instabilities . . . . . . . . . . . . . . . 3 1.1.2 Laboratory generation of Kelvin-Helmholtz instabilities 8 1.1.3 Kelvin-Helmholtz instabilities in the field . . . . . . . 10 1.1.4 Numerical simulations of Kelvin-Helmholtz instabilities 12 1.2 Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.1 Mixing transition . . . . . . . . . . . . . . . . . . . . . 14 1.3 Length scales . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 iv Table of Contents 2 The evolution of large and small-scale structures in Kelvin- Helmholtz instabilities . . . . . . . . . . . . . . . . . . . . . . . 21 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.1 Governing equations . . . . . . . . . . . . . . . . . . . 24 2.2.2 Relevant dimensionless numbers . . . . . . . . . . . . 25 2.2.3 Semi-analytical model of Corcos & Sherman (1976) . 25 2.3 Direct numerical simulation . . . . . . . . . . . . . . . . . . . 28 2.4 Pre-saturation structure of the instability . . . . . . . . . . . 30 2.4.1 Dependence on Re0, Pr = 0.72 . . . . . . . . . . . . . 31 2.4.2 Dependence on Pr, Re0 = 400 . . . . . . . . . . . . . 34 2.5 Temporal evolution of the billow . . . . . . . . . . . . . . . . 36 2.5.1 Vorticity depletion . . . . . . . . . . . . . . . . . . . . 38 2.5.2 Billow height . . . . . . . . . . . . . . . . . . . . . . . 39 2.5.3 Saturation . . . . . . . . . . . . . . . . . . . . . . . . 41 2.6 Similarity solutions . . . . . . . . . . . . . . . . . . . . . . . . 42 2.7 Small length scales . . . . . . . . . . . . . . . . . . . . . . . . 48 2.8 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 56 3 The effect of Reynolds number on mixing in Kelvin-Helmholtz instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.2.1 Relevant dimensionless parameters . . . . . . . . . . . 62 3.2.2 Temporal evolution of a KH billow . . . . . . . . . . . 62 3.2.3 Energy partitions and transfer rates . . . . . . . . . . 64 3.3 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . 68 3.4 Mixing transition in DNS . . . . . . . . . . . . . . . . . . . . 70 3.5 Density structure of the billow . . . . . . . . . . . . . . . . . 73 3.6 Energy partitions . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.7 Mixing rate and efficiency . . . . . . . . . . . . . . . . . . . . 80 v Table of Contents 3.8 Comparison with laboratory experiments . . . . . . . . . . . 84 3.9 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 87 4 The effect of Prandtl number on mixing in Kelvin-Helmholtz instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . 92 4.3 Two-dimensional evolution of KH billows . . . . . . . . . . . 93 4.4 Three-dimensional evolution of the billow . . . . . . . . . . . 97 4.5 Mixing and energy partitions . . . . . . . . . . . . . . . . . . 99 4.6 Overall mixing . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.7 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 108 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.2 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . 111 5.3 Future research . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 vi List of Tables 2.1 Description of the two-dimensional simulations carried out in chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 Description of the three-dimensional simulations carried out in chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Description of the simulations carried out in chapter 3 . . . . 69 4.1 Description of the two and three-dimensional simulations car- ried out in chapter 4 . . . . . . . . . . . . . . . . . . . . . . . 94 vii List of Figures 1.1 Density field of a series of Kelvin-Helmholtz billows. . . . . . . 3 1.2 Profiles of the background velocity and density . . . . . . . . . 6 1.3 Stability diagrams . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Density structure of turbulent billows . . . . . . . . . . . . . . 15 2.1 Configuration of braid and core regions . . . . . . . . . . . . . 27 2.2 Vorticity structure for Pr = 0.72, and different Re0 . . . . . . 35 2.3 Vorticity structure for Re0 = 400, and different Pr . . . . . . 37 2.4 Comparison of DNS and CS76 model ratio of braid to total circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5 Comparison of DNS and CS76 model half billow heights nor- malized by the wavelength . . . . . . . . . . . . . . . . . . . . 40 2.6 Dimensionless billow height, braid shear, and strain rate at sat- uration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.7 Variation of the ratio of the asymptotic half density interface thickness at saturation to the wavelength with Re0 and Pr . . 45 2.8 Comparison of the similarity solution and DNS braid centerline vorticity at saturation . . . . . . . . . . . . . . . . . . . . . . 47 2.9 Snapshot of log10χ, the scalar variance dissipation rate . . . . 50 2.10 Spectra of scalar variance at maximum intensity of turbulence 52 2.11 Variation of different length scales with the Reynolds number . 54 2.12 The percentage of scalar dissipation not captured when exclud- ing length scales below LB and LCS . . . . . . . . . . . . . . . 55 2.13 Snapshot of secondary instabilities . . . . . . . . . . . . . . . 55 viii List of Figures 3.1 The density structure of the KH billow at different stages . . . 64 3.2 A schematic of energy partitions and transfers between them . 67 3.3 Variation of the amount of mixing gained at the end of the KH life-cycle with the Reynolds number . . . . . . . . . . . . . . . 70 3.4 Comparison of the time evolution of the gain in the background potential energy . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.5 Density field at t∗3, maximum intensity of turbulence . . . . . . 74 3.6 Time evolution of the energy partitions normalized by K0 . . . 78 3.7 Time evolution of the rate of mixing and turbulent viscous dis- sipation rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.8 Mixing efficiency . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.9 Variation of mixing with the outer-scale Reynolds number from laboratory measurements and DNS . . . . . . . . . . . . . . . 85 4.1 Snapshots of the density structure of the billow . . . . . . . . 95 4.2 Snapshots of the span-wise density structure . . . . . . . . . . 99 4.3 Time evolution of the energy partitions normalized by K0 for Re0 = 300 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Time evolution of the energy partitions normalized by K0 for Re0 = 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.5 Variation of tD with Pr . . . . . . . . . . . . . . . . . . . . . 104 4.6 Time variation of φM from two- and three-dimensional simula- tions at Re0 = 300 . . . . . . . . . . . . . . . . . . . . . . . . 106 4.7 Time variation of φM from two- and three-dimensional simula- tions at Re0 = 400 . . . . . . . . . . . . . . . . . . . . . . . . 106 4.8 Variation of the total amount of mixing with Pr . . . . . . . . 107 4.9 Variation of M/K0 with Re0 for Pr = 1 and Pr = 9 . . . . . 108 ix List of Symbols A Domain area α Wavenumber of the instability mode c Complex phase speed ci Imaginary part of the complex phase speed cr Real part of the complex phase speed δ Shear layer thickness δ0 Initial shear layer thickness δb Vorticity interface thickness of the braid δP Chemical product thickness in experiments ε Rate of viscous dissipation of kinetic energy ε Rate of viscous dissipation of mean kinetic energy ε ′ Rate of viscous dissipation of turbulent kinetic energy f Perturbation quantity f̂ Perturbation quantity eigenfunction F Arbitrary function for similarity solution of density φ Perturbation streamfunction φ̂ Perturbation streamfunction eigenfunction φ ′ b Turbulent buoyancy flux φi Rate of diffusion of the background stratification φM Rate of mixing g Gravitational acceleration g ′ Reduced gravitational acceleration G Arbitrary function for similarity solution of vorticity γ Strain rate x List of Symbols γ∗ Dimensionless strain rate γK Kolmogorov strain rate Γ Total circulation Γb Braid circulation Γc Core circulation (chapter 2) Γc Cumulative flux coefficient (chapter 3) Γi Instantaneous flux coefficient h Density interface thickness h0 Initial density interface thickness hb Density interface thickness of the braid H Billow height η Direction normal to the braid i Imaginary unit J Bulk Richardson number k Streamwise scalar spectra wavenumber k∗ Dimensionless streamwise scalar spectra wavenumber k̂ Unit vertical vector K Kinetic energy K Mean kinetic energy K ′ Turbulent kinetic energy K0 Initial kinetic energy K2d Two-dimensional kinetic energy K3d Three-dimensional kinetic energy κ Diffusivity of density L Thickness of a diffusing sheet LB Batchelor length scale LCS Corcos-Sherman length scale LD Minimum dissipation length scale LK Kolmogorov length scale LM Maximum dissipation length scale xi List of Symbols LT Taylor length scale L ′ T Equivalent Taylor length scale in scalar field Lx Domain length in the streamwise direction Ly Domain length in the spanwise direction Lz Domain length in the vertical direction M Mixing n Number of rolls in the billow N2 Squared buoyancy frequency Nx Number of grid points in streamwise direction Ny Number of grid points in spanwise direction Nz Number of grid points in vertical direction ν kinematic viscosity P Potential energy Pa Available potential energy Pb Background potential energy P Pressure field P ′ Perturbation pressure P Hydrostatic background pressure Pr Prandtl number R Scale ratio Re Reynolds number Re0 Initial Reynolds number Reδ Outer-scale Reynolds number Ri Gradient Richardson number r Velocity ratio of streams in mixing layer experiments ρ Density field ∆ρ Density difference ρ ′ Perturbation density ρ0 Reference density ρ(z) Background density stratification xii List of Symbols S Shear S Strain tensor σ Direction aligned with the braid ψ Rate of shear production/destruction (chapter 3) ψ Magnitude of scalar spectra (chapter 2) ψ∗ Dimensionless magnitude of scalar spectra t Time t∗ Dimensionless time t∗1 Dimensionless time of end of the first phase t∗2 Dimensionless time of end of the second phase t∗3 Dimensionless time of end of the third phase t∗4 Dimensionless time of end of the fourth phase t∗M Dimensionless time of end of the unmixed phase t∗s Dimensionless time of saturation t∗u Dimensionless time of first roll-up θ Angle of inclination of the braid U(z) Background shear flow ∆U Velocity difference u Velocity component in the streamwise direction u Velocity field vector u Mean one-dimensional velocity u2d Two-dimensional velocity u3d Three-dimensional velocity u ′ Perturbation velocity in the streamwise direction uK Kolmogorov velocity scale v Velocity component in the spanwise direction Ω Two-dimensional scalar vorticity w Velocity component in the vertical direction w ′ Perturbation velocity in the vertical direction x Streamwise direction xiii List of Symbols ∆x Grid spacing in streamwise direction χ Rate of dissipation of scalar y Spanwise direction ∆y Grid spacing in spanwise direction z Vertical direction ∆z Grid spacing in vertical direction xiv Acknowledgments During the course of my studies at UBC, I had the great opportunity to interact with many wonderful people, to whom I am indebted for their remarkable support and encouragement. First, I would like to thank my supervisors Greg Lawrence and Brian Seymour for their limitless support and advice, plenty of patience and benevolent personality. Over these years my scientific approach towards research problems has greatly improved owing to their advice. I am also grateful to my supervisory committee Neil Balmforth and Bernard Laval. The discussions I had with them and their valuable comments has enriched this thesis. I would also like to thank Roger Pieters for his kind help with different aspects of the work, and feedback on my research. The numerical code used in this thesis was developed by Kraig Winters and Bill Smyth, to whom I am deeply grateful. I had the opportunity to learn about many fascinating Fluid Mechanics research projects and get feedback on my research through my interaction with the Environmental Fluid Mechanics group at UBC. A special thanks goes to Jeff Carpenter, Ted Tedford, Anirban Guha, Andrew Hamilton, Violeta Martin and Joel Atwater for their help with different matters during my stay at UBC. Attending the Complex Fluids group meetings directed by Neil Balmforth gave me the chance to hear about many neat problems in Fluid Mechanics and get feedback from brilliant people. Finally, the support and love from my family and friends have been a great source of energy and encouragement for me. I cherished the love from my parents, my sister’s intelligence, my brother’s sense of humor, and the good company of my friends after many long days. xv Dedication To my parents, sister and brother xvi Chapter 1 Introduction In the dynamics of the oceans, atmosphere, lakes and rivers a large num- ber of physical phenomena are present. These phenomena occur at a wide range of scales, from oceanic circulation down to small-scale mixing events that can occur over a few millimeters. The processes at different scales have a close interaction with each other, so that small-scale mixing relies on the large-scale dynamics of the flow, and large-scale processes are influenced by mixing and dissipation at small scales (Wunsch & Ferrari, 2004). Therefore, for fluid environments it is fundamental to understand the processes that result in small-scale mixing. Such understanding addresses many scientific and en- gineering concerns, such as transport of contaminants, climate change effects, and changes in water quality in lakes and oceans. An important process that generates a significant amount of mixing in the interior ocean and atmosphere is the generation of shear instabilities in density stratified flows. Through the interaction of the instability with the background shear, energy is extracted from the large-scale motions to induce small-scale mixing. The best known type of shear instabilities are Kelvin-Helmholtz (KH) waves. These instabilities form by rolling up the interface between different layers of fluid and creating billows of rotating fluid. Breaking of these billows generates turbulent patches and mixing, clearly observed in nature and the laboratory (Thorpe, 2005). While different aspects of the generation, growth, and breaking of KH waves have been extensively studied before (Thorpe, 2005), the process of transition of the initially laminar flow to a turbulent stage and the behavior of mixing through this transition is still not completely understood. Especially 1 1.1. Kelvin-Helmholtz instabilities quantification of mixing properties in turbulent flow fields has been a chal- lenging problem (Ivey et al., 2008). Recent improvements in computational resources and laboratory and field measurement techniques have allowed de- tailed investigations of the small-scale processes. The present thesis is a con- tribution to this subject using numerical simulations. In this chapter, first a background on KH instabilities is given by reviewing previous studies. Then some ideas on mixing and small scales are introduced. Finally, an overview of the thesis is stated. 1.1 Kelvin-Helmholtz instabilities KH instabilities form as a series of vortices that structure periodic rolling billows, as shown in figure 1.1. The region between billow cores forms a tilted interface with a thin thickness, called the braid. The billows naturally tend to amalgamate together due to vortex pairing. In the process of vortex pairing, two vortices interact to first roll around each, and then form a single vortex with double the strength of each initial vortex (Winant & Browand, 1974; Pierrhumbertt &Windall, 1982). The two-dimensional evolution of KH billows is followed by the growth of three-dimensional secondary instabilities, and then transition to turbulence. Different methods have been used to study KH instabilities. Theoretically, the onset of KH waves and their initial stages of growth are described by linear stability analysis. Laboratory generation of KH waves allow detailed examination of the evolution of KH billows through the entire life-cycle. Field measurements describe KH billows as part of a larger-scale process. Numerical simulations provide a tool for accurate quantification of KH flow fields, however under idealized conditions. Hence, these approaches complement each other in KH instability studies. In the following, previous studies using each technique are briefly reviewed. 2 1.1. Kelvin-Helmholtz instabilities vortex pairing braid billow core Figure 1.1: The evolution of KH billows and configuration of braid, core and vortex pairing. The snapshot shows the density field of a series of KH billows from a numerical simulation. 1.1.1 Theory of shear instabilities Hydrodynamic stability of sheared density stratified flows is a classic problem of fluid mechanics (e.g. Drazin & Reid, 1981). In this problem a background parallel shear flow of U(z) and vertical density stratification of ρ(z) are con- sidered. This flow is perturbed with a small two-dimensional perturbation velocity of (u ′ , w ′ ), density of ρ ′ , and pressure of P ′ , so that the total two- dimensional velocity field, (u, w), density field, ρ, and pressure field, P, are u = U + u ′ , w = w ′ , ρ = ρ+ ρ ′ , P = P + P ′ , (1.1) where P is the hydrostatic background pressure, and dP/dz = −ρg. For this two-dimensional flow, in a Cartesian system of coordinates of (x, z) and time frame of t, the Boussinesq Navier-Stokes equations express the conservation of horizontal and vertical momentum as ut + uux + wuz = −Px ρ0 + ν(uxx + uzz), (1.2) wt + uwx + wwz = −Pz ρ0 − g ρ ρ0 + ν(wxx + wzz), (1.3) 3 1.1. Kelvin-Helmholtz instabilities where ρ0 is a reference density, ν the kinematic viscosity, g the gravitational acceleration, and the subscripts denote partial differentiation. In these equa- tions, using Boussinesq approximations (e.g. see Kundu & Cohen, 2004), it is assumed that the variation in density is negligible except when ρ is multiplied by g. For incompressible fluids the conservation of mass is stated as ux + wz = 0. (1.4) Assuming a linear equation of state, that relates scalars such as temperature and salinity to density, the conservation of a scalar(s) is expressed in terms of ρ ρt + uρx + wρz = κ(ρxx + ρzz), (1.5) with κ being the molecular diffusivity of the scalar. The stability of the background flow to the perturbations is linearly an- alyzed by substituting equations 1.1 in equations 1.2 to 1.5, and neglect- ing the terms with the product of perturbations. A conventional way for solving the resulted linear equations is the method of normal modes. In this method each perturbation quantity, f ′ , is assumed to have the from f ′ (x, z, t) = f̂(z)eiα(x−ct), with i being the imaginary unit, α the real wavenum- ber of the mode, and c the complex phase speed. Substituting these forms in linear equations of motion gives an eigenvalue problem, with the eigenvalue c and eigenfunction f̂ . For each mode, c = cr + ici, and αci is the exponen- tial growth rate and cr the phase speed of the mode. Therefore, the mode is unstable if ci > 0. For a non-diffusive and inviscid fluid, the eigenvalue equations obtained using the steps described above, simplify to Taylor-Goldstein equation (Taylor, 1931; Goldstein, 1931) d2φ̂ dz2 + [ N2 (U − c)2 − 1 (U − c) d2U dz2 − α2 ] φ̂ = 0, (1.6) 4 1.1. Kelvin-Helmholtz instabilities where φ̂ is the eigenfunction for the perturbation streamfunction, and N2 is the squared buoyancy frequency, defined as N2(z) = −g/ρ0dρ/dz. By solving the Taylor-Goldstein equation for a general background velocity of U(z) and ρ(z), Miles (1961) and Howard (1961) showed that a necessary condition for instability of the flow is that the gradient Richardson number, Ri = − g ρ0 ρz U2z , (1.7) of the flow is below 1/4 somewhere in the flow field. The gradient Richardson number measures the ratio of the strength of density stratification to that of shear, and Ri = 1/4 is known as Miles-Howard critical value. This criterion indicates that if Ri > 1/4 everywhere in the flow, the shear is not strong enough to overcome the stratification, and therefore the flow is stable. The exact solution of the Taylor-Goldstein equations depends on the form of the background shear and stratification. In this thesis we are interested in two streams of different velocities, forming a smooth shear layer between them, with the top layer having lighter density, forming a smooth stable density stratification, as shown in figure 1.2. This flow is a good representation of geophysical shear flows (e.g. Farmer & Armi, 1999) and mixing layers in the laboratory (e.g. Koop & Browand, 1979), and is commonly expressed as (e.g. Hazel, 1972) U(z) = ∆u 2 tanh( 2z δ0 ), and ρ̄(z) = −∆ρ 2 tanh( 2z h0 ), (1.8) with δ0 being the vorticity interface thickness , h0 the density interface thick- ness, ∆u the horizontal velocity difference and ∆ρ the density difference. This flow forms four important dimensionless numbers Re0 = δ0∆u ν , J = ∆ρgδ0 ρ0∆u2 , P r = ν κ , and R = δ0 h0 , (1.9) the initial Reynolds number, bulk Richardson number, Prandtl number and 5 1.1. Kelvin-Helmholtz instabilities Figure 1.2: Profiles of the background velocity and density. scale ratio of the thickness of the shear layer to that of the diffuse density layer, respectively. The Taylor-Goldstein equation for these profiles has been solved semi-analytically by Hazel (1972); Haigh (1995); Alexakis (2005); Carpenter et al. (2010a). Also piecewise profile approximations of the smooth profiles were used by Holmboe (1962); Lawrence et al. (1991); Baines & Mitsudera (1994); Haigh & Lawrence (1999); Carpenter et al. (2010a) for the analytical solution. These studies show that variations in J and R significantly influence the nature and characteristics of the instability. As J increases from J = 0, the case of homogeneous flow, the flow becomes more stable and the growth rate of the KH instability decreases. For sufficiently high J , there is a transition form KH to Holmboe waves (Holmboe, 1962). The Holmboe instability consists of two oppositely travelling waves, reaches smaller heights and generates much less vigorous turbulence compared to KH instability. The increase in the scale ratio, R, can also cause a transition from KH to Holmboe waves. Alexakis (2005) found that for smooth hyperbolic tangent profiles, R > 2 is a necessary condition for the formation of Holmboe waves. 6 1.1. Kelvin-Helmholtz instabilities 0 0 0 0. 04 0.04 0. 08 0. 12 0.16 α δ0 J • (a) 0.5 1 1.5 0.02 0.04 0.06 0.08 0 0. 04 0.0 4 0. 08 0. 08 0. 12 0.16 α δ0 J • (b) 0.5 1 1.5 0.02 0.04 0.06 0.08 Figure 1.3: Stability diagrams for Re0 = 300, and (a) R = 3 and (b) R = 26.5. The contours show the growth rate, the dashed line the boundary of transition from KH to Holmboe wave, and the dot the fastest growing mode at J = 0.03. In this thesis a single bulk Richardson number of J = 0.03 (and 0.035 for some cases) is considered, and the scale ratio varies from R = 1 to R = 26.5. High values of R represent stratification with higher Pr, for which the rate of diffusion of the density interface is considerably smaller than that of the shear layer. Scaling of the diffusion times shows R = Pr1/2. In figure 1.3, the growth rates and the boundary of transition from KH to Holmboe waves are shown on J − αδ0 planes (αδ0 is the dimensionless wavenumber) for Re0 = 300 and R = 3 and 26.5. This analysis is obtained by numerically solving the linear stability equations which also include viscosity and diffusivity. The figure shows that J = 0.03 is small enough to ensure KH instabilities even at high values of R. The KH mode region significantly shrinks as R increases from 3 to 26.5. The implications of increasing R are also a decrease in the growth rate (from 0.13 at R = 3 to 0.12 at R = 26.5) and a decrease in the wavenumber (from 0.74 at R = 3 to 0.64 at R = 26.5) for the fastest growing modes at J = 0.03 shown in figure 1.3. Linear stability analysis provides valuable information about the type of the instability, its growth rate and wavelength. However, as the instability grows, the amplitude of the perturbation quantities become large, the nonlinear terms 7 1.1. Kelvin-Helmholtz instabilities cannot be neglected any more. To study the nonlinear stage of the growth of the instability semi-analytical models or numerical simulations are used. In chapter 2 of the thesis we study a semi-analytical model for the nonlinear growth of KH instability, proposed by Corcos & Sherman (1976) based on an earlier work by Batchelor (1972). 1.1.2 Laboratory generation of Kelvin-Helmholtz instabilities KH instabilities have been generated in the laboratory using different experi- mental setups. Thorpe (1968, 1971) and Atsavapranee & Gharib (1997) gener- ated temporally evolving KH instabilities by tilting a tube filled with a dense fluid in the bottom and light fluid in the top. Spatially evolving KH instabili- ties formed downstream of a splitter plate separating two streams of flow with different velocities in the experiments of Brown & Roshko (1974); Breidenthal (1981); Lawrence et al. (1991); Schowalter et al. (1994). Spatially accelerating KH billows were generated in the flow of a dense fluid over a sill in a contracted channel in the experiments of Pawlak & Armi (1998) and Hogg & Ivey (2003). The primary KH instability grows as a wave on the interface between the two layers, then rolls up the interface and forms a series of rotating billows (e.g. Schowalter et al., 1994; Atsavapranee & Gharib, 1997). The sheared interface becomes unstable to KH waves when the gradient Richardson number is below the Miles-Howard critical value of 0.25 (Thorpe, 1973). The initial growth rate and wavelength of laboratory KH waves are in good agreement with the predictions of linear theory (Thorpe, 1971; Scotti & Corcos, 1972; Lawrence et al., 1991). After growing to a maximum height, the primary KH instability becomes susceptible to different types of secondary instabilities, vividly observed in the laboratory. Streamwise vortices on the braid and in the convectively unstable regions inside the billow are the most important types of secondary instabil- ities and the main mechanism for the transition to turbulence (e.g. Bernal 8 1.1. Kelvin-Helmholtz instabilities & Roshko, 1986; Schowalter et al., 1994). In homogeneous or weakly strati- fied KH billows, counter-rotating streamwise vortices are caused by the vortex stretching of the braid region, a mechanism explained in theory by Lin & Cor- cos (1984). This streamwise vortex structure was first indicated by the vortex streaks in experiments of Brown & Roshko (1974); Konrad (1976); Breidenthal (1981), and studied in more detail in the experiments of Jimenez (1983); Bernal & Roshko (1986); Lasheras & Choi (1988); Schowalter et al. (1994). As the stratification increases the convective instability inside the core of the billow becomes more dominant in forming the streamwise convective rolls (Schowalter et al., 1994). Secondary streamwise vortices significantly enhance the entrain- ment and mixing (Brown & Roshko, 1974; Bernal & Roshko, 1986). Due to the baroclinic generation of vorticity, the braid region becomes unstable to secondary shear instabilities. These instabilities have the same structure as the primary KH billow and were observed in the experiments of Atsavapranee & Gharib (1997). Vortex pairing is another important feature in the series of KH billows generated in the laboratory, and enhances the entrainment inside the billow region (Winant & Browand, 1974; Patterson et al., 2006). Pairing is severely suppressed by the increase in stratification (Atsavapranee & Gharib, 1997). Three-dimensional vortex tubes and knots are other features observed by Thorpe (1987). The growth of these different secondary structures and their interaction with each other and with the primary KH billow, creates chaotic small-scale structures in the flow and a transition to the turbulent phase. While the small-scale motions develop thoroughly through the billow, the experiments of Brown & Roshko (1974) and Dimotakis & Brown (1976) suggest that the structure of the primary billow plays an important role in the transition to turbulence and dynamics of the turbulent flow. Mixing caused through different phases of KH instabilities has been quan- tified in the laboratory in a limited number of studies. Density and velocity profile measurements of Thorpe (1973) indicate that only about 25% of the 9 1.1. Kelvin-Helmholtz instabilities energy extracted from the mean flow by the instability is converted to mixing. Koop & Browand (1979) measured the spatial evolution of mixing downstream a splitter plate. Their results showed a slow initial increase in mixing, a rapid rise in mixing in the turbulence decay region of the flow, and an asymptotic final amount of mixing at late stages of the evolution. Mixing measurements of Konrad (1976); Breidenthal (1981); Koochesfahani & Dimotakis (1986) in the laboratory, at a fixed location near the final stages of the evolution of KH billows, revealed a transitional range of Reynolds numbers for the amount of mixing. The time evolution of entrainment and mixing in shear flows in the laboratory measured by Patterson et al. (2006) showed a strong dependence on the stage of evolution of the KH billow. 1.1.3 Kelvin-Helmholtz instabilities in the field KH billows have been observed in many field studies in the ocean and at- mosphere, and found to be an important mechanism for the generation of turbulence and mixing in geophysical flows. The vertical extent of billows can reach as high as 20 m in the ocean (DeSilva et al., 1996), and 1 km in the atmosphere (Luce et al., 2010). The destabilizing shearing of the density stratification can be provided by internal waves, flow over topography or tidal flows in salt wedge estuaries. KH roll-up induced by internal waves and the turbulent patches gener- ated by their breakdown was first observed by Woods (1968) on the summer thermocline in the Mediterranean Sea. Thorpe (1978) describes KH billow structure generated by wind forcing on a thermocline near the surface of Loch Ness. High frequency temperature fluctuations on the internal waves observed by Marmorino (1987) and Herbert et al. (1992) are attributed to KH instabili- ties. In deep ocean, Haren & Gostiaux (2010) found KH instabilities caused by internal waves on a sloping sea floor. KH billows formed on internal solitary waves were observed by Moum et al. (2003) and Lamb & Farmer (2011). These billows normally appeared near the trough of solitary waves, where shearing 10 1.1. Kelvin-Helmholtz instabilities was maximum, and their turbulent breakdown generated a significant amount of mixing and dissipation. Farmer & Armi (1999) illustrate how tidal flow over topography in density stratification generates KH instabilities, and explain how the mixing caused by these instabilities influences the establishment of the large-scale flow. Shear instabilities, with KH-like wave structure, on the thermocline of salt wedges were detected by Geyer & Smith (1987); Geyer & Farmer (1989); Tedford et al. (2009); Geyer et al. (2010). In these studies the maximum shearing occurred during the ebb phase of the tidal flow, therefore the thermocline was most unstable during this phase. Echo sounding images of Geyer et al. (2010) clearly showed the occurrence of secondary KH instabilities on the braids at very high Reynolds numbers. In the atmosphere, the evidence of billowing caused by KH instabilities was shown by Ludlam (1967); Gossard et al. (1970); Browning (1971); Luce et al. (2010). A great number of field measurements suggest that shear instabilities form in regions of the flow with the gradient Richardson number, Ri, close to or below 0.25 (e.g. Geyer & Smith, 1987; Tedford et al., 2009; Geyer et al., 2010), in accordance with Miles-Howard criteria. DeSilva et al. (1996) showed the agreement between non-dimensional length and time scales of KH instabil- ities observed in the field and laboratory. Smyth et al. (2001) showed the consistency between statistics of turbulence generated by the KH billows from measurements on a thermocline and DNS. Therefore laboratory and DNS gen- eration of KH instabilities, if performed at sufficiently high Reynolds numbers, can be realistic representations of instabilities in the ocean. Mixing caused by KH billows has been quantified in the field measurements of Seim & Gregg (1994, 1995); Smyth et al. (2001); Moum et al. (2003). A key question to answer has been how much of the total energy extracted by the KH billow from the mean flow is spent to induce mixing. A large number of observed turbulent events in the ocean (reviewed by Inoue & Smyth, 2009) suggest that this ratio has a universal value. In chapter 3 of this thesis we 11 1.1. Kelvin-Helmholtz instabilities examine this ratio, along with other energy budgets of the evolution of the flow. 1.1.4 Numerical simulations of Kelvin-Helmholtz instabilities Numerical simulations of the roll-up of a vortex sheet by KH instability were first performed by Rosenhead (1931). Since then numerical simulations have been used to simulate large and small-scale features of the evolution of KH in- stabilities and quantify the characteristics of generated turbulence and mixing. In direct numerical simulations (DNS) the domain is resolved to the smallest scale of motion, and therefore simulation of the exact behavior of turbulent motions is possible. The growth of the primary KH instability and the pre-turbulent struc- tures developed in the flow has been simulated numerically in several studies (e.g. Patnaik et al., 1976; Klaassen & Peltier, 1985a; Staquet, 1995). The susceptibility of the primary KH billow to secondary instabilities, as observed in the laboratory, has been shown in numerical simulations (e.g. Klaassen & Peltier, 1985b, 1991; Caulfield & Peltier, 1994; Staquet, 1995; Smyth, 2003). The nature of secondary instabilities, their evolution and contribution to trig- gering turbulence is examined in detail by Caulfield & Peltier (2000); Peltier & Caulfield (2003); Smyth (2003). Simulating the turbulence and mixing in KH flow field is an important ap- plication of DNS. Smyth (1999) examined the process of dissipation of scalars at small scales in turbulence generated by KH instabilities. The evolution of important length scales that characterize turbulence in KH flow field were simulated by Smyth & Moum (2000) and Smyth et al. (2001). Time evo- lution of mixing, and quantities that characterize mixing behavior, such as mixing efficiency, are examined over entire or partial life-cycle of KH billows in numerous numerical simulations (for instance Scinocca, 1995; Cortesi et al., 1999; Caulfield & Peltier, 2000; Staquet, 2000; Carpenter et al., 2007; Inoue & 12 1.2. Mixing Smyth, 2009). However, these studies either focus on the effect of stratification on mixing (e.g. Caulfield & Peltier, 2000), or consider a very limited range of Re0 and Pr in mixing events (e.g. Smyth et al., 2001; Cortesi et al., 1999; Sta- quet, 2000). So a comprehensive investigation of the effect of Reynolds and Prandtl number on mixing induced by KH instabilities is still missing from the literature. Chapters 3 and 4 of this thesis focus on studying the effects of Reynolds and Prandtl number on mixing. In this thesis the DNS code described by Winters et al. (2004) and Smyth et al. (2005) is used. This code solves Boussinesq Navier-Stokes, continuity, and diffusion-advection equations in three space dimensions and time. Full Fourier transforms in horizontal directions and half-range sine or cosine Fourier transforms in the vertical direction are used to convert the set of partial dif- ferential equations to ordinary differential equations. The integration in time is performed using a third-order Adams-Bashforth method. The code is run on distributed memory multiprocessor computers. 1.2 Mixing The last stage of mixing in any flow field is diffusion of scalars at molecu- lar scales. The preceding stages, as described by Dimotakis (2000), are the entrainment, and stirring. Through the entrainment process, fluid from non- turbulent regions of the flow are “engulfed”or “entangled”by large scale eddies (Brown & Roshko, 1974). Through the stirring process, the entrained fluid is distributed in the turbulent region by intermediate size eddies, generating more interfacial area between the scalars. Studying mixing in density stratified flows has the advantage that the density serves as a label to measure mixing. Therefore to study mixing it suffices to obtain the density field of the flow. However, as pointed out by Dimotakis (2005), the density stratification couples to the dynamics of the flow and makes the mixing behaviour more complicated compared to passive 13 1.2. Mixing mixing. The dynamics of large, intermediate and small scale eddies that drive the mixing process are significantly under the influence of the Reynolds number. In an extensive number of studies, this effect has indicated a transition in mixing, a topic briefly reviewed in the following. 1.2.1 Mixing transition The “mixing transition” is an interesting feature of turbulent flows, previously observed in different types of flows such as shear layers, wakes, jets and pipe flows. This transition refers to a rapid increase in the amount of mixing in the turbulent flow due to an abrupt change in the characteristics of the flow over a narrow range of transitional Reynolds numbers. Beyond this transition zone the dependence of mixing on the Reynolds number is weak. The transition starts at the Reynolds number, at which the energizing effect of inertial forces in the flow just overcomes the dissipating effect of viscous forces to develop chaotic three-dimensional motion throughout the disturbed region of the flow. As the Reynolds number increases, these motions become more energetic and more cascade of energy to small-scale eddies occurs. The development of a wide range of sizes of eddies in the flow facilitates the stirring process of scalars and therefore enhances the mixing. This accounts for the rapid increase in the amount of mixing in the transitional range of Reynolds numbers. After sufficiently increasing the Reynolds number, the spectrum of length scales includes all sizes of eddies from the largest scale of motion to the smallest scales below which the action of molecular diffusion dissipates all motions. At this limit the turbulence is fully developed and further increases in the Reynolds number do not cause more mixing. Several studies have shown that Reynolds number independent mixing is related to a fully developed turbulence (see Dimotakis, 2000, 2005). An example of under developed and fully developed turbulence in shear layers is illustrated in figure 1.4 by showing the density structure of turbulent billows. By reviewing a wide 14 1.2. Mixing (a) (b) Figure 1.4: Density structure of turbulent billows for: (a) an under developed turbulence, (b) a fully developed turbulence. number of studies of the transition in different turbulent flows, Dimotakis (2000) concludes that a fully developed turbulence is reached at an outer-scale Reynolds number (explained in more detail in section 2) greater than order of 104. In shear layers, the mixing transition was first noted by Konrad (1976) in a spatially evolving gas-phase shear layer, formed by two parallel layers flowing at different velocities downstream of a splitter plate. The amount of mixing was measured by the amount of chemical product between two reacting chemicals added to the two streams. He observed a rapid 25% increase in the amount of mixing beyond a critical Reynolds number where three-dimensional instabilities grew in the braid region of the billow and caused small-scale disor- ganized motions to develop in the flow. The mixing transition in a liquid-phase mixing layer was measured by Breidenthal (1981) and Koochesfahani & Dimo- takis (1986) by measuring the amount of chemical product between the two layers. For the same velocity ratio of the two streams, the mixing transition observed in these two studies occured over the same range of Reynolds num- 15 1.3. Length scales bers observed by Konrad (1976). The one order of magnitude increase in the amount of mixing in the measurements of Breidenthal (1981) and Koochesfa- hani & Dimotakis (1986) was associated with the development of small-scale turbulent motions in the flow. The effect of secondary streamwise vortices on the mixing transition was also confirmed in the experiments of Bernal & Roshko (1986). Understanding the transition in mixing and ranges of Reynolds numbers over which this transition occurs is important for mixing studies in turbulent flows in the ocean and atmosphere. As these flows usually have very high Reynolds numbers (Re > 105, Re defined based on a typical large length scale of the flow), they can safely be assumed to exhibit post-transition mixing behaviour. Hence, any numerical simulation or laboratory experiment that is performed at sufficiently high Reynolds number to ensure a post-transition mixing, can be a reasonable representation of the real-life flow. 1.3 Length scales To study mixing in turbulent flows, it is essential to understand different length scales of the flow. Associated with each length scale, there is a velocity scale, based on which a turn-over time scale can be defined for each size of the eddy. The large length scale is a measure of the size of the largest eddy existing in the flow. In KH billows this large scale is usually measured by the vertical extent of the billow, or the wavelength of the instability. Following the former convention, δ is considered as the large scale here. The total velocity difference in the shear layer, ∆U , is the large velocity scale. The largest scales of motion are also referred to as “outer scales” (e.g. Buch & Dahm, 1996). In turbulent flows large scale eddies break to smaller scale eddies, and this break-down or “cascade” of energy continues until viscosity dissipates the energy of small scales. For small scales, first flows with Pr = 1 are considered, for which the 16 1.3. Length scales viscous dissipation and scalar diffusion occur at the same rate. For such cases the smallest length scale of motion is the Kolmogorov length scale, LK , first proposed by Kolmogorov in 1941. Below this scale, kinetic energy of eddies is severely suppressed by viscosity. We refer to the smallest velocity scale as uK . These scales are also called the “inner scales” in the literature (Tennekes & Lumley, 1972). It is conventional to relate LK to the rate of viscous dissipation of kinetic energy, ε, which is defined as ε = 2ν〈SijSij〉, (1.10) where S is the strain tensor, the indicial notation on the right hand side indicates summation of all terms resulted from the internal product of the two strain tensors, and 〈〉 denotes volume averaging. At the scale LK , sheets of strained fluid is formed over which dissipation (or diffusion) is active (Buch & Dahm, 1996). For such a structure the main component of the strain tensor is the inner strain rate γK = uK/LK , and ε ≈ νu2K/L2K . It is also known that at the dissipation scale Re = LKuK/ν ≈ 1. Hence, combining the latter two, the Kolmogorov length scale is estimated to be LK ≈ ( ν3 ε )1/4 . (1.11) Now we turn to flows with Pr > 1. For these flows the diffusivity of the scalar field is smaller than the viscosity in the vorticity field, and so the finest length scale in the density field is smaller than LK . The finest length scale in the density field is called the Batchelor scale, LB, after Batchelor (1959). Similar to the Pr = 1 flows, for Pr > 1 the molecular diffusion at the smallest scales takes place at sheet-like structures in the scalar field (Batchelor, 1959; Buch & Dahm, 1998). The thickness, L, of a diffusing sheet under the strain rate γ is estimated by (Davidson, 2005) L = (κ/γ)1/2. (1.12) 17 1.3. Length scales Hence, LB = (κ/γK) 1/2, as the smallest scale strain rate the sheets in the scalar field feel is γK. Combining this with the estimation for ε above, LB = (κ 2ν/ε)1/4 = LKPr −1/2. (1.13) This is the length scale below which any density gradient is smeared out by the diffusion. As seen from the equation 1.12, estimating the inner strain rate plays a key role in estimating the smallest scale the motion. Also equation 1.12 is based on the assumption that the strain rate has acted on the sheet long enough for the thickness to reach an equilibrium. In turbulent flows however the time scale of the smallest motion are sometimes too short to satisfy this condition. We discuss an alternative way of estimating γ in equation 1.12 that also ensures the equilibrium state in chapter 4 of the thesis. Another small length scale of interest in a turbulent flow field is the Taylor length scale, LT , first introduced by Taylor in 1921. This is an intermediate scale between the Kolmogorov and the large scale. The Taylor scale is large enough to be influenced by an outer-scale strain rate, ∆U/LT . Assuming an isotropic turbulence, substitution in equation 1.10 gives (Tennekes & Lumley, 1972) ε = 15ν ( ∆U LT )2 . (1.14) The coefficient 15 appears because at Taylor scale all terms of the strain ten- sor contribute to the deformation of the fluid structure, unlike a sheet-like structure at the smallest scales with only diagonal non-zero components in the strain tensor. The Taylor scale is large enough that it is not affected by the viscous dissipation. So in that sense it is a useful measure of the largest scale in the dissipation range. However, the location of the Taylor scale on the spectrum of scales between LK and δ has not been universal in different studies (Lewis & Swinney, 1999; Su & Clemens, 2003; Ishiharaa et al., 2005). Small scales can be related to the large scale, δ, given that the rate of 18 1.4. Overview destruction of energy at large scales should be equal to the rate of viscous dissipation of kinetic energy at small scales (Davidson, 2005), or ε ≈ ∆U 3 δ . (1.15) Substituting ε from equations 1.11, 1.13 and 1.14 gives LK ∼ Re−3/4δ (1.16) LB ∼ Re−3/4Pr−1/2δ (1.17) LT ∼ √ 15Re−1/2δ. (1.18) To accurately capture mixing down to the molecular level it is important to resolve the numerical domain to scales close to LK or LB. The relevance of these small scales in shear layer mixing is discussed in chapter 4 of the thesis. 1.4 Overview The objective of this thesis is to use numerical simulations to explore different aspects of the generation, evolution and turbulent breakdown of KH instabil- ities that are not yet well understood. These aspects include the evolution of large and small-scale structures of the flow and mixing characteristics for a wide range of Re0 and Pr. In stratified shear flows in the ocean and atmosphere the Reynolds number typically varies between 103 . Re . 108. The stratification usually results from gradients in temperature or salinity. The Prandtl numbers corresponding to these stratification are: Pr ≈ 1 for thermal diffusion in atmosphere, Pr ≈ 9 for thermal diffusion in water and Pr ≈ 700 for salt water diffusion in water (usually referred to as Schmidt number in this case). While the variation of Re0 and Pr has important consequences for the evolution of instabilities and the induces mixing, these effects are still not well studied. Especially mixing 19 1.4. Overview events at high Pr has not received much attention in previous studies. This is because of the extremely challenging requirements for numerical simulations and field measurements at high Re0 and Pr. In chapter 2, the effect of Re0 and Pr on the evolution of the primary KH billow is studied. The semi-analytical model of Corcos & Sherman (1976) is used to examine the vorticity structure of the evolving billow for a range of Re0 and Pr representing real oceanic flows, as an extension to previous work of Patnaik et al. (1976); Staquet (1995); Smyth (2003). This chapter also examines the small-scale length scales in turbulent KH flow field and the possibility of using a newly defined, Corcos-Sherman length scale as a guide for numerical simulations. Chapter 3 describes the mixing induced through the life-cycle of KH insta- bilities for a range of Reynolds numbers. As Re0 increases there is a transition in the characteristics of the mixing, associated with the enhanced entrain- ment and mixing by the growth turbulent motions (Dimotakis, 2000). This transition in mixing is examined in detail with the goal to explain how the more energetic flow field at higher Re0 modifies the entrainment and therefore mixing. Mixing characteristics are quantified and the transitional range of Reynolds numbers are compared to those observed by Breidenthal (1981) and Koochesfahani & Dimotakis (1986) in the laboratory. The effect of Prandtl number on the evolution of KH billows and their mixing behavior is studied in chapter 4. Slower turbulent mixing due to lower molecular diffusivity was observed in the experiments of Turner (1968). The energizing effect of higher Pr on turbulent motions was indicated in numerical simulations of Cortesi et al. (1998) and Staquet (2000). The results of chapter 4 shows how these effects influence the mixing behavior in different ways. 20 Chapter 2 The evolution of large and small-scale structures in Kelvin-Helmholtz instabilities 1 2.1 Introduction Sheared, density-stratified, flows frequently present in the ocean and atmo- sphere, can become unstable to Kelvin-Helmholtz (KH) type instabilities under conditions predicted by linear stability analysis (Drazin & Reid, 1981). The mixing induced by the nonlinear evolution and turbulent breakdown of these instabilities is a major source of mixing in the ocean and atmosphere (Woods, 1968; DeSilva et al., 1996; Smyth et al., 2001; Moum et al., 2003). The evo- lution of the primary instability, and therefore the transition to turbulence depend on the strength of the shear and stratification, as well as properties of the initial flow and the stratifying agent (Klaassen & Peltier, 1985a, 1991; Lawrence et al., 1991; Caulfield & Peltier, 2000; Peltier & Caulfield, 2003; Sta- quet, 2000; Smyth & Winters, 2003; Smyth et al., 2005). During its life-cycle a KH instability experiences two-dimensional growth of the primary billow to a saturation point followed by two-dimensional post- saturation evolution, where in the absence of vortex pairing the billow merely exchanges energy with the background flow without growing, subsequently 1A version of this chapter has been submitted for publication. M. Rahmani, B. R. Seymour and G. A. Lawrence (2011) The evolution of large and small-scale structures in Kelvin-Helmholtz instabilities. 21 2.1. Introduction there is a three-dimensional transition to turbulence and eventually turbulent breaking of the billow (Thorpe, 1973; Bernal & Roshko, 1986; Caulfield & Peltier, 2000; Rahmani et al., 2011a). The primary two-dimensional insta- bility starts by rolling up the density interface and forming a series of two- dimensional elliptical core vortices, also referred to as billows, connected by thin braids containing the remnant of core vorticity, with hyperbolic stagnation points in the middle of the braids (Brown & Roshko, 1974; Corcos & Sherman, 1976). Although mixing occurs mainly during the three-dimensional turbulent stages, understanding the growth of the primary instability is important since it is the two-dimensional pre-turbulent flow that determines characteristics of the transition and consequently the total mixing induced (Corcos & Lin, 1984; Smyth & Peltier, 1991; Schowalter et al., 1994; Caulfield & Peltier, 2000; Peltier & Caulfield, 2003; Smyth & Winters, 2003). Specifically, braid stag- nation points can be the origin of three-dimensional instabilities that lead to turbulence (Smyth & Peltier, 1991; Caulfield & Peltier, 2000). Here we study the two-dimensional evolution of the instability up to saturation in more depth to gain insight into the physics of the primary instability. Our investigation includes KH instabilities with Prandtl numbers and Reynolds numbers varying over three decades. The initial growth of the instability can be described by linear theory (Hazel, 1972), but as the instability grows the evolution becomes highly non- linear. All analytical models of the non-linear stage inevitably make use of simplifying assumptions and semi-analytical techniques (e.g. see del Castillo- Negrete, 2000). The semi-analytical model of Corcos & Sherman (1976), CS76 hereafter, describes the evolution process up to saturation using vorticity dy- namics of the rolling-up billow. The model benefits from its simplicity and computational efficiency as only the braid region is examined. Across the braid, properties are described by similarity solutions. Here we examine the large-scale vorticity structure of the core and braid, as well as the small-scale structures inside the core and on the braid. The former are compared to the 22 2.1. Introduction time evolution model of CS76, and the latter to the similarity solution of CS76 for the braid. The predictions of CS76 for the saturation state are in good agreement with the low Reynolds number, salt stratified (high Prandtl number) labora- tory experiments of Koop & Browand (1979) and the numerical simulations of Staquet (1995) for high Reynolds numbers and thermal atmospheric stratifica- tion (low Prandtl number). Smyth (2003) validated and extended the similar- ity solutions of CS76 for the saturation state and applied them to predicting the occurrence of secondary instabilities on the braid. Our direct numerical simulations (DNS) of KH billows extend these studies to a wider range of Reynolds and Prandtl numbers and examine further details of the structure of the vorticity field. The evolution of KH instabilities has been studied using DNS (e.g. Scinocca, 1995; Staquet, 1995, 2000; Caulfield & Peltier, 2000; Peltier & Caulfield, 2003; Smyth & Winters, 2003; Inoue & Smyth, 2009), where the numerical domain is resolved to the finest scale of motion. Although this technique provides the most accurate simulation of all flow fields, its main drawback is the require- ment of large computational resources. The minimum length scale is estimated by the Batchelor length scale (Batchelor, 1959), and hypothetically is small- est at the maximum intensity of turbulence. Here we discuss the role of the hyperbolic stagnation point in determining the pre-turbulent minimum length scale that locally occurs on the braid, based on which we define the Corcos- Sherman length scale. The similarity solution of the braid and the large-scale structure prediction of CS76 are used to predict this length scale. While the scale defined by Batchelor (1959) gives a bulk measure of the smallest scale below which scalar gradients are suppressed by diffusion, the minimum braid thickness prediction approach seeks the local fine scales. Finest scale predic- tion is valuable in determining the required resolution for DNS. This chapter is outlined as follows. A background on KH instabilities and the model of CS76 is given in section 2, while the DNS method is described 23 2.2. Background in section 3. The pre-saturation structure of the billow is discussed in sec- tion 4. The large-scale temporal evolution of the billow and the predictions of the CS76 model are discussed in section 5. In section 6, we discuss the similarity solutions for the braid steady state at saturation. Minimum length scales of the flow field that should be considered in DNS, and the relevance of the Corcos-Sherman length scale are discussed in section 7. Conclusions are presented in section 8. 2.2 Background 2.2.1 Governing equations Using the Boussinesq approximation the Navier-Stokes equations in a Carte- sian coordinate system of (x, y, z) take the form Du Dt = −∇P ρ0 − ρ ρ0 gk̂+ ν∇2u, (2.1) where u = (u, v, w) denotes the three-dimensional velocity field, P the pressure field, and ρ the density field. ρ0, g, and ν are the reference density, gravitational acceleration and viscosity, respectively. D/Dt denotes a material derivative and k̂ the unit vertical vector. Incompressibility of the fluid reduces the continuity equation to ∇.u = 0, (2.2) while the evolution of the density field is described by the advection-diffusion equation Dρ Dt = κ∇2ρ, (2.3) where κ is the molecular diffusivity of the stratifying agent. Equations 2.1 to 3.12 are solved using DNS as described in section 2.3. 24 2.2. Background 2.2.2 Relevant dimensionless numbers The initial flow has horizontal velocity and density profiles, with a vorticity interface thickness of h0 and density interface thickness of δ0. The horizontal velocity and the density differences are denoted by ∆U and ∆ρ, with initial profiles expressed as U(z) = ∆U 2 tanh( 2z δ0 ), and ρ̄(z) = −∆ρ 2 tanh( 2z h0 ). (2.4) The flow is governed by four dimensionless numbers Re0 = ∆Uδ0 ν , J = ∆ρgδ0 ρ0∆U2 , P r = ν κ , and R = δ0 h0 , (2.5) the initial Reynolds number, bulk Richardson number, Prandtl number and scale ratio. By scaling the thickness of the density and vorticity interfaces, which grow in time due to the diffusion, it can be shown that the Prandtl number is related to the scale ratio, R, through R ∽ (Pr)1/2. 2.2.3 Semi-analytical model of Corcos & Sherman (1976) CS76 describe the evolution of the instability using vortex dynamics, as follows. The background flow with a strip of vorticity concentration evolves to periodic regions of vorticity concentration connected by thin tilted vorticity containing strips, see figure 2.1. The former are referred to as cores, and the latter as braids. Stagnation points are formed in the middle of a braid separating two cores and in the center of the cores. Concentration of vorticity in the core evokes a strain rate on the braid that advects more vorticity to the core. The Boussinesq vorticity equation for the two-dimensional scalar vorticity, 25 2.2. Background Ω = ∂u/∂z − ∂w/∂x, is expressed as DΩ Dt = ν∇2Ω + g ρ0 ∂ρ ∂x , (2.6) where the first term on the right hand side is the viscous diffusion and the second term the baroclinic generation of vorticity. In a domain with periodic boundary conditions and sufficiently far top and bottom boundaries (so that vorticity fluxes vanish), the total circulation, Γ = ∫ A ΩdA , with A being the domain area, is conserved and remains equal to its initial value ∆ULx, with Lx being the wavelength of the instability. The two mechanisms of baroclinic vorticity generation and vorticity advection act in competition and determine how Γ is distributed between the braid circulation Γb and core circulation Γc, while Γ = Γc + Γb remains constant. As Γc rises the core gets stronger, and the braid experiences a higher strain rate. By examining the rate of advection of vorticity from the braids, the model of CS76 predicts the evolution of the distribution of the total circulation and therefore the characteristics of the instability. A locally orthogonal moving frame of reference is used with σ denoting the direction along the braid and η the direction of the normal, with the stagnation point set as the origin, see figure 2.1. The height of the billow, H , is defined as the vertical distance between the highest and lowest point of the braid. The η dependence of equation 2.6 is eliminated by describing the braid in terms of the total local shear across the interface defined as S = ∫ Ωdη. By integrating the vorticity equation 2.6 across the interface, making the assump- tion that ∂/∂σ ≪ ∂/∂η and noting that the ∂/∂η terms vanish far away from the braid, the equation for the time evolution of S when moving along the braid is dS dt = −γS + g′sinθ, (2.7) where the material derivative d/dt here is ∂/∂t + uσ∂/∂σ, with uσ being the braid centerline velocity along σ. γ(t) is the strain rate, with its com- 26 2.2. Background σ η H L xcore center Braid Figure 2.1: Configuration of braid and core regions and braid coordinates. The color scheme shows the vorticity distribution, with the blue shading showing a higher concentration of vorticity. The braid centerline is shown by the black dashed line and the density interface by the green line. The blue dashed line shows the central tangential line. pressive component aligned with η and expansive component aligned with σ. g ′ = ∆ρg/ρ0 denotes the reduced gravitational acceleration, and θ the local inclination angle of the braid. This equation states when travelling at the speed uσ with a particle along the σ axis on the braid center, shear is lost due to the action of the strain field and gained baroclinically because of the local inclination of the braid. It should be noted that equation 2.7 is not directly dependent on Re0 or Pr. However, γ and θ, which are both measures of the core strength, vary depending on Re0 and Pr. The billow reaches a saturation or self-limiting state when the core has consumed all the available vorticity in the braids and dS/dt = 0, also implying ∂S/∂σ = 0. At this stage the core stops gaining more vorticity. Consequently, the core strength and therefore the height by which the core rolls up the density interface, H , do not grow any further. For stratified flows positive vorticity is constantly generated on the braids due to the tilting of the interface. So in this case saturation is attained when the baroclinic generation of vorticity is in balance with its advection to the core. The model of CS76 proceeds to the saturation steady state from an initial state of weakly perturbed vorticity strip, taking the steps that follow. Know- ing S at any time, the braid circulation can be found as: Γb = ∫ Lx/2 −Lx/2 Sdσ. 27 2.3. Direct numerical simulation Equation 2.7 along with the braid streamline equation, dσ/dt = uσ are solved numerically as an initial value problem using Lagrangian markers. However, the velocity field is needed at any time step, and this is generated using stream- functions introduced by Stuart (1967). These generate a close approximation to the true vorticity structure and can be quantified in terms of the total core circulation. For more details of this substitute flow the reader is referred to CS76. The main assumptions of the CS76 model are that no vorticity re-circulates from the core to the braid and that baroclinically generated vorticity can be treated as a passive scalar. By comparing the analytical predictions to the numerical simulations of Patnaik et al. (1976), CS76 concluded that the first assumption requires the waves to be sufficiently long. Long waves were defined as: α = kδ0/2 < 0.2, with k = 2pi/Lx, which is equivalent to: Lx/δ0 > 15.7. They showed that by satisfying this condition, the predictions of the model remain accurate up to the saturation point even for small Reynolds numbers. The above solution gives the temporal characteristics of the braid and defines the large-scale billow structure. Across the braid similarity solutions are used for Ω and ρ, which give the braid structure dependence on η. The similarity solutions are discussed in more details in section 6. 2.3 Direct numerical simulation The equations of motion are solved using the DNS spectral model described by Winters et al. (2004). To extend simulations to three dimensions a mod- ified version of this model described by Smyth et al. (2005) is used. In this model the resolution of the density field is two times finer than that of the velocity field. Two-dimensional simulations at different Reynolds and Prandtl numbers were carried out, see table 4.1. The rationale for the choices of non- dimensional parameters in these simulations are as follows. In all simulations 1 to 10 J = 0.035, the same J as in simulations of Patnaik et al. (1976). Firstly, 28 2.3. Direct numerical simulation two of the numerical simulations of Patnaik et al. (1976) used by CS76 for comparison to their model are duplicated. These simulations are at Re0 = 200 and Re0 = 400 with Pr = 0.72 (simulations 1 and 2). Then, a sequence of Re0 = 800, 1600, 3200, 10,000, 40,000 and 100,000 with Pr = 0.72 is chosen to extend the results to higher Reynolds numbers (simulations 3 to 8). Next, for a constant value of Re0 = 400, Prandtl numbers of 0.72 (atmospheric), 9 (thermal) and 700 (salt water) are considered (simulations 2, 9 and 10, respec- tively). Simulation 11 is performed to simulate the occurrence of secondary instabilities on the braid. To capture secondary instabilities we had to increase Re0 to 240,000 and J to 0.2 for this simulation. Linear theory predicts KH in- stability in all these cases. In all two-dimensional simulations one wavelength of the fastest growing instability is simulated and the vertical height of the domain is Lz = 1.25Lx. Therefore, we are neglecting pairing (as captured by Winant & Browand, 1974; Smyth et al., 2005; Carpenter et al., 2010b) to particularly focus on the evolution of one billow. In addition to the two-dimensional simulations, we are examining a series of three-dimensional simulations (simulations 12 to 21 in table 2.2) to check the applicability of the Corcos-Sherman simulations to three-dimensional tur- bulent flow fields. These three-dimensional simulations have slightly different configurations from the two-dimensional ones as they were performed indepen- dently and borrowed from another study (Rahmani et al., 2011a, submitted to JFM). In these simulations J = 0.03, Lz = 9δ0 and the span-wise width of the domain is Ly = Lx/2. For Pr = 9 the Reynolds number ranges between 400 and 1200. For Pr = 1 a wider range of Reynolds numbers of 400 to 4800 are simulated. The required resolution for the DNS of stratified flows with Prandtl num- bers greater than 1 is estimated from the Batchelor length scale, LB = (νκ 2/ε)1/4, (2.8) with ε being the viscous energy dissipation. For Pr 6 1, the Kolmogorov 29 2.4. Pre-saturation structure of the instability length scale LK = (ν 3/ε)1/4, (2.9) is used to estimate the required resolution, so that LK = LBPr 1/2. Grid spacings of 2.5LB for Pr > 1 and 4LK for Pr 6 1 were found to be sufficient (Smyth & Winters, 2003; Moin & Mahesh, 1998). In all our simulations the domain is resolved to at least 2.5LB or 2.5LK throughout the simulation. Periodic horizontal and free slip vertical boundary conditions are applied. Basic velocity and density profiles in equations 4.1 are perturbed using the eigenfunctions from the numerical solution of the linear stability equation. The magnitude of this perturbation is large enough to cause a maximum dis- placement of 0.1δ0 of the interface, and small enough for the perturbation to be linear. For three-dimensional simulations, a random noise with the mag- nitude of ±0.05∆U is added to the velocity field to initiate three-dimensional secondary instabilities. 2.4 Pre-saturation structure of the instability In this section we discuss the vorticity structure of the braid and core from DNS. The focus is on the pre-saturation behavior of the billow. The satura- tion state from the CS76 model was previously defined in section 2.3 using equation 2.7. In DNS the core and braid circulation reach their saturation values asymptotically, so defining saturation based on circulation introduces a level of arbitrariness. We found dH/dt = 0, which is a consequence of sat- uration, a more robust definition. H is numerically found as the difference between the highest and lowest points of the density interface, see figure 2.1. The dimensionless saturation times based on this definition, t∗s, are presented in table 4.1, where time is non-dimensionalized by t∗ = t∆U/δ0. 30 2.4. Pre-saturation structure of the instability Simulation Re0 Pr R Lx/δ0 Nx h(ts)/2δ0 t ∗ u t ∗ s ns 1 200 0.72 0.85 7.5 256 0.40 16 36 0.9 2 400 0.72 0.85 7.3 256 0.24 14 35 1.3 3 800 0.72 0.85 7.1 384 0.12 12 33 1.5 4 1600 0.72 0.85 7.1 768 0.085 10 32 1.6 5 3200 0.72 0.85 7.1 1200 0.063 10 32 1.7 6 10,000 0.72 0.85 7.1 768 0.035 9 32 1.8 7 40,000 0.72 0.85 7.1 2400 0.022 8 32 1.9 8 100,000 0.72 0.85 7.1 3200 0.015 8 32 1.9 9 400 9 3 8.5 512 0.048 14 42 1.7 10 400 700 26.5 10.1 2000 0.0066 19 46 1.6 11 240,000 1 1 6.7 3200 - - - - Table 2.1: Description of the two-dimensional simulations carried out. In all these simulations Lz = 1.25Lx and J = 0.035, except for simulation 11, for which J = 0.2. Also, the saturation half braid thickness divided by the initial vorticity interface thickness, h(ts)/2δ0, the dimensionless time the billow first becomes statically unstable, t∗u, the dimensionless saturation time, t ∗ s, and the number of rolls at saturation, ns, are presented for simulations 1 to 10. 2.4.1 Dependence on Re0, Pr = 0.72 The vorticity structure of the billow at different non-dimensional times along with velocity vectors and density interface are shown in figure 2.2 for differ- ent Re0. As explained by Batchelor (1972) the perturbation of the uniform vorticity concentration in the initial shear flow creates a velocity field that accumulates vorticity in the core by advecting vorticity away from the braid centers (see the t∗ = 10 panels). This accumulation of vorticity in the core further helps the roll-up of the instability and accelerates the advection of vor- ticity from the braid center to the core. A thin braid and core of vorticity have clearly formed by t∗ = 20. As seen from the velocity vectors, this is accom- panied by the formation of braid hyperbolic stagnation points (at mid-depth on the left and right side of the domain). At t∗ = 30 the billow is close to its saturation state, where the braid is thinnest and the core has maximum 31 2.4. Pre-saturation structure of the instability Simulation Re0 Pr R Lx/δ0 Nx ×Ny ×Nz h(ts)/2δ0 Pr = 9 12 400 9 3 8.2 384×192×384 0.0445 13 600 9 3 8.1 384×192×384 0.0365 14 900 9 3 8.0 480×240×480 0.0296 15 1200 9 3 8.0 512×256×512 0.0257 Pr = 1 16 400 1 1 7.4 256×128×320 0.0763 17 600 1 1 7.3 256×128×320 0.0618 18 900 1 1 7.2 256×128×320 0.0501 19 1200 1 1 7.2 256×128×320 0.0434 20 2400 1 1 7.1 384×192×512 0.0306 21 4800 1 1 7.1 512×256×640 0.0217 Table 2.2: Description of the three-dimensional simulations carried out. In all these simulations J = 0.03, Lz = 9δ0 and Ly = Lx/2. The presented mesh sizes are for the density field, the velocity field mesh sizes are half of these. Also, the saturation half braid thickness divided by the initial vorticity interface thickness, h(ts)/2δ0, are presented. strength. At t∗ = 40, the billow height and core vorticity concentration relax back after passing the saturation state. At this stage the flow starts to become too complex to be described by the core and braid structure. As Re0 increases in figure 2.2, less viscous diffusion of vorticity occurs dur- ing the evolution of the core and braids. Hence, at all stages more concentrated core and braid vorticity prevails for higher Re0. The stronger core at higher Re0, due to reduced diffusion, invokes a higher γ on the braid and hastens the growth of the core. As a result the saturation time decreases from t∗s = 36 at Re0 = 200 to t ∗ s = 32 at Re0 = 1600, and remains almost the same for higher Reynolds numbers. The excessive viscous diffusion at low Re0 results in an early expansion of the core to the braid (e.g. see panel b in figure 2.2), involving re-circulation of vorticity from the core to the braid. Due to this vorticity re-circulation the flow is highly dependent on Re0 for low Reynolds 32 2.4. Pre-saturation structure of the instability numbers. However, for sufficiently high Re0 clearly distinct core and braid vorticity structures are retained prior to saturation times and no vorticity re-circulation occurs. While the large-scale core and braid structures reach an almost Re0 inde- pendent regime for Re0 > 1600, the local concentration of vorticity in certain regions increase as Re0 rises. On the braid, while vorticity is advected along the σ axis (see figure 2.1) and away from the the stagnation point, the vorticity concentration across the braid (along the η axis in figure 2.1) is determined by the competitive action of the compressive strain rate and diffusion. As Re0 increases the diffusion becomes less dominant and the braid centerline vortic- ity (i.e. Ω at η =0) rises more. The density field remains less diffuse during the evolution of the instability for higher Re0, resulting in more concentrated generation of baroclinic vorticity due to sharper density gradients. This fur- ther helps the formation of less diffuse core and braid vorticity structures. On the braid, positive baroclinic generation of vorticity intensifies at higher Re0. Inside the core, baroclinic effects generate positive vorticity at interfaces with ∂ρ/∂x > 0 and negative vorticity at interfaces with ∂ρ/∂x < 0. Negative gen- eration of vorticity counter-balances the core’s positive vorticity on the latter interfaces and appears as red patches in panels h, l, p and t of figure 2.2. However, no negative vorticity is formed prior to saturation as all these panels correspond to post-saturation times. The core structure and its baroclinically generated vorticity is determined by the progress of the rolling inside the billow. While the roll-up proceeds faster as Re0 increases, the evolution of the rolling layer is fairly similar for the high Re0 billows. The number of rolls, n, can be defined as the number of times that the line tangential to the density interface in the center of the core, see figure 2.1, winds around the center of the core. At n = 1/4 the density interface becomes statically unstable for the first time. We denote this time by t∗u. The number of rolls at saturation, ns, and t ∗ u are presented in table 4.1. Both ns and t ∗ u asymptote to a constant number at high Re0, indicating an 33 2.4. Pre-saturation structure of the instability inviscid regime. The number of rolls is an indication of the evolution of the primary billow and strength of the core vorticity. A higher n indicates a more density structured core, more baroclinically generated vorticity inside the core and greater potential for mixing. 2.4.2 Dependence on Pr, Re0 = 400 The vorticity structure of the billows at different non-dimensional times are presented in figure 2.3 for different Pr. Viscous diffusion of the vorticity occurs at the same rate for all Pr, but the density field diffuses less during the evolution of the instability as Pr increases. As discussed in section 4.1, Re0 = 400 falls in the viscous diffusion dominant range of Re0. However, increasing Pr has similar influences as increasing Re0 on both the large-scale and local concentration of vorticity: less re-circulation of vorticity from the core to the braid and more pronounced baroclinic generation of vorticity. At high Pr the billow grows higher (compare panels d, h and l in figure 2.3) and the braid reaches steeper slopes, indicating a stronger core. The higher concentration of baroclinic generation of vorticity at higher Pr creates a more structured core, rather than a patch of diffusing vorticity. Regions of negative vorticity appear within the core at earlier times and grow stronger as Pr increases. The accumulation of negative vorticity in the inner and outer edges of the outermost rolling layer wraps the core positive vorticity (panels h and l of figure 2.3). The significance of this feature is that it limits the diffusion of core vorticity. Hence, in addition to varying the local small-scale structure of the vorticity field, the baroclinic generation of vorticity can also influence the large-scale structure of the vorticity field by inhibiting the vorticity re- circulation. Less vorticity re-circulation at higher Pr is partly because of a less diffuse core vorticity structure and partly due to longer wavelengths at higher Pr. By increasing Pr from 0.72 to 700, saturation time is delayed from t∗s = 36 to t∗s = 46, and the time the billow first becomes unstable from t ∗ u = 34 2.4. Pre-saturation structure of the instability a b c d e f g h i j k l m n o p q r s t Figure 2.2: Vorticity structure at t∗ = 10, 20, 30 and 40 from left to right for Pr = 0.72, and different Reynolds numbers : Re0 = 200, 800, 1600, 3200 and 10,000 from top to bottom. The color shading represents the non-dimensional vorticity. In each panel the aspect ratio is preserved. The solid line shows the density interface. The last panel shown for each Reynolds number (i.e. t∗ = 40) is after saturation. 35 2.5. Temporal evolution of the billow 14 to t∗u = 19. This is attributed to longer wavelengths at higher Pr that reduce the strain rate on the braid and therefore slow down the vorticity advection. Despite the slower growth, the core remains strong enough to advect the baroclinically generated vorticity on the braid and stronger cores at saturation are obtained for higher Pr. The rolling evolution does not have a straightforward dependence on Pr. The slower evolution of the roll-up at higher Pr is due to the slower growth of the instability. However, ns shows less dependence on either Re0 or Pr when vorticity re-circulation is small. The DNS shows that the process of evolution of vorticity from the braids to the cores is consistent with, but more complicated than, the simplified model of CS76. Viscous diffusion especially in low Re0 instabilities makes the no re-circulation assumption not fully satisfied so that the distinction between the core and the braid is rather arbitrary. Baroclinic generation of vorticity in high Pr instabilities results in a more complicated vorticity structure than the simple Stuart’s vortex used in the CS76 model. We will concede these differences when comparing the predictions of the CS76 model to the DNS results in section 5. 2.5 Temporal evolution of the billow In this section we quantify the large-scale evolution of the billow using the concept of core and braid structure discussed in the previous section. We compare the DNS results to the predictions of the model of CS76. The CS76 model starts from an initial perturbed condition where a great portion of vorticity is contained within the braid, and evolves by predicting the advection and baroclinic generation of vorticity on the braid using Lagrangian markers. In DNS, the numerical braid is found as the loci of maximum vorticity passing through the stagnation point in the left middle of the domain (the DNS domain shown in figures 2.2 and 2.3). As the CS76 model and DNS are initially perturbed differently, the model results are normally time shifted from those 36 2.5. Temporal evolution of the billow a b c d e f g h i j k l Figure 2.3: Vorticity structure at t∗ = 10, 20, 30 and 40 from left to right for Re0 = 400, and different Prandtl numbers : Pr = 0.72, 9 and 700 from top to bottom. See figure 2.2’s caption for the rest of details. The last panel shown for each Prandtl number (i.e. t∗ = 40) is after saturation for Pr = 0.72 and before saturation for Pr = 9 and Pr = 700. 37 2.5. Temporal evolution of the billow of DNS. To match these two the CS76 model results were shifted in time to get the same initial total core (or braid) circulation. We first examine the DNS vorticity depletion during the instability evolution in section 5.1, then the height of the billow in section 5.2. Lastly, we discuss how the variation in Re0 and Pr influences the saturation state in section 5.3. 2.5.1 Vorticity depletion Time variation of the ratio of braid circulation to total circulation is shown in figure 2.4. Braid circulation, initially close to total circulation, depletes to the core until saturation. A small portion of the initial circulation remains in the braid at saturation as baroclinic vorticity is constantly generated on the braid. The vorticity depletion process is well predicted by the CS76 model, except at Re0 = 200. The DNS braid circulation evolution becomes closer to the predictions of the CS76 model as Re0 and Pr increase. In DNS, as we discussed in section 4.1, the core expands to the braids due to its viscous diffusion and growth. At low Re0 and Pr this expansion results in the re-circulation of positive vorticity to the braid before saturation. Therefore, the DNS braid circulation does not reach its expected asymptotic value from the CS76 model at saturation and starts to rise after saturation. At higher Re0 and Pr there is less vorticity re-circulation and it occurs later. So, the saturation state is reached asymptotically and at values closer to the predictions of CS76. In these cases the growth of the core brings negative baroclinically generated vorticity to the braid (see panel t in figure 2.2 for instance). This accounts for the slightly decreasing trend of Γb close to and after saturation in panels (c) and (d) of figure 2.4. In addition to the vorticity re-circulation, the approximation of the flow field by the Stuart solution also contributes to the difference between the DNS and CS76 model braid circulation. Velocity and strain rate on the braid re- motely induced by core vorticity are the key parameters that determine the transport of vorticity to the core. Similar to the comparisons by CS76, by 38 2.5. Temporal evolution of the billow (b) Pr = 0.72, Re0=400 0 0.5 1 (a) Pr = 0.72, Re0=200 Γ b /Γ 0 10 20 30 40 50 (d) Pr = 700, Re0=400 t* 0 10 20 30 40 50 0 0.5 1 (c) Pr = 0.72, Re0=10,000 t* Γ b /Γ Figure 2.4: Comparison of DNS (dashed) and CS76 model (solid) ratio of braid to total circulation for increasing Reynolds numbers (a) and (c) and increasing Prandtl numbers (b) and (d). The vertical dashed lines in each panel show saturation. matching the braid circulation between the real and approximated flow the CS76 model slightly under-estimated the braid centerline velocity, but pro- vided reasonable approximations for the strain rate. After comparing the real and approximated flow fields, we found that Stuart’s approximation improves as the instability evolves and more vorticity accumulates in the core. In gen- eral the error introduced by the flow field approximations is small compared to the vorticity re-circulation effects. 2.5.2 Billow height The DNS time evolution of billow heights are compared to the predictions of CS76 in figure 2.5. The billow height is a proper measure of the largest scale in the two-dimensional flow. In DNS, the billow height grows until saturation as more vorticity accumulates inside the core, then reduces as the core strength diminishes after saturation. The DNS heights are in agreement with those of the original simulations of Patnaik et al. (1976). As expected from the earlier 39 2.5. Temporal evolution of the billow (b) Pr=0.72, Re0=400 0 0.1 0.2 (a) Pr=0.72, Re0=200 H /2 L x 0 10 20 30 40 50 t* (d) Pr=700, Re0=400 0 10 20 30 40 50 0 0.1 0.2 (c) Pr=0.72, Re0=10,000 t* H /2 L x Figure 2.5: Comparison of DNS (dashed) and CS76 model (solid) half billow heights normalized by the wavelength for (a) and (c): increasing Reynolds numbers, (b) and (d): increasing Prandtl numbers. The vertical dashed lines in each panel show the saturation. Black triangles in panels (a) and (b) show the original numerical simulations of Patnaik et al. (1976) used by CS76 for comparison to their model. Γb discussion, the CS76 model predictions for H improve for increasing Re0 and Pr. The saturation heights are very well predicted by the CS76 model for high Re0 and Pr, even though none of our simulations comply with CS76’s condition of long waves (Lx/δ0 > 15.7). In none of the simulations is the CS76 model capable of predicting accurate heights at all times. The real heights are almost always over-estimated by the CS76 model, partly as a result of the mismatch between the DNS and CS76 model braid circulations in figure 2.4 caused by the vorticity re-circulation. However, in section 5.1 we showed that the re-circulation of vorticity is greatly overcome by increasing Re0 and Pr. Hence, the over-estimation of H by the CS76 model at times when Γb’s match, see figure 2.4, cannot be associated with the vorticity re-circulation. Instead, the difference stems from the approxi- mate braid configuration considered in Stuart’s vortices. This approximation improves as the instability evolves, and achieves great accuracy at saturation. 40 2.5. Temporal evolution of the billow 2.5.3 Saturation The self-limiting state of the instability is reached when the vorticity advec- tion and generation on the braid balance. At this state the vorticity advection towards the core also balances the more dominant negative generation of vor- ticity inside the core, and the large-scale strength of the core reaches a steady state. However, as seen in section 4, the local concentrations of vorticity con- tinue to vary. The saturation state is characterized by the CS76 model by predicting how much vorticity at this stage resides in the core and the braid. As seen in figures 2.4 and 2.5, this stage is reached asymptotically as a steady state in the CS76 model. At saturation times predicted by DNS, the CS76 model has just attained a state with small variations in Γb and H . Although the DNS values do not quite achieve a steady state, their time variations are small at saturation. The dimensionless billow height, braid shear and strain rate at saturation, from both DNS and the CS76 model, are presented in figure 2.6 for varying Re0 and Pr. For both DNS and the CS76 model, the shear and strain rates are averaged over x = 0 to x = Lx/8 on the braid to avoid the small numerical noise in DNS. Since non-dimensionalizing equation 2.7 eliminates its depen- dence on Re0 and Pr, none of the non-dimensional parameters shown in figure 2.6 are directly related to Re0 and Pr in the CS76 model. The CS76 model dependence on Re0 and Pr in this figure however comes from different dimen- sionless strain rates at different wavelengths. Higher strain rates generated at shorter wavelengths hasten the vorticity advection on the braid and result in stronger cores at saturation. A stronger core leaves less shear on the braid and generates a higher roll-up. This explains the trends of variation of H(ts) and S(tS) with Re0 and Pr in figure 2.6. In DNS, at low Re0 and Pr the domination of vorticity re-circulation limits the core strength at saturation and causes significant difference between the CS76 model and DNS. This difference decreases as Re0 and Pr increase. Both panels (a) and (c) of figure 2.6 suggest that stronger saturation cores are formed 41 2.6. Similarity solutions for higher Re0. For Re0 > 300, which is the start of the transition in mixing as found by Rahmani et al. (2011a) (although for slightly different Pr and J), there is considerably less dependence on Re0 for H , S and γ and improvement in the predictions by the CS76. This is the critical Reynolds number beyond which the flow overcomes the viscous forces and sustains three-dimensional motions. Here, we also observe the correspondence of this Reynolds number to a transition in the two-dimensional billow structure. For Re0 > 900, after the completion of the mixing transition (Rahmani et al., 2011a), there is very little dependence of H , S and γ on Re0. Hence, the transition from viscous to inviscid two-dimensional flows seen in figure 2.6 is similar to the transition in mixing in three-dimensional flows. The CS76 model predicts weaker cores at saturation for longer wavelengths at higher Pr. However, this prediction is suppressed in DNS, as seen in panels (b) and (d), by the counterbalancing effect of vorticity re-circulation. The similarity solutions that will be discussed in section 6, rely on S and γ predictions at saturation. Figure 2.6 shows that S can be predicted very well for high Re0 and Pr. Predictions of γ are accurate for all values of Re0 and Pr studies. 2.6 Similarity solutions In the previous section we discussed the large-scale properties of the billow, which revealed a Reynolds number independent behaviour for sufficiently high Re0. Local small-scale braid properties are described by the similarity solu- tions. Across the braid, the similarity forms for density and vorticity proposed by CS76 are ρ(η) = −∆ρG(η/hb(t)), and Ω(η) = 2S(t) δb(t) F (η/δb(t)), (2.10) where G and F are arbitrary functions, η denotes the axis normal to the braid, 42 2.6. Similarity solutions 0.1 0.2 0.3 H (t s )/2 L x (a) 0.1 0.2 0.3 S( t s) / ∆ u (c) 102 103 104 105 0.15 0.2 0.25 γ(t s) δ 0 / ∆ u Re0 (e) (b) (d) 100 101 102 103 Pr (f) Figure 2.6: Dimensionless billow height, braid shear, and strain rate at sat- uration for varying Reynolds numbers with Pr = 0.72 and varying Prandtl numbers with Re0 = 400. The solid lines indicate the model of CS76 and the triangles the DNS. S and γ are averaged over x = 0 to x = Lx/8. The vertical dashed lines show the range of the mixing transition Reynolds numbers found by Rahmani et al. (2011a). and hb and δb are measures of the density and vorticity interface thickness on the braid. The CS76 model assumes hb = Rδb, which remains reasonably valid throughout the braid evolution in DNS. Similarity solutions of the density field and thickness of the density interface (equation 3.12) are obtained with relative ease (see Corcos & Sherman, 1976): ρ(η) = −∆ρ 2 erf( pi1/2 2 η hb(t) ). (2.11) For sufficiently large time the diffusion of the interface balances the thinning 43 2.6. Similarity solutions by the strain field and hb(t) ∼ ( 2piκ γ(t) )1/2 as t→∞. (2.12) As the strain rate has reached its asymptotic value at saturation, the saturation density interface thickness, hb(ts), is described by equation 2.12. Equation 2.12 can be rearranged to give hb(ts) δ0 = ( 2pi PrRe0γ∗(ts) )1/2 , or hb(ts) h0 = ( 2pi Re0γ∗(ts) )1/2 , (2.13) where γ∗ = γδ0/∆U is the non-dimensional strain rate. The γ ∗ dependence on Re0 and Pr is nearly compensated by the Lx/δ0 dependence on Re0 and Pr and a linear proportionality of hb(ts)/2Lx to Re −1/2 0 and Pr −1/2 results (figure 2.7). As γ∗ revealed an almost Re0 independent behavior and a dependence on Pr through different wavelengths (figure 2.6), it is expected that ratio by which the asymptotic braid density interface reduces compared to its initial value, hb(ts)/h0, varies very closely with Re −1/2 0 and be weakly dependent on Pr. The asymptotic density thicknesses from DNS and the CS76 model, non- dimensionalized by the domain length, are compared in figure 2.7 for varying Re0 and Pr. Using the density similarity solution, equation 2.11, the DNS density thickness is found as ρ(η = hb/2) = −∆ρ erf(pi1/2/2)/2. The DNS and CS76 model thicknesses are averaged over x = 0 to x = Lx/8. Except for the two lowest Re0, the DNS thicknesses follow the CS76 model thickness lines closely in both panels. The same transition as in section 5.3. is seen in the behaviour of hb beyond the Reynolds number of the beginning of the mixing transition, Re0 = 300. In general, we found the density profiles a good fit to the similarity structure given by equation 2.11. This good fit combined with the accurate predictions of γ(ts) (see figure 2.6) result in close predictions of hb(ts) by the CS76 model. The differences at the two low Re0 correspond to 44 2.6. Similarity solutions 102 103 104 105 10−3 10−2 10−1 Re0 h b (t s )/2 L x (a) 100 101 102 103 Pr (b) Figure 2.7: Variation of the ratio of the asymptotic half density interface thick- ness at saturation to the wavelength with (a) Reynolds number for a constant Prandtl number of 0.72 and (b) Prandtl number for a constant Reynolds num- ber of 400. The solid lines indicate the model of CS76 and the triangles the DNS. h is averaged over x = 0 to x = Lx/8. The vertical dashed lines delineate the range of the mixing transition Reynolds numbers found by Rahmani et al. (2011a). the diffusive expansion of the core density to the braid. The similarity solution for vorticity takes a simple form only for Pr = 1, where δb(t) = hb(t): Ω(η) = S(t) δb(t) exp (−pi 4 ( 2η δb(t) )2 ) . (2.14) To extend this solution to higher Pr but only for the steady state at satura- tion, Smyth (2003) proposed a similarity form which in terms of dimensional parameters is expressed as Ω(η) = 2g ′ sin(θ) γ(ts)δb(ts)Pr F (η/δb(ts)). (2.15) 45 2.6. Similarity solutions Substituting this form in the vorticity equation gives d2F dξ2 + 2ξ Pr dF dξ = −e−ξ2 , ξ = pi1/2 η hb . (2.16) For Pr =1, and g ′ sin(θ) = γS at saturation, which can be seen from equation 2.7, this similarity form converts to the similarity form of CS76, and the above equation has the solution F = e−ξ 2 /2. For Pr 6= 1, F can be solved numerically and the equilibrium vorticity structure predicted. The DNS braid centerline vorticity, i.e. Ω(η = 0), at saturation is compared to the similarity solution vorticity in figure 2.8 for varying Re0 and Pr. The similarity solution vorticity is obtained by numerically solving equation 2.16 and substituting F in equation 2.15, using saturation values from the CS76 model for each parameter. Similar to the rest of comparisons in figures 2.6 and 2.7, both similarity and DNS vorticity are averaged over x = 0 to x = Lx/8. The correct trend of variation of vorticity with Re0 and Pr is predicted by the similarity solution, with less than 10% error except for Re0 = 200 and 100,000. At low Re0 and Pr despite large differences in predictions of S in figure 2.6 due to the re-circulation of vorticity, the error in braid centerline vorticity predictions in figure 2.8 remain relatively small. This indicates that at satuation the vorticity recirculation has caused the core vorticity to intrude into the braid vorticity structure, but left the near centerline braid vorticity mainly intact. Figure 2.8 shows that the equilibrium similarity form suggested by Smyth (2003) can be extended to Pr far above 1. However, unlike the other predictions in figures 2.6 and 2.7, the voticity predictions in figure 2.8 do not improve by the increase in Re0. This is most likely because the similarity form given by equation 2.15 has its limitations in predicting the vorticity structure dominated by a very thin density interface. The increase in Pr, even to high values, creates a marginal increse in braid vorticity, while the increase in Re0 has a more significant effect. Prediction of secondary instabilities on the braid is an important applica- 46 2.6. Similarity solutions 102 103 104 105 0.25 0.5 1 2 4 6 Re0 Ω (t s ) δ 0/ ∆ u (a) 100 101 102 103 Pr (b) Figure 2.8: Comparison of the similarity solution (solid line) and DNS (trian- gles) braid centerline vorticity at saturation for (a) different Reynolds numbers and a constant Prandtl number of 0.72 and (b) different Prandtl numbers and a constant Reynolds number of 400. Ω is averaged over x = 0 to x = Lx/8. The vertical dashed lines show the range of the mixing transition Reynolds numbers found by Rahmani et al. (2011a). tion of the vorticity similarity solution. KH-like secondary instabilities origi- nating from the braid stagnation point have been found by Smyth (2003) and Staquet (1995) on the braid of vortex pairing billows, where vorticity is suffi- ciently high. Smyth (2003) show that Ω/γ > 35 ∼ 40 on the stagnation point of the braid is a sufficient condition for the generation of secondary instabili- ties. This critical ratio is 54 in simulations of Staquet (1995). In the absence of pairing and for the low value of J , none of the simulations 1 to 10 devel- oped secondary instabilities as Ω/γ < 12 for the range of Reynolds numbers studied. We used the CS76 model to predict that Ω/γ > 35 on the braid of a single wavelength is achievable for Re0 = 240, 000 and J = 0.2 (simulation 12 of table 2.2). Secondary KH-like instabilities develop for this case (as shown in section 7). This further supports the critical Ω/γ value of 35 ∼ 40 suggested by Smyth (2003). 47 2.7. Small length scales 2.7 Small length scales One of the most important aspects of DNS is choosing an appropriate grid spacing, ∆x, with the goal of reproducing the smallest relevant scale of motion, which is not known in advance. In this section we investigate the possible use of CS76 as a guide for choosing ∆x. Here, we consider the following length scales: LB, L ′ T , LD, LM and LCS, defined as follows. The Batchelor length scale, LB, is the scale below which diffusion represses scalar gradients. We compute this length scale globally for the whole domain using equation 2.8. L ′ T is the largest length in the scalar dissipation range that is significantly influenced by diffusion. L ′ T is the equivalent definition of the Taylor scale, LT (Tennekes & Lumley, 1972), in the scalar field, and found from: L ′ T = √ 15Re 1/4 0 LB. While LB and L ′ T are defined theoretically, LD and LM are measured directly from spectra of the stream-wise scalar field. The dissipation length scale, LD, is defined as the length below which 0.1% of the total scalar dissipation is not captured. The length scale of maximum dissipation, LM , is the scale at which maximum dissipation of scalar variance occurs. For a pure strain field with strain rate γ, viscous dissipation rate is ε = 4νγ2. Substitution in the definition of LB, defined by equation 2.8, yields a length scale, that we shall refer to as the Corcos-Sherman scale LCS(t) = ( κ 2γ(t) )1/2 . (2.17) The novelty of this definition is that it uses the strain rate invoked by the large-scale core vortex on the braid and estimated by the CS76 model as γ, rather than estimating the order of magnitude of the strain rate as in LT and LB. Therefore, the Corcos-Sherman length scale is a more deterministic estimation of the small scale of the motion. Conceptually, LCS is similar to LT as they both use the large-scale or “outer-scale ”strain rate (e.g. see Su & Clemens, 1999; Dimotakis, 2005). However, LCS is derived based on a pure two-dimensional strain field and is equal to the local Batchelor scale of the 48 2.7. Small length scales braid since it integrates the local dissipation on the braid into the definition of the Batchelor scale. As also pointed out by Corcos & Sherman (1984) and Lin & Corcos (1984), the special properties of the braid create a unique mechanism for the generation of small scales in this region. The strain on the braid, with the highest rate in the vicinity of the stagnation point, is perpendicular to the density gradient. This produces a highly efficient thinning of the density interface (Smyth, 1999). This condition is persistent during the saturation of the billow, making the saturation braid density thickness the smallest scale of the whole domain. At saturation, LCS is related to the braid thickness as hb(ts) = 2pi 1/2Lcs(ts) (compare equations 2.12 and 2.17). Hence, LCS , which is defined based on the thinning of the braid and proportional to the braid thickness, is a relevant measure of small scales in the two-dimensional pre-saturation flow. While the braid thickness can confidently be considered as the smallest length scale before saturation, it is not necessarily the case after saturation and during the turbulent phases due to the more complex billow structure and the appearance of smaller scale motions. However, as discussed by Smyth (1999) these small-scale motions are non-persistent and their strain field is not perfectly aligned with their density gradient direction. Thus, it still can be speculated that braid thickness reached at saturation under the most efficient thinning condition is comparable to the smallest scale attained during turbu- lent phases. The evidence for this speculation can be qualitatively seen in the snapshots of the scalar dissipation rate, χ, in figure 2.9, where χ is defined as χ = 2κ 〈∇ρ.∇ρ〉 , (2.18) with 〈〉 denoting volume averaging. At saturation, the dissipation occurs on the interface of the rolling billow, with its majority in the braid region. At maximum intensity of turbulence, where the viscous dissipation rate of kinetic energy, ε, reaches a maximum, the dissipation is concentrated on thin disor- ganized dissipation layers, studied in detail by Buch & Dahm (1996, 1998); 49 2.7. Small length scales x (m) z (m ) 0.2 0.4 0.6 0.8 1 1.2 1.4 0.4 0.6 0.8 1 1.2 x (m) z (m ) 0.2 0.4 0.6 0.8 1 1.2 1.4 0.4 0.6 0.8 1 1.2 Figure 2.9: Snapshot of log10χ, the scalar variance dissipation rate, for Re0 = 1200 and Pr = 9 at (a) t∗ = 35, the saturation time, and (b) t∗ = 116, the maximum intensity of turbulence. Su & Clemens (1999); Kuthnur & Clemens (2005). These studies show that the structure of these layers is in agreement with the similarity solutions of the diffusion of a strained layer, explained in section 6. The thickness of the braid dissipation layer at saturation is close to the thickness of the most dom- inant dissipation layers at maximum intensity of turbulence. This hints that the effective strain rate multiplied by its time of action has been almost the same in the formation of the dissipation layers at the two stages shown. As pointed out by Buch & Dahm (1996) and Buch & Dahm (1998) the thickness of the dissipation layers do not vary over a wide range. Hence, LCS, which is h/2pi1/2, is a relative measure of the fine dissipation length scales at the turbulent stage in figure 2.9. We examine this indication more quantitatively by examining LD and LM , representing the minimum and maximum scales of the scalar dissipation, at turbulent stage of three-dimensional simulations. To get LD and LM , scalar spectra are built by applying a fast Fourier transform to the stream-wise scalar field at each (y, z) point and then aver- aging in the vertical and span-wise directions. In terms of the stream-wise 50 2.7. Small length scales wavenumber, k, this scalar spectrum is ψ(ki) = (ρ̂) 2, 2pi Lx Nx/2∑ i=1 ψ(ki) = 〈ρ2〉xyz, (2.19) withˆdenoting the Fourier transform. The dissipation rate spectrum is then 4pik2i κψ(ki), and the volume averaged scalar dissipation rate is found from χ = 4piκ Lx Nx/2∑ i=1 k2iψ(ki), (2.20) Following Smyth (1999), we non-dimensionalize the spectra as ψ∗ = ψκk3B χ , k∗ = k kB , (2.21) where kB = 1/LB is the Batchelor wavenumber. The non-dimensional spectra in a variance preserving form are presented in figure 2.10 for the lowest and highest tested Reynolds numbers at Pr = 9. Also, the theoretical Kraichnan spectrum is fitted to each DNS spectrum (see Smyth, 1999, for details of the fitting). The spectra from DNS are cut at 0.8 of the Nyquist frequency to eliminate aliased wavenumbers (Moser & Moin, 1987; Smyth & Moum, 2000). The area under the spectrum between the lowest wavenumber and a given wavenumber indicates the total scalar dissipation obtained down to that scale. As pointed out by Moin & Mahesh (1998), the adequacy of a grid resolution can be evaluated by the amount of dissipation it captures. Therefore, LD, the length scale at which 99.9% of dissipation is recovered, is a stringent resolution requirement. When the spectra are cut before this scale, the constructed Kraichnan part of the spectra is used for LD calculation instead. LM is measured from the theoretical spectra. Different length scales are compared in figure 2.11. LD and LM , the scales measured from the spectra, reach the Batchelor scaling slope of Re −3/4 0 for Re0 > 600. This scaling is consistent with the results of Buch & Dahm (1998) 51 2.7. Small length scales 10−2 10−1 100 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 ↓ k*Bk * CS ↓ k*M ↓ k*D k* Ψ * k* 3 10−2 10−1 100 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 ↓ k*B ↓ k*CS ↓ k*M ↓ k*D k* Ψ * k* 3 Figure 2.10: Spectra of scalar variance at maximum intensity of turbulence for (a) Re0 = 400 and Pr = 9 and (b) Re0 = 1200 and Pr =9. The constructed Kraichnan spectra are shown by dashed lines. The arrows indicate the location of non-dimensional wave-numbers corresponding to different length scales. 52 2.7. Small length scales and Su & Clemens (1999) that show the ratio of the average layer thickness to the outer length scale, i.e. δ0 here, scales with Re −3/4 0 Pr −1/2. LM gives a measure of the largest scale in the dissipation range. As expected from the definition of the Taylor scale, L ′ T is close to LM , however it varies at a slope of Re −1/2 0 . As also remarked by Moin & Mahesh (1998), LB is a conservative resolution requirement since it remains at least three times smaller than LD. In this sense, LCS is also considered a stringent measure as it is well below LD, although it varies at the Taylor scaling slope of Re −1/2 0 . Since LCS is computed based on the outer scale strain rate, this scaling is expected (Su & Clemens, 1999). The proportionality of LCS to L ′ T is almost: L ′ T = 16LCS. Hence, LCS is related to the Taylor scale as: LT = 48LCS. Based on the range of Reynolds numbers examined, the figure suggests that LCS will not exceed LD for several decades of increase in Re0. This is also confirmed by the higher Reynolds numbers examined at Pr = 1. Hence, for a significantly wide range of Reynolds numbers LCS is within the dissipation range scales, i.e. LB < LCS < LM , and close to the lower end of the dissipation range, below which the scalar dissipation is small. Resolving the domain to LCS leaves out a very small portion of the total scalar dissipation, as shown in figure 2.12. This portion is greater than the one left out below LB, but less than 0.1% for all Reynolds numbers at Pr = 1 and 9. The percentages of the dissipation left out below each length scale do not show a significant dependence on the Prandtl number. Considering that the examined Reynolds numbers are sufficiently large to generate small scale turbulent motions (they are above the mixing transition), it is safe to use LCS as a reasonable measure of fine scales even at a turbulent stages. The significant advantage of LCS is its predictability without the need to perform computationally demanding numerical simulations or rely on empiri- cal scaling for estimations of ε or γ. As seen in sections 5 and 6 this prediction is very accurate for high Re0 and Pr. The calculation of LCS however involves some approximations of the flow field, yet the more deterministic represen- 53 2.7. Small length scales 103 10−3 10−2 10−1 Re0 sc a le (m ) L,T LM LD LCS LB Figure 2.11: Variation of different length scales with the Reynolds number for simulations 12 to 15 of table 2.2, with Pr = 9. Length scales are: LT (thick dotted line), LM (solid line), LD (thin dashed line), LCS (dashed-dotted line) and LB (dotted line). tations of the two-dimensional laminar flow makes the error in LCS minimal. The interesting feature of LCS is that a scale found from a two-dimensional laminar state of the flow is still applicable to a three-dimensional turbulent state. It should be noted that LCS, in the form presented here, has its limitations in predicting the smallest scale for all two-dimensional KH flow fields. An ex- ample is the occurrence of secondary instabilities, as illustrated in figure 2.13. For this case the secondary instabilities generate scales smaller than the braid thickness. However, the density and vorticity fields suggest that the struc- ture of these secondary instabilities resembles those of the primary instability. Hence, their length scales can be predicted by the CS76 model. Our presen- tation of LCS also excludes vortex pairing. Smyth (2003) shows that braid thickness reaches its lowest value before pairing. So, the pre-pairing LCS is 54 2.7. Small length scales 102 103 104 10−12 10−10 10−8 10−6 10−4 10−2 10−1 100 Re Pe rc en ta ge o f s ca la r d iss ip at io n (a) Pr = 9 102 103 104 Re (b) Pr = 1 Figure 2.12: The percentage of scalar dissipation not captured when excluding length scales below LB (solid line) and LCS (dashed line) for (a) Pr = 9 and (b) Pr = 1. Figure 2.13: Snapshot of secondary instabilities illustrated by (a) the non- dimensional density and (b) the non-dimensional vorticity for Re0 = 240,000, Pr = 1 and J = 0.2 at t∗ = 47. still applicable, yet pairing triggers secondary instabilities easier. Investigating these problems is beyond the scope of this paper. 55 2.8. Summary and conclusions 2.8 Summary and conclusions The large and small-scale structure of Kelvin-Helmholtz instabilities are ex- amined through a series of two-dimensional direct numerical simulations ex- tending the simulations of Patnaik et al. (1976) to Reynolds numbers as high as 100,000 and Prandtl numbers as high as 700. The evolution of instabil- ities until saturation, the time at which the billow stops gaining height, is investigated in detail. As proposed by CS76, the structure of the instability is studied by dividing the domain to core and braid regions. The DNS results show that at high Re0 and Pr more distinct core and braid vorticity form and less vorticity re-circulates from core to braid. A transition in large and small-scale billow structure occurs in the same range of mixing transition Reynold numbers found by Rahmani et al. (2011a). The large-scale structure of the core and braid has little dependence on Re0 for Re0 > 1600. Very similar number of rolls and therefore billow structure occur for high Re0. The dependence of the large-scale structure on Pr is less straightforward, but the high Pr structure reveals similar characteristics to those of high Re0 structure. The small-scale structure is dependent on Re0 and Pr. Higher concentrations of vorticity are baroclinically generated as Re0 and Pr increase and this results in more structured vorticity fields. Negative vorticity appear inside the core at high Re0 and Pr. The evolution of the large-scale billow structure is described by the CS76 model by predicting the rate of vorticity depletion from the braid to the core. DNS results show that the CS76 model provides useful estimates of the evo- lution of the instability, even though the real vorticity structure is more com- plicated. Predictions of the model for the braid circulation, shear and strain rate and the billow height are very accurate at saturation for high Re0 and Pr as less vorticity re-circulates from core to braid as Re0 and Pr increase. Simulations at both high Re0 and Pr were not possible due to computational constraints. However, the results of this study suggest that the CS76 model is 56 2.8. Summary and conclusions capable of accurate predictions of such conditions. Given the strain rate and shear on the braid from the CS76 model, sim- ilarity solutions are used to describe the vorticity and density profiles across the braid. These solutions accurately predict the braid density thickness and maximum braid vorticity on the centerline except for low Re0 at low Pr. At very high Re0 and Pr the accuracy of braid centerline vorticity prediction marginally decreases as the braid thickness gets very thin. These predictions of the similarity solutions describe the small-scale structure of the braid. On a single wavelength of KH billow, secondary instabilities occur for very high Reynolds and bulk Richardson numbers. The braid thickness is a significant scale as it is the finest length scale in the domain prior to saturation. Con- sidering that the billow height is a large-scale measure of the billow, the CS76 accurately provides estimates of the smallest and largest scales in the two- dimensional flow. Based on the speculation that the braid thickness is a guide to find the smallest relevant scale, the Corcos-Sherman scale is introduced. This scale is predicted from the CS76 model at any time and asymptotes to the local Batchelor scale on the braid at saturation. The action of the high strain rate on the braid and its persistency makes the braid thickness of the same order of the thickness of the dissipative layers in turbulent flow fields and the Corcos- Sherman scale comparable to Batchelor scale in turbulent flows. For the range of Reynolds numbers tested here, LCS is always between the Batchelor scale, LB, and LM , the scale at which the maximum scalar dissipation occurs, and is proportional to the Taylor scale, LT . By examining the scalar variation spectra of three-dimensional flow fields at the maximum rate of viscous dissipation, we show that a very small portion of the total scalar dissipation is associated with the length scales smaller than LCS. This suggests that LCS is a useful scale in providing a measure of small scales in the dissipation range in turbulent flows and can be used as a guide in DNS. The main advantage of LCS over the other turbulence length scales is that it can be estimated prior to the performance 57 2.8. Summary and conclusions of numerical simulations. The high accuracy of the CS76 model at high Pr and Re0 makes LCS a reliable prediction for such conditions. However, this length scale might have its limitation in application to more complex KH flow fields. Further investigation is required to study these special cases. 58 Chapter 3 The effect of Reynolds number on mixing in Kelvin-Helmholtz instabilities 1 3.1 Introduction The turbulent breakdown of Kelvin-Helmholtz (KH) instabilities is responsible for a significant portion of the vertical transport of scalars and momentum in the ocean (Woods, 1968; DeSilva et al., 1996; Smyth et al., 2001; Moum et al., 2003). KH instabilities form when density stratified flows are sufficiently sheared (see Drazin & Reid, 1981, for the details of the instability conditions). During its entire life-cycle a KH instability goes through different phases. The primary KH instability consists of a series of periodic spanwise vortices that roll up the density interface and form billows (Corcos & Sherman, 1976). After sufficient amplification, the primary instabilities become susceptible to the development of small-scale, three-dimensional, secondary instabilities. These appear as streamwise vortices and precede the turbulent stage of the flow. In the turbulent phase, fluid from both layers is entrained inside the billow and mixes while the primary structure of the billow disintegrates. Finally, turbulent motions are dissipated by viscosity and the flow becomes laminar again (Caulfield & Peltier, 2000; Peltier & Caulfield, 2003). Through this life- 1A version of this chapter has been submitted for publication. M. Rahmani, G. A. Lawrence and B. R. Seymour (2011) The effect of Reynolds number on mixing in Kelvin- Helmholtz instabilities. 59 3.1. Introduction cycle, large-scale motions are mainly responsible for entrainment, while mixing occurs at the smallest scales (Dimotakis, 2000). An important feature of mixing in turbulent flows is the “mixing transi- tion,” which manifests itself by an abrupt increase in mixing beyond a critical Reynolds number and a weak dependence of mixing on the Reynolds num- ber after the transition (Dimotakis, 2000, 2005). For KH billows in mixing layers, this transition is associated with the development of small scale three- dimensional motions that facilitate the mixing by increasing the interfacial area between the scalars (Konrad, 1976; Breidenthal, 1981; Koochesfahani & Dimotakis, 1986; Bernal & Roshko, 1986; Dimotakis, 2000). Konrad (1976) first noted that the rapid increase in three-dimensionality downstream of a splitter plate occurs at the local-outer scale Reynolds number, based on the vertical extent of the billow, of order of 104. Laboratory mixing measurements of Breidenthal (1981) and Koochesfahani & Dimotakis (1986) in a mixing layer showed a transitional range of outer-scale Reynolds numbers between approx- imately 103 and 104. The amount of mixing caused by KH billows was only weakly dependent on the Reynolds number outside this range. In this chapter, our primary objective is to investigate the mixing transi- tion in a single KH billow using direct numerical simulations (DNS). The key to understanding the mixing transition is in quantifying the mixing process through different phases of the lifetime of KH instability. For this purpose DNS is a powerful tool due to its accurate representation of small-scale gradi- ents and mixing. The total mixing and time varying mixing rate through the whole life-cycle of a KH billow have been predicted using DNS by Cortesi et al. (1999); Staquet (2000); Smyth et al. (2001, 2007); Carpenter et al. (2007); In- oue & Smyth (2009). Examination of mixing efficiency has provided valuable insight into the mixing process in DNS studies of Scinocca (1995); Caulfield & Peltier (2000); Smyth et al. (2001); Peltier & Caulfield (2003); Inoue & Smyth (2009). Here we use the framework suggested by Winters et al. (1995) to study the evolution of mixing rate and efficiency throughout the lifetime 60 3.1. Introduction of the instability. Our investigation is a more comprehensive quantification of the dependence of mixing properties on the Reynolds number compared to previous DNS mixing studies. The results enable us to describe the mixing transition phenomenon and the post-transition behavior of the mixing in more detail. We compare the transitional range of Reynolds numbers over which the final amount of mixing in DNS rises rapidly to those found by Breidenthal (1981) and Koochesfahani & Dimotakis (1986). The behavior and amount of mixing caused through different phases of KH instabilities have been investigated in numerous laboratory studies (e.g. Thorpe, 1973; Koop & Browand, 1979; Patterson et al., 2006) and field mea- surements (e.g. Seim & Gregg, 1994, 1995; Smyth et al., 2001; Moum et al., 2003). These studies show strong dependence of mixing properties on the phase of evolution of the KH billow. However, measuring the evolution of mixing during different phases of KH billows in the laboratory and field has limita- tions due to the high resolution requirements, the large spatial dimensions of the evolution of the billows, and three-dimensional nature of the flow. For these reasons DNS is a valuable complementary technique in mixing studies. Previous mixing studies indicate that there are several factors that can po- tentially influence the transitional Reynolds numbers, as well as the behavior and amount of mixing. These factors include: vortex pairing of two adjacent billows (Moser & Rogers, 1991), the strength of stratification (Atsavapranee & Gharib, 1997), the Prandtl number (Konrad, 1976; Breidenthal, 1981), and the initial conditions of the shear flow (Breidenthal, 1981). Although we briefly discuss the effects of some these factors when comparing the DNS and labora- tory results, the detailed investigation of their effects is the subject of future work. This chapter is organized as follows. In section 2 we present a brief overview of the background on KH instabilities. In section 3 the numerical model and choice of simulations are described. The mixing transition in DNS is presented in section 4. Then, we examine the density structure of the billow as it goes 61 3.2. Background through the transition in section 5. Energy partitions are discussed in section 6, and mixing rates and efficiencies in section 7. The DNS mixing is compared to the results of Breidenthal (1981) and Koochesfahani & Dimotakis (1986) in section 8, with the conclusions stated in section 9. 3.2 Background 3.2.1 Relevant dimensionless parameters The initial flow has density and horizontal velocity profiles, with a vorticity interface thickness of δ0 and density interface thickness of h0, representing the initial thickness for general thickness of h(t) and δ(t). The horizontal velocity and the density differences are denoted by ∆U and ∆ρ, with initial profiles expressed as U(z) = ∆U 2 tanh( 2z δ0 ), and ρ̄(z) = −∆ρ 2 tanh( 2z h0 ). (3.1) The flow is governed by four dimensionless numbers Re0 = δ0∆U ν , J = ∆ρgδ0 ρ0∆U2 , P r = ν κ , and R = δ0 h0 , (3.2) the initial Reynolds number, bulk Richardson number, Prandtl number and scale ratio, where ν denotes the kinematic viscosity, κ the molecular diffusivity, ρ0 the reference density and g the gravitational acceleration. 3.2.2 Temporal evolution of a KH billow The life-cycle of a KH instability is divided to four phases, as suggested by Caulfield & Peltier (2000). The end of each phase is marked by a transition time, to which we will refer as t∗1, t ∗ 2, t ∗ 3, t ∗ 4 hereafter, with t ∗ = t∆U/δ0 being the non-dimensional time. The three-dimensional density structure of a simulation 62 3.2. Background with Re0 = 900 and Pr = 1 at these four transition stages is shown in figure 3.1. The instability starts by re-distributing the vorticity of the initial vortex sheet and accumulating vorticity in the core of the billow (Corcos & Sherman, 1976). The core vortex of the billow drives the roll-up of the density inter- face as exhibited in panel (a) of figure 3.1. This process entrains fluids from both layers inside the core, while creating a strained interface on the braids connecting the cores. As a result the billow gains height and the potential energy of the flow increases. This phase ends at t∗1 when the total potential energy of the billow, P , reaches a maximum. We refer to this stage as the “saturation”of the two-dimensional growth of the instability. In the absence of pairing, the two-dimensional billow continues to evolve by exchanging en- ergy with the mean flow and generating more mixed fluid inside the billow through molecular diffusion of the rolled interface. This phase ends with the onset of growth of three-dimensional secondary instabilities at t∗2. We define t ∗ 2 as the time when K3d/K0 first exceeds 10 −3, where K3d is the kinetic energy of the three-dimensional motions andK0 the initial kinetic energy of the flow. As a result of the growth of the three-dimensional instabilities (studied in detail by Klaassen & Peltier, 1985b; Caulfield & Peltier, 2000) small-scale chaotic motions develop in the flow field and more fluid from both layers is entrained inside the billow region. This state is referred to as the growing turbulence phase here. However, it should be noted that at low Reynolds numbers the flow lacks a wide spectrum of energy-containing to energy-dissipating length scales to exhibit developed turbulence. As the turbulent motions grow, they generate higher velocity gradients and therefore more viscous dissipation of the turbulent kinetic energy. After reaching a maximum intensity at the peak of K3d/K0 at t ∗ 3, the decay of the three-dimensional motions starts. In this turbulence decaying phase, the mixing of the density field continues and a re- laminarized stable flow state is reached at the end of the life-cycle, t∗4. At this time the turbulent motions are very small. We define t∗4 as the time K3d/K0 63 3.2. Background P K3d (a) (b) (c) (d) t*1 t * 2 t * 3 t * 4 Figure 3.1: The density structure of the KH billow at different stages at t∗1, t∗2, t ∗ 3 and t ∗ 4 in panels (a), (b), (c) and (d) for a Re0 = 900 instability. These times correspond to the maximum of the potential energy, the onset of three- dimensional instability, the maximum intensity of three-dimensional motions and re-laminarization, respectively. The definition of the different stages is illustrated in the bottom panel using the total potential energy, P , and the three-dimensional turbulent kinetic energy, K3d. drops below 10−3. 3.2.3 Energy partitions and transfer rates Quantitative description of the instability is assisted by an analysis of the exchange rates between different components of the energy budget. Based on earlier work of Lorenz (1955), Winters et al. (1995) discuss a systematic description of the energetics in density stratified flows for a closed system, which will be followed here. A schematic of energy partitions and the transfer rates between them are shown in figure 3.2. All energy partitions are non- 64 3.2. Background dimensionalized by dividing the dimensional energy partitions by ρ0(∆U) 2; en- ergy transfer rates are non-dimensionalized by dividing the dimensional energy transfer rates by ρ0(∆U) 3/h0. The total volume averaged non-dimensional ki- netic energy, K, is K = 〈u.u〉V 2(△u)2 , (3.3) where 〈 〉V denotes an average over the domain volume. Contributions to the total kinetic energy are usually partitioned to three components: the two- dimensional and the three-dimensional turbulent motions, and the mean flow. These are defined respectively as (Caulfield & Peltier, 2000) K2d = 〈u2d.u2d〉xz 2(△u)2 , K3d = 〈u3d.u3d〉xyz 2(△u)2 and K = 〈u.u〉z 2(△u)2 , (3.4) where u(z, t) = 〈u〉xy, u2d(x, z, t) = 〈u− u〉y, and u3d(x, y, z, t) = u− u− u2d, (3.5) and 〈 〉i shows an average in the direction(s) i. Then the turbulent kinetic energy is K ′ = K2d + K3d. The total non-dimensional potential energy is defined as P = g〈ρz〉V ρ0(△u)2 . (3.6) The background potential energy obtained through adiabatic rearrangement of fluid particles is Pb = g〈ρz∗〉V ρ0(△u)2 , (3.7) where z∗ is the new position of the fluid parcel located at z in the adiabati- cally redistributed state with minimum potential energy. In a closed system the background potential energy can change only due to irreversible mixing process. The difference between the total potential energy, P , and the back- ground potential energy is defined as the available potential energy, Pa. Then 65 3.2. Background the non-dimensional available potential energy is Pa = P − Pb. (3.8) The available potential energy represents the portion of the potential energy that can be exchanged with the kinetic energy reversibly. The time rates of changes of the different components of the energy budget provide a tool to investigate the evolution process. As shown in figure 3.2, through the shear production or destruction of turbulent kinetic energy, ψ, energy is extracted from or given back to the mean kinetic energy partition. Meanwhile the mean and turbulent kinetic energy irreversibly lose energy to the internal energy through the mean and turbulent viscous dissipation, ε and ε ′ . The total rate of dissipation of kinetic energy is then ε = ε + ε ′ . The energy gained by K ′ through the shear production can be exchanged reversibly with the available potential energy through the turbulent buoyancy flux, φ ′ b. Because of the flow condition here, the mean kinetic energy cannot exchange energy with the available potential energy (i.e. the mean buoyancy flux is zero). Therefore, the total potential energy can only be increased by the turbulent kinetic energy and the internal energy. The latter occurs due to the molecular diffusion of the background density profile and at the rate φi. Mixing is an irreversible process through which the flow uses the available potential energy to enhance the background potential energy at the rate φM . Therefore the rate of change of the background potential is dPb/dt ∗ = φM+φi. The amount of mixing at any time, M , is measured by the increase in the background potential energy from its initial state minus the diffusion that would occur if there were no mean flow, i.e. M = Pb − Pb(0)− φit∗. The mixing efficiency is an important diagnostic for characterizing mixing. A commonly used measure of the efficiency of mixing is the instantaneous flux coefficient, defined as (Caulfield & Peltier, 2000; Smyth et al., 2001; Carpenter 66 3.2. Background Figure 3.2: A schematic of energy partitions and transfers between them based on the framework suggested by Winters et al. (1995). et al., 2007) Γi = φM ε′ . (3.9) In fully developed turbulent mixing caused by shear instabilities Γi is close to 0.2 (Caulfield & Peltier, 2000; Smyth et al., 2001; Carpenter et al., 2007; Inoue & Smyth, 2009). Also of interest is the cumulative flux coefficient (suggested by Caulfield & Peltier, 2000; Peltier & Caulfield, 2003) Γc = ∫ T φMdt ∗/ ∫ T ε ′ dt∗, (3.10) over the duration T of the mixing event. 67 3.3. Numerical modeling 3.3 Numerical modeling Using the Boussinesq approximation and incompressible fluid assumption the Navier-Stokes and continuity equations in a Cartesian coordinate system take the form Du Dt = −∇P ρ0 − ρ ρ0 gk̂+ ν∇2u, and ∇.u = 0, (3.11) where P denotes the pressure field. D/Dt denotes a material derivative and k̂ the unit vertical vector. These equations along with the advection-diffusion equation Dρ Dt = κ∇2ρ, (3.12) are solved using a DNS pseudo-spectral model described by Winters et al. (2004) and modified by Smyth et al. (2005). In this model the resolution of the density field is twice that of the velocity field. Three-dimensional DNS of the full life cycle of KH billows are performed at different initial Reynolds numbers, see table 3.1. These simulations are chosen to include a wide range of Reynolds numbers to cover the pre- and post- transition range. In all simulations J =0.03, Pr = 1, and R = 1. The chosen value for Pr is close to the Prandtl number of the thermal diffusion in gases. One wavelength of the fastest growing instability, found from linear stability analysis, is simulated. Therefore, we are neglecting pairing (as captured by Winant & Browand, 1974; Smyth et al., 2005; Carpenter et al., 2010b) to particularly focus on the evolution of one billow. The height of the domain Lz = 9δ0, is sufficiently high to ensure negligible vertical fluxes on top and bottom boundaries (Haigh, 1995; Smyth et al., 2001; Inoue & Smyth, 2009). The spanwise length of the domain is Ly = Lx/2, which accommodates at least one wavelength of the three-dimensional secondary instability (Thorpe, 1985). The required resolution for the DNS of stratified flows with Pr 6 1 is estimated from the Kolmogorov length scale, LK = (ν 3/ε)1/4. (3.13) 68 3.3. Numerical modeling Simulation Re0 Lx/δ0 Nx ×Ny ×Nz ∆z/LK 1 100 8.3 64× 32× 64 0.7 2 200 7.7 128× 64× 128 0.5 3 300 7.5 128× 64× 128 0.6 4 350 7.5 128× 64× 160 0.6 5 400 7.4 128× 64× 160 0.8 6 600 7.3 128× 64× 160 1.1 7 900 7.2 128× 64× 160 1.5 8 1200 7.2 160× 80× 192 1.7 9 2400 7.1 192× 96× 256 2.0 9∗ 2400 7.1 256× 128× 320 1.6 10 4800 7.1 256× 128× 320 2.5 Table 3.1: Description of the simulations carried out. In all simulations J = 0.03, Pr = 1 and Lz = 9δ0. The presented mesh sizes are for the velocity field; the density field mesh sizes are two times larger. ∆z represents the vertical grid spacing of the velocity field. In all our simulations the grid spacing of the velocity field (the coarse resolu- tion) is smaller than 2.5LK throughout the entire duration of the simulation, which satisfies the requirements suggested by Smyth & Winters (2003); Moin & Mahesh (1998). Also, a resolution sensitivity test was performed for simu- lation 9 (Re0 = 2400), as described by simulation 9 ∗. Increasing the number of grid points by 122% in simulation 9∗, resulted in a less than 4% increase in the final amount of mixing. Periodic horizontal and free slip vertical boundary conditions were ap- plied. Basic velocity and density profiles in equations 4.1 are perturbed using the eigenfunctions from the numerical solution of the linear stability equation. The magnitude of this perturbation is large enough to cause a maximum dis- placement of 0.1δ0 of the interface, and small enough for the perturbation to be linear. The amplitude of the two-dimensional perturbation had a negligible effect on mixing. Also, a random noise with the magnitude of ±0.05∆U is added to the velocity field to initiate three-dimensional secondary instabili- 69 3.4. Mixing transition in DNS 102 103 104 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Re0 M /K 0 Figure 3.3: Variation of the amount of mixing gained at the end of the KH life-cycle normalized by the initial kinetic energy, i.e. M(t∗4)/K0, with the Reynolds number (solid line). The dashed line shows the mixing gained during the two-dimensional phases, i.e. M(t∗2)/K0. ties. Lowering the amplitude of the initial three-dimensional perturbation by a factor of 4 resulted in small, but inconsequential, changes in mixing. A de- tailed discussion of the effect of the initial perturbations on mixing is given by Smyth et al. (2007). 3.4 Mixing transition in DNS In this section we describe the mixing trend seen in DNS and its special fea- tures. The variation in the total amount of mixing that the evolution of KH instability has produced until the re-laminarization stage at t∗4 (note that this time is different for different Re0) with the initial Reynolds number is pre- sented in figure 3.3. This figure also shows the variation of the amount of mixing induced in the two-dimensional phases up to t∗2 with Re0. 70 3.4. Mixing transition in DNS The total mixing shows an abrupt increase for Re0 > 300. The total background potential energy gained due to the mixing rises to 6% of the initial kinetic energy at Re0 = 400 from 2% at Re0 = 300. This trend of increase in total mixing continues until a maximum mixing ratio of near 9% is reached at Re0 = 900. Beyond this Reynolds number of maximum mixing, a weaker but non-monotonic dependence of mixing on Re0 occurs. Similar to the laboratory mixing measurements of Breidenthal (1981) and Koochesfahani & Dimotakis (1986) the mixing trend here shows a transition that starts at Re0 = 300 and completes at Re0 = 900. After the completion of the transition however, the DNS mixing does not reveal an Re0 independent behaviour as suggested by these laboratory measurements. Before the start of the transition, i.e. for Re0 6 300, the total mixing consists of the mixing generated during the two- dimensional evolution of the billow; the beginning of the transition is marked by the significant rise of the total mixing beyond the two-dimensional mixing obtained by t∗2. The mixing behavior in DNS can be better understood by examining the time evolution of the gain in the background potential energy for different Re0, presented in figure 3.4. Below the transition (Re0 ≤ 300), the background potential energy rises due to the two-dimensional growth of the billow, then reaches a plateau as the billow is completely mixed after the saturation of the two-dimensional growth and no more mixing follows. Above the transition, the increase in the background potential energy due to the two-dimensional growth is succeeded by a greater rise associated with the mixing caused by the growth of three-dimensional instabilities. The former rise is almost the same for all billows with Re0 > 300, while the latter varies significantly with Re0. This suggests that the variation in total mixing with Re0 most significantly originates from the difference in mixing generated after the appearance of three-dimensional motions. The rise in the background potential energy is enhanced by the growth of small scale three-dimensional motions until Re0 = 900. At Re0 = 1200 the 71 3.4. Mixing transition in DNS 0 100 200 300 400 500 600 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 300 400 600 900 1200 2400 4800 t* M /K 0 Figure 3.4: Comparison of the time evolution of the gain in the background potential energy due to mixing at different Reynolds numbers. The number next to each curve shows Re0. rise in Pb becomes slower than that at Re0 = 900 after t ∗ = 120, a time in the decay phase. The two-step rise in the background potential energy vanishes for Re0 > 1200 to form a smooth gradual growth. A considerable increase in the lifetime of the instability and slower gain in the background potential energy result. These effects lower the total amount of mixing at Re0 = 2400 and 4800. However, at Re0 = 4800 the persistent growth compensates for the slow growth of the background potential energy and the total mixing rises slightly above that of Re0 = 2400. The time variation of the background potential energy suggests that the mixing caused at the end of the life-cycle of a KH instability is a function of 72 3.5. Density structure of the billow the flow condition at different stages, without a straightforward dependence on Re0. Examining the density structure of the billow in the next section links the mixing behaviour revealed in figures 3.3 and 3.4 to the qualitative description of the flow. 3.5 Density structure of the billow The effect of the Reynolds number on the density structure of the billow is significant through all different phases of the evolution of the instability. Prior to the generation of three-dimensional instabilities at t∗2, a higher Re0 produces a less diffuse density structure of the billow with sharper density gradients. Hence, at higher Re0 the billow is more structured and the local baroclinic torques are stronger when the three-dimensional instabilities start to grow. The more significant effect of Re0 is evident after the development of three- dimensional motions; the smaller scale vortices at higher Re0 result in a more developed turbulence inside the billow where mixing is facilitated by smaller scale density gradients. As this effect is displayed well at t∗3, we show the density structure at this stage for increasing Reynolds numbers in figure 3.5. Before the start of the mixing transition, i.e. Re0 < 300 the flow is not energetic enough to overcome viscous effects and sustain three-dimensional in- stabilities. Klaassen & Peltier (1985b) suggest that there is a critical Reynolds number below which secondary instabilities do not form, and this critical Reynolds number depends on J and Pr. In our DNS, for Re0 > 300 the three- dimensionality starts by a combination of stream-wise vortices on the braid and convective instabilities on the edges of the billow. This type of secondary instability is expected as the bulk Richardson number, J , is small (Caulfield & Peltier, 2000). At the stage shown in figure 3.5 the three-dimensional instabil- ities have led the flow to a state of chaotic small scale motions. The smallest scale of these motions is estimated by the Kolmogorov scale, defined by equa- tion 3.13, and can be related to the Reynolds number by: LK ∼= CKδ0Re−3/4, 73 3.5. Density structure of the billow Figure 3.5: Density field at t∗3, maximum intensity of turbulence, for Re0 = 300, 400, 600, 900, 2400, 4800. with Ck being a constant (Tennekes & Lumley, 1972). Therefore, the smallest scale at Re0 = 4800 should be approximately 6.4 times smaller than that at Re0 = 400. As the small scales become finer for higher Re0 they generate more mixing in the billow region and smooth out the sharp density gradients. The mixing inside the billow destroys the density structure of the billow. Nonethe- less, the density field sustains its three-dimensional structure although with different intensity at different Re0. In the process of the transition, e.g. at Re0 = 400 and Re0 =600, the three- dimensionality of the density structure is strong, yet length scales of variations in the span-wise direction are very limited. The entrainment in these cases is governed by the large-scale eddy that still carries the structure of the initial stream-wise vortices. AtRe0 = 900, the Reynolds number of maximum mixing, a range of different scale vortices are generated, indicating more developed turbulence. Some residues of the initial three-dimensional instabilities that 74 3.5. Density structure of the billow triggered the turbulence are still preserved. Through the action of the larger vortices unmixed fluid is entrained inside the turbulent billow region. At Re0 = 900 these vortices are strongly engaged in the entrainment process. Interestingly, the span-wise two-dimensional vortices still play an important role in the entrainment. As the trend of cascade of large to smaller sized eddies continues with the increase in Re0, the entrainment reduces for Re0 > 900. At Re0 = 2400 the two-dimensional structure of the billow has been mainly destroyed and the entrainment occurs by means of stream-wise and span-wise vortices that are significantly diminished in size. We associate the reduction in the vertical extent of the billow and therefore the entrainment with a comprehensive development of small-scale eddies in the billow that destructs the large-scale structure of the billow. Also, we speculate that due to the disappearance of large density gradients the energizing effect of the baroclinic torques (Cortesi et al., 1999; Staquet, 2000) on three-dimensional motions weakens. Hence, the strength of large-scale two and three-dimensional structures is suppressed by fully development of small-scale motions at high Re0. Based on an analysis by Broadwell, Breidenthal (1981) suggests that the mixing may be limited by either the diffusion of small scales or the entrain- ment. The ratio of the small-scale diffusion to the large-scale entrainment time scales varies with PrRe−1/2. Hence, he proposes that for large values of PrRe−1/2 mixing is diffusion limited, while for small values of PrRe−1/2 mixing is entrainment limited. It is clear from figure 3.5 that the three-dimensional snapshots encompass the diffusion limited to entrainment limited mixing be- haviour from Re0 = 400 to Re0 = 4800, with the ratio of PrRe −1/2 dropping from 0.05 to 0.014. At Re0 = 400, while the large-scale eddy is strong enough to entrain sufficient fresh fluid inside the billow, the smallest-scale eddies are too large in size and therefore have long diffusion time scales, limiting the mixing. Contrarily, at Re0 = 4800, the small-scale eddies are much smaller than the large-scale eddies and have short diffusion time scales. In this case 75 3.5. Density structure of the billow the large-scale eddy cannot provide the turbulent region with fresh fluid fast enough. A special property of the Re0 = 900 density field is that the Reynolds number is high enough to generate small-scale eddies essential for an enhanced mixing, yet low enough not to induce rapid mixing by fully developed small scales and limit the entrainment. Theoretical turbulent mixing models, proposed by Broadwell & Breidenthal (1982) and Broadwell & Mungal (1990), indicate that mixing occurs at two stages. During the first stage, unmixed fluid is entrained inside the billow and stirred by streamwise and spanwise motions. The relevant length scale of the diffusive layers at this stage is the Taylor length scale, LT ∼ √ 15δ0Re −1/2 0 . Through the stirring process, the interface between fluids is stretched until the Kolmogorov (or Batchelor for Pr > 1) scale, LK ∼ δ0Re−3/40 , is reached. At this scale, the second stage of mixing or molecular diffusion ensues. Broadwell & Breidenthal (1982) and Broadwell & Mungal (1990) show that the second stage of mixing occurs much faster than the first stage. For fully turbulent flows, the two stages of mixing occur at distinguishably different length and time scales. At low Reynolds numbers, however there is a narrow spectrum of length scales, and the second stage takes place before the first stage is complete. Therefore, a high ratio of LT/LK indicates a separate two-stage mixing, while a low ratio of LT/LK shows a one-stage mixing with a single important length scale. For simulations in figure 3.5, LT/LK varies between 16 and 32, when increasing Re0 from 300 to 4800. This suggests that even at Re0 = 4800, LK and LT are only one order of magnitude different, and it is possible that turbulence is not yet fully developed. An intersting speculation of the model of Broadwell & Mungal (1990) is that mixing goes to zero in the limit of Re0 → ∞. While figure 3.5 suggests a reduction in mixing for high Re0, a generalized conclusion for Re0 → ∞ cannot be drawn here. Examining the energy partitions in the next section provides a more quan- titative evaluation of our entrainment and mixing observations so far. 76 3.6. Energy partitions 3.6 Energy partitions In this section we discuss the evolution of the energy partitions through differ- ent phases of the lifetime of the billow. Figure 3.6 shows the time variation of the background and available potential energy, and two and three-dimensional turbulent kinetic energy partitions added together respectively. The gain in the background potential energy due to the diffusion because of the top and bottom density difference, i.e. φit ∗, is excluded. The summation of all these four energy partitions shows the total amount of energy that is extracted by the instability from the background flow and not lost to the internal energy through the viscous dissipation. The summation of the background and avail- able potential energy partitions indicates the portion of this extracted energy that is given to the total potential energy through the turbulent buoyancy flux. The increase in the total potential energy is due to the entrainment by buoyancy fluxes (see figure 3.2), while the increase in the background potential energy is owing to the mixing of the entrained fluid. Therefore, small available potential energy implies entrainment limited mixing. By the saturation time, t∗1, the two-dimensional entrainment by the roll-up of the billow is complete and P reaches a peak. This is a diffusion limited process as Pa remains high for all Re0. The two-dimensional entrainment produces very little mixing until t∗1 at high Re0 as the roll-up occurs much faster than the diffusion. Conversely, below the transition, at Re0 = 300, a large portion of the total mixing is gained by this time. No significant variation occurs at later stages for this Re0. For all Re0, K2d reaches a peak before t∗1, indicating that the entrainment continues to increase for a short time after the two-dimensional motions have started to decay. The peak in K2d monotonically rises with Re0 as the two-dimensional flow generates more shear production with increasing Reynolds number. The second phase shortens as Re0 rises due to the faster growth of the three-dimensional instabilities. In this phase a portion of the energy extracted by K ′ is returned to K. Since the mixing does not act fast enough to consume 77 3.6. Energy partitions 0 0.05 0.1 Re0 = 300 Re0 = 400 0 0.05 0.1 0.15 Re0 = 600 Re0 = 900 0 100 200 300 400 500 600 0 0.05 0.1 0.15 t* Re0 = 2400 0 100 200 300 400 500 600 t* Re0 = 4800 Figure 3.6: Time evolution of the energy partitions normalized byK0: the lines from bottom to top show the background potential energy, Pb, the available potential, Pa, the two-dimensional turbulent kinetic energy,K2d, and the three- dimensional turbulent kinetic energy, K3d partitions added together. Pa and K3d partitions are shaded. The vertical lines show the transition times. 78 3.6. Energy partitions Pa, a part of the gained P is also returned to K ′ . So the process of energy extraction from the background shear flow reverses. Meanwhile the billow mixes the entrained fluid by the two-dimensional roll-up through the enhanced diffusion, while no fresh fluid is entrained inside the billow. As a result Pa diminishes to small values, indicating an entrainment limited pre-turbulent mixing. The third phase, the turbulence growing phase, also shortens as Re0 rises, indicating a faster saturation of turbulent motions at higher Re0. However, the turbulent motions are most energetic at Re0 = 900. In this phase as K3d starts to grow, it first acts as a sink for K2d. Later when the three-dimensional motions intensify, they give rise to K2d through the shear production they generate. At its peak, K3d reaches about 10% of K0 and is the largest energy partition. The largest portion of the total extracted energy is in the form of turbulent kinetic energy. The new entrainment starts in the third phase by the growth of three-dimensional motions and causes a slow gain in P . In this phase the already completely mixed billows with Re0 < 900 start to gain Pa, while little mixed billows with Re0 > 900 spend more Pa for mixing. In spite of these differences in behaviour, the total mixing at the end of stage 3 shows a weak dependence on Re0. The majority of the entrainment and mixing occurs in the last phase, the turbulence decay phase. While K ′ dissipates, more of the total extracted en- ergy is transformed to P and thereafter Pb. The amount of increase in P depends on the strength of three-dimensional motions, measured by K3d. The occurrence of the highest peak in K3d and also the widest K3d curve at Re0 = 900 indicates the strongest growth and sustainability of three-dimensional motions at this Reynolds number. At Re0 = 2400 and 4800, K3d drops rapidly after reaching a peak. As seen in section 5, we relate this drop to the major destruction of large-scale structure of the billow by the development of small- scale motions that reduces the shear production of turbulent kinetic energy. However, later the interaction of the energetic small and intermediate scale 79 3.7. Mixing rate and efficiency eddies on the boundaries of the turbulent region with the mean flow generates shear production and the three-dimensional motions remain energetic for a prolonged lifetime. Therefore, for Re0 < 900 the three-dimensional entrain- ment process is fast, while for Re0 > 900 it is slow and long. This long lifetime at high Re0 is in agreement with the field observations of Seim & Gregg (1995). The fast entrainment generates a diffusion limited mixing centered around t∗3 for Re0 < 900. Contrarily, the Re0 > 900 billows spend most of their turbulent mixing phase in entrainment limited state with small Pa. In the case of Re0 = 4800, the decay phase mostly consists of near zero Pa. This indicates that the entrained fluid mixes almost immediately. The two-dimensional structure of the billow (indicated by K2d), in phase 4 is best maintained for Re0 = 900, while it is significantly destroyed for Re0 = 2400 and 4800. As seen in figure 3.5, the sustained two-dimensional structure substantially enhances the entrainment. In summary, examining the energy partitions suggests that the maximum mixing occurs when K2d and K3d are the highest and best sustained in a period shortly before and after t∗3. This is the period in which the maximum entrainment of fresh fluid takes place. This entrained fluid is mixed inside the billow while the turbulence is decaying in phase 4. In the next section we discuss how the mixing rate varies between different phases and how efficiently it occurs. 3.7 Mixing rate and efficiency In the previous sections we showed that the increase in the Reynolds number has some significant, yet non-monotonic, effects on the evolution of the KH billow. Now, we focus more on two important rates of energy exchange (see figure 3.2): the mixing rate, φM , and the turbulent viscous dissipation rate, ε ′ . Figure 3.7 shows the time evolution of these two rates. The instantaneous flux coefficient, Γi, also discussed in this section, is found from the ratio of the 80 3.7. Mixing rate and efficiency two rates, as defined by equation 3.9. The entire mixing event consists of two major parts, each having a peak in the rate of mixing. The first peak corresponds to the two-dimensional roll-up of the primary instability and occurs at or after the saturation, depending on Re0. The second peak occurs either at or slightly after t ∗ 3. The pre-transition flow at Re0 ≤ 300 lacks the second peak as no significant three-dimensionality is generated. As Re0 increases the two peaks of the mixing rate gradually move closer together. This is because for higher Re0 the mixing due to the roll-up is delayed more, and also the three-dimensional entrainment grows earlier. At Re0 = 4800 the two peaks in mixing have become so close that there is almost no distinction between the two-dimensional mixing caused by the roll-up and the three-dimensional mixing caused by the growth of secondary instabilities. As Re0 rises the two-dimensional mixing that occurs during phases 1 and 2 becomes less important compared to the three-dimensional mixing that occurs during phases 3 and 4. As already seen in figure 3.6, the most significant portion of the mixing takes place in phase 4 for all Re0 > 300. For high Re0, i.e. Re0 = 2400 and 4800, the mixing generated during the growth of the turbulence is also a large portion of the total mixing. The occurrence of the highest φM at or slightly after t∗3 indicates that the stretching of the density field by three-dimensional motions is an important mechanism in providing the interfacial area between the scalars to enhance the rate of mixing. For some Re0 the maximum φM occurs slightly after t∗3. This suggests that, although the small-scale motions have started to decay, their interaction and stretching of the density field is still strong enough to enhance the rate of mixing. Through the mixing transition, higher maximum mixing rates are generated for higher Re0. While the mixing rates near t∗3 drop considerably for Re0 > 900, the turbulence remains almost as dissipative as that at Re0 = 900. This is an indication of entrainment limited mixing: small-scale motions cause high dissipation rates while the rate of mixing is suppressed by the reduced extent of large-scale motions. 81 3.7. Mixing rate and efficiency 0 0.5 1 1.5 x 10−4 φ M , ε, /5 Re0 = 300 Re0 = 400 0 0.5 1 1.5 x 10−4 φ M , ε, /5 Re0 = 600 Re0 = 900 0 100 200 300 400 500 0 0.5 1 1.5 x 10−4 t* φ M , ε, /5 Re0 = 2400 0 100 200 300 400 500 t* Re0 = 4800 Figure 3.7: Time evolution of the rate of mixing (solid line) and turbulent viscous dissipation rate divided by 5 (dashed line) for increasing Reynolds numbers. The vertical lines show the transition times. 82 3.7. Mixing rate and efficiency The turbulent kinetic energy dissipation rate follows the trend of the mixing rate in the three-dimensional phases. The ε ′ /5 and φM curves become very close shortly after t∗3. This match, intentionally assisted by dividing ε ′ by 5, remarks the universal value of 0.2 for Γi observed in previous studies near the end of KH life-cycle (Smyth et al., 2001, 2007; Inoue & Smyth, 2009). The time variation of Γi for three selected Re0 are presented in figure 3.8(a). Before the rise in ε ′ , the two-dimensional event is highly efficient with high Γi. Smyth & Winters (2003) and Smyth et al. (2007) discuss the significance of this highly efficient mixing phase. The rise in ε ′ is followed by a rise in φM in figure 3.7. However, the time lag in the increase in φM causes a drop of Γi. This lag is associated with the time required for the development of three-dimensional coherent structures and the action of molecular diffusion. For some Reynolds numbers this effect reverses after t∗3 and high fluctuations of Γi occur in the decay phase. As pointed out by Smyth et al. (2001), Γi does not reveal a clear sensitivity to Re0. However, a better indication of the dependence of the entire mixing event on Re0 is the cumulative flux coefficient, described by equation 3.10, and shown in panel (b) of figure 3.8. When only the decay phase is considered, the most efficient mixing with the highest Γc occurs at Re0 = 900. Over the entire three-dimensional phase of mixing, the trend of Γc is very similar to that of the total mixing, with the same rise of Γc at Re0 = 4800. This high Γc is due to the shifting of the two-dimensional pre-turbulent mixing to the beginning of phase 3 at Re0 = 4800. The trend of Γc suggests that the total amount of the induced mixing correlates with the cumulative efficiency of the mixing. Through the transition in mixing, the efficiency of the mixing improves as the more energetic three-dimensional motions generate more entrainment and therefore mixing with little increase in the total dissipation. The maximum mixing occurs when this mechanism is the most efficient. The post-transition mixing diminishes in the total amount of mixing and overall efficiency. This non-monotonic variation of Γc with Re0 is also observed by Smyth et al. (2001). Also, Shih 83 3.8. Comparison with laboratory experiments 100 200 300 4000 0.2 0.4 0.6 0.8 1 t* Γ i Re0 = 600 Re0 = 900 Re0 = 2400 102 103 104 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 Re0 Γ c Figure 3.8: Mixing efficiency: time evolution of the instantaneous flux coeffi- cient, Γi, in (a) and the cumulative flux coefficients, Γc, over the time duration of t∗2 to t ∗ 4 (solid line) and t ∗ 3 to t ∗ 4 (dashed line) in (b). et al. (2005) show a non-monotonic variation of the mixing efficiency (although using a different definition) with the variation of turbulence intensity. 3.8 Comparison with laboratory experiments In this section we compare the mixing transition from DNS to previous lab- oratory studies. The mixing transition in shear layers was illustrated quan- titatively by the product thickness measurements of Breidenthal (1981) and Koochesfahani & Dimotakis (1986). In the literature the mixing transition is commonly stated in terms of the outer-scale Reynolds number Reδ = δ∆U ν , (3.14) where the outer scale, δ, describes the large scale of the flow. Following Brown & Roshko (1974), we find δ at any time from δ = (∂U/∂z)−1max∆U , with U(z) being the averaged vertical velocity profile. In our DNS, δ/δ0 varies between 3 and 6, while the maximum δ/δ0 reached near 12 in the experiments of Koop & Browand (1979) as pairing occurred. Variation of mixing with Reδ from lab- 84 3.8. Comparison with laboratory experiments 102 103 104 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Reδ δ p /δ 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 M /K 0 Figure 3.9: Variation of mixing with the outer-scale Reynolds number from laboratory measurements and DNS. The product thickness, δP , normalized by the outer scale, δ, from the laboratory measurements is shown for: Breidenthal (1981)’s shear layer experiment with r = 0.76 (circles) and r = 0.38 (squares) and Koochesfahani & Dimotakis (1986)’s shear layer experiment with r = 0.38 (triangles). The total mixing normalized by the initial kinetic energy (solid line), M/K0, from DNS is shown by the solid line. oratory experiments, measured by δP/δ, and DNS, measured by M/K0, are compared together in figure 3.9, where δP is the chemical product thickness (see Breidenthal, 1981, for definition). When stated in terms of the outer-scale Reynold number, the laboratory mixing transition depends on the ratio of ve- locity of slower to faster flowing streams, r. The laboratory mixing transition trend of Koochesfahani & Dimotakis (1986) with r = 0.38 is very close to that of Breidenthal (1981) with the same value of r. At r = 0.38, the laboratory mixing transition starts at an outer-scale Reynolds number of 3000 and com- pletes at 104. In DNS, the transition starts at Reδ = 900 and ends at Reδ = 5200, and is between the transition trends at r = 0.76 and r = 0.38. 85 3.8. Comparison with laboratory experiments In summary, the DNS mixing transition shows a similar trend to those of previous laboratory experiments. In DNS, the outer-scale transitional Reynolds numbers are within the range found by Breidenthal (1981) and Koochesfahani & Dimotakis (1986) for different initial conditions. However, the differences between the DNS and laboratory experiments should be noted in this com- parison. Initial and boundary conditions in DNS are different from the exper- imental setup, we expect this difference to account for the shifting of the DNS mixing transition from the experimental data in figure 3.9. The dimensionless parameters, J and Pr are 0.03 and 1 in DNS, while they are 0 and 600 in the experiments. From the experiments of Konrad (1976) with Pr ∼ 0.7 it can be concluded that the decrease in Pr increases the amount of mixing, but does not affect the transition range. It is known that the increase in J decreases the amount of mixing (e.g. Atsavapranee & Gharib, 1997), and from the experi- ments of Atsavapranee & Gharib (1997); Pawlak & Armi (1998) it can roughly be inferred that J does not significantly change the transition range as billows with different J all fully develop turbulence at Reδ ∼ 104. Therefore, we do not expect the difference in Pr and J to significantly influence our compar- ison of the transitional Reynolds numbers in figure 3.9. Another difference between DNS and laboratory mixing transition is the stage at which mixing is measured. We address this difference in the following. To further check the relevance of our mixing transition comparison in figure 3.9, we compare the stage at which mixing is measured in DNS and laboratory. The final DNS mixing is computed at t∗4, at which the mixing is complete by M/K0 asymptoting to a constant (figure 3.4). This stage is reached at later times as Re0 rises due to the longer lifetime at higher Re0. However, in laboratory experiments of Breidenthal (1981) and Koochesfahani & Dimotakis (1986) mixing is measured at fixed locations. As pointed out by Breidenthal (1981), this might mean measuring mixing at a stage of incomplete mixing. Assuming that the location in experiments is related to the time in DNS by: x = Ut, with U being the mean velocity, we have: t∗ = 2x(1 + r)/δ0(1 − r). 86 3.9. Summary and conclusions The DNS results in figure 3.4 show that the mixing is complete at t∗4 = 50 for Re0 = 100, and t ∗ 4 = 350 for Re0 = 900. Testing the location of mixing measurement in shear layers in the experiments of Breidenthal (1981) shows that these locations represent a state of complete mixing in DNS. Therefore, the comparisons between DNS and experiment are made at the same stages. Breidenthal (1981) speculates that the mixing is complete at lower x/δ0 for higher Re0. However, the DNS does not support this speculation. 3.9 Summary and conclusions Using direct numerical simulations, we have studied the effect of the initial Reynolds numbers in the range of 100 to 4800 on the evolution and final amount of mixing in KH instabilities. A monotonic transition in mixing for Re0 = 300 to 900, and a non-monotonic variation of the post-transition mixing for Re0 > 900 is observed. We explain these two phenomena by examining the properties of the flow at different stages of the life-cycle. An abrupt increase in mixing induced at the end of the life-cycle of a KH billow occurs when the initial Reynolds number increases above 300. As dis- cussed by previous studies, this transition in mixing is due to the growth and sustainability of three-dimensional instabilities. While the two-dimensional roll-up of the billow entrains fluids from both layers inside the billow and generates mixing through enhancing the scalar interfacial area, the more sig- nificant portion of the mixing is caused by the three-dimensional turbulent motions. The growth of three-dimensional secondary instabilities generates large-scale eddies that entrain substantial amount of fluid inside the billow region. Small-scale motions that follow these secondary instabilities enhance the molecular mixing. The majority of the entrainment and mixing caused by this mechanism occurs in the turbulence decay phase. In the process of the mixing transition, the entrainment by large-scale eddies increases as the three-dimensional motions become more energetic, and the mixing increases 87 3.9. Summary and conclusions as smaller scale motions develop in the flow. This trend of an increase in entrainment and mixing increases until Re0 = 900. So, our mixing transition range is Re0 = 300 to 900. For Reynolds numbers above the transition, i.e. Re0 > 900, the total amount of mixing first reduces, then slightly increases for Re0 = 4800. This non-monotonic variation of mixing with Re0 has not been observed in previous laboratory experiments. The reduction in the amount of mixing for Re0 > 900 is due to a significant decrease in the entrainment by the three-dimensional motions. As Re0 increases above 900, the development of small-scale motions destroys the large-scale structure of the billow. This significantly decreases the size and strength of large-scale two and three-dimensional motions. Therefore, the development of small-scale motions in the turbulent region acts in favor of entrainment and mixing until an optimum Re0 of 900, with a reverse effect be- yond this Reynolds number. However, the lifetime of the billow is significantly prolonged after the mixing transition. This longer lifetime compensates for the lower mixing rates to some extent and enhances the final amount of mixing. Therefore, the overall mixing induced varies little with Re0 once the transition is complete: the total gain in the background potential energy due to the mix- ing varies between 7 to 9% of the initial kinetic energy for the post-transition mixing. Comparison of our DNS mixing transition to the experimental mixing tran- sition obtained by Breidenthal (1981) and Koochesfahani & Dimotakis (1986) shows a DNS transition trend similar to the laboratory observations. The outer-scale Reynolds number for the completion of the mixing transition in DNS is almost 2 times smaller than the outer-scale transition Reynolds num- ber of 104 often referred to in the literature. However, as found by Breidenthal (1981), the transitional Reynolds numbers vary depending on the flow condi- tions. It should also be noted that the transition in mixing takes place over a range of Reynolds numbers that can take almost one decade of the outer-scale Reynolds numbers. 88 3.9. Summary and conclusions The results of this study also show some interesting general properties of the turbulent mixing. The cumulative flux coefficient over the entire three- dimensional stage shows a close correlation with the total amount of mixing, with an optimal value at the Reynolds number of maximum mixing. At max- imum efficiency, only around 15% of the total energy extracted by the billow during three-dimensional phases converts to mixing (for a maximum cumula- tive flux coefficient of 0.18). The diffusion and entrainment-limited mixing, exhibited by our low and high Reynolds number billows, both have low cu- mulative flux coefficients. Therefore, in addition to the ability of the flow to extract energy from the background shear the efficiency of the mixing is important in determining the total amount of mixing. 89 Chapter 4 The effect of Prandtl number on mixing in Kelvin-Helmholtz instabilities 1 4.1 Introduction In geophysical flows the Prandtl number, the ratio of viscosity to molecular diffusivity of scalars, varies over a wide range from nearly 1 in the atmosphere up to 1000 for salt water in the ocean (usually referred to as the Schmidt number for salt water). However, the effect of Prandtl number on mixing is still not well understood. This is mainly because high Prandtl number mixing occurs at very small scales and requires high resolution numerical simulations or field and laboratory measurements. The importance of molecular diffusiv- ity in turbulent mixing has been demonstrated in the experiments of Turner (1968); Jackson & Rehmann (2003) and numerical simulations of Merryfield et al. (1998); Gargett et al. (2003); Smyth et al. (2005). The results of these studies indicate that a higher Prandtl number induces lower mixing rates. However, this effect becomes less significant as turbulence intensifies at higher Reynolds numbers. Turbulent flows manifest a transition in mixing at sufficiently high Reynolds numbers (Dimotakis, 2000). In Kelvin-Helmholtz (KH) billows this transition 1A version of this chapter is in preparation for publication. M. Rahmani, G. A. Lawrence and B. R. Seymour (2011) The effect of Prandtl number on mixing in Kelvin-Helmholtz instabilities. 90 4.1. Introduction starts at Reynolds numbers high enough to sustain three-dimensional motions and ends at or above Reynolds numbers at which turbulence is fully developed (Konrad, 1976; Breidenthal, 1981; Koochesfahani & Dimotakis, 1986; Rah- mani et al., 2011a). Comparison of the mixing transition observed by Konrad (1976) at Prandtl number of 0.72 with that observed by Breidenthal (1981) and Koochesfahani & Dimotakis (1986) at Prandtl number of 600 suggests that increasing the Prandtl number does not affect the range of the transition Reynolds numbers, but does reduce the amount of mixing. The effect of Prandtl number on different aspects of the evolution of KH billows has been examined in different studies, mainly focusing on the lower range of Prandtl numbers (1. Pr .10). The results of two-dimensional simu- lations by Klaassen & Peltier (1985a) show that increasing the Prandtl num- ber from 0.5 to 2 significantly influences the energy budget of the primary KH instability, and also makes the two-dimensional billow more unstable to secondary convective instabilities inside the core of the billow. In numerical simulations of Smyth (2003), increasing the Prandtl number from 1 to 7 en- hanced the instability of the braid to secondary shear instabilities. Cortesi et al. (1998, 1999) and Staquet (2000) point out that the higher baroclinic generation of vorticity at higher Prandtl numbers increases the susceptibility of the two-dimensional KH billow to three-dimensional instabilities and en- hances the three-dimensional turbulent motions. However, there has been less agreement on the effect of Prandtl number on mixing in previous studies (in all these studies the Reynolds number is close to or above the mixing tran- sition). The results of Staquet (2000) suggests the increase in the Prandtl number from 0.7 to 1.4 has an inconsistent effect on mixing at different levels of density stratification. Cortesi et al. (1999) discuss that the enhanced three- dimensionality at higher Prandtl numbers, by increasing the Prandtl number from 0.005 to 2.2, results in higher entrainment and mixing rates. However, Smyth et al. (2005) show that in a double diffusive KH billow the turbulent breaking of the billow induces less mixing in the scalar field with the Prandtl 91 4.2. Numerical simulation number of 50 compared to the field with the Prandtl number of 7. Here, our objective is to examine the effect of Prandtl number on the amount of mixing for Reynolds numbers below and within the mixing tran- sition using direct numerical simulations. Our investigation includes Prandtl numbers as high as 64 in three-dimensional KH billows and 700 in two-dimensional KH billows. To achieve numerical simulation of such high Prandtl numbers we were limited to Reynolds numbers that cover the lower and middle range of the mixing transition based on our previous mixing study (Rahmani et al., 2011a). The paper is organized as follows. Section 2 presents the details of our numerical simulations. The two and three-dimensional evolution of the KH billow is described in sections 3 and 4. Time evolution of energy parti- tions and mixing are discussed in section 5. The overall amount of mixing is examined in section 6, and the conclusions stated in section 7. 4.2 Numerical simulation As in Rahmani et al. (2011a), we choose an initial flow with density and horizontal velocity profiles defined by U(z) = ∆U 2 tanh( 2z δ0 ), and ρ̄(z) = −∆ρ 2 tanh( 2z h0 ), (4.1) with δ0 being the vorticity interface thickness , h0 the density interface thick- ness, ∆U the horizontal velocity difference and ∆ρ the density difference. The dimensionless numbers for this flow are Re0 = δ0∆U ν , J = ∆ρgδ0 ρ0∆U2 , P r = ν κ , and R = δ0 h0 , (4.2) the initial Reynolds number, bulk Richardson number, Prandtl number and scale ratio, where ν denotes the kinematic viscosity, κ the molecular diffusivity, ρ0 the reference density and g the gravitational acceleration. The scale ratio is related to the Prandtl number by R = Pr1/2. 92 4.3. Two-dimensional evolution of KH billows The Boussinesq Navier-Stokes equations are solved using the model de- scribed by Winters et al. (2004) for two-dimensional simulations and Smyth et al. (2005) for three-dimensional simulations. The choice of simulation pa- rameters is presented in table 4.1. In these simulations Re0 = 100 and 300 are below the mixing transition, Re0 = 400, 600 and 900 are in the mixing transition range, and Re0 = 1200 is above the mixing transition (Rahmani et al., 2011a). In all simulations one wavelength of the fastest growing mode of the instability, Lx, is considered. The height of the domain is Lz = 9δ0, and the width of the domain in three-dimensional simulations is Ly = Lx/2. Initial and boundary conditions are the same as the ones used in Rahmani et al. (2011a) for three-dimensional simulations and Rahmani et al. (2011b) for two-dimensional simulations. The streamwise, spanwise and vertical num- ber of grid points are Nx, Ny and Nz. In all simulations the grid spacing is smaller than 2LB, where LB is the Batchelor length scale (Batchelor, 1959). 4.3 Two-dimensional evolution of KH billows Snapshots of the time evolution of the KH billow from the two-dimensional simulations are shown in figure 4.1 for Pr = 1, 9, 64 and 700, and Re0 = 300. This Reynolds number is below the mixing transition according to our earlier investigation, Rahmani et al. (2011a), for Pr = 1. At Re0 = 300 and Pr = 1, diffusion dominates and the two-dimensional billow does not develop. As Pr rises, the billow is structured since diffusion is less important. The dimensionless times, t∗ = t∆U/δ0, presented in figure 4.1 belong to a phase after the billow reaches its maximum height and before the time of growth of any potential three-dimensional instabilities (as will be seen in the three- dimensional simulations). As Pr increases, the primary billow is less diffuse and mixed areas of fluid generated until t∗ = 100 reduce significantly. The two high Pr billows exhibit some density structure features that are absent in the density structure of 93 4.3. Two-dimensional evolution of KH billows Simulation dimensions Pr R Lx/δ0 Nx ×Ny ×Nz Re0 = 100 1 3D 1 1 8.3 128×64×128 2 3D 9 3 9.1 128×64×128 3 3D 16 4 9.5 128×64×128 4 3D 25 5 9.8 128×64×128 Re0 = 300 5 3D 1 1 7.5 256×128×256 6 3D 9 3 8.0 384×192×384 7 3D 16 4 8.6 384×192×384 8 3D 25 5 8.8 512×256×512 9 3D 64 8 9.1 512×256×512 10 2D 700 26.5 9.8 1600×0×1600 Re0 = 400 11 3D 1 1 7.4 256×128×320 12 3D 9 3 8.2 384×192×384 13 3D 16 4 8.5 384×192×384 14 3D 25 5 8.7 640×320×640 Re0 = 600 15 3D 1 1 7.3 256×128×320 16 3D 9 3 8.1 480×240×480 17 3D 16 4 8.4 512×256×512 18 3D 25 5 8.6 640×320×640 Re0 = 900 19 3D 1 1 7.2 256×128×320 20 3D 9 3 8.0 512×256×512 Re0 = 1200 21 3D 1 1 7.2 256×128×320 22 3D 9 3 8.0 640×320×640 Table 4.1: Description of the two and three-dimensional simulations carried out. In all these simulations J = 0.03, Lz = 9δ0, and in three-dimensional simulations Ly = Lx/2. The presented mesh sizes are for the density field, the velocity field mesh sizes are half of these. Simulations 5 to 9, and 11 to 14 have also been performed in two dimensions. 94 4.3. Two-dimensional evolution of KH billows Figure 4.1: Snapshots of the density structure of the billow at Re0 = 300 at different non-dimensional times t∗ = 50, 70, 80 and 100 from left to right and for Pr = 1, 9, 64 and 700 from top to bottom. 95 4.3. Two-dimensional evolution of KH billows the two low Pr billows due to the rapid diffusion. At Pr = 1 and 9 the density structure of the billow has been mainly destroyed by diffusion before t∗ = 100. The highly structured density fields of high Pr billows resemble the density field of the billows observed in the experiments of Schowalter et al. (1994) with Pr ∼ 1500 and Atsavapranee & Gharib (1997) with Pr ∼ 600. Among the features common to our high Pr simulations and the mentioned experimental studies are: oscillations in the size of the billow, thinning of the top rolling layer, projection of a small wave from the left corner of the core and development of small-scale distortions inside the billow. We discuss these features in the following. The top rolling layer is stretched by the strain field and thinning continues until it is in balance with the opposing effect of diffusion. This layer forms a very small length scale in the billow, which is close to the Corcos-Sherman length scale (Rahmani et al., 2011b). This thinning process almost results in a detachment of the core fluid from the bottom layer fluid in the Pr = 64 and 700 billows at t∗3 = 70. At this time, a small amplitude wave is gener- ated in the corners of the billow near the braid center. We speculate that this wave projection is caused by the interaction of baroclinically generated vorticity layers of opposite signs. This wave subsequently plays an influen- tial role in the evolution of the billow; after growing and getting entrained by the large scale vortex, the wave replaces the previously existing braid and entrains more fluid inside the billow (t∗ = 80). After this stage, the inter- action of highly concentrated baroclinic vorticity at sharp density interfaces generates small-scale distortions in the flow (t∗ = 100). These instabilities create thin stretched dissipating layers that can enhance mixing. Oscillations of the billow that are caused by oscillations of the core elliptical vortex (Miura & Sato, 1978), are more pronounced at higher Pr. The stronger core vortex at higher Pr (Rahmani et al., 2011b) is the reason for stronger oscillations. These oscillations make the billow periodically reach its maximum and mini- mum horizontal extent, e.g. at t∗ = 70 and 80, respectively. As will be seen in 96 4.4. Three-dimensional evolution of the billow section 5, the stronger oscillations at higher Pr cause higher rates of energy exchange between the instability and the mean flow. As suggested by Klaassen & Peltier (1985b,a, 1991); Caulfield & Peltier (2000); Peltier & Caulfield (2003) the structure of the two-dimensional flow in the pre-turbulent stage determines the nature of secondary three-dimensional instabilities and the transition to turbulence. Therefore, the significantly dif- ferent billow structures in the pre-turbulent stage for different Pr indicates different transition mechanisms to turbulence. Specifically, the small-scale disorganized motions in the two-dimensional flow at high Pr can play a sig- nificant role in enhancing turbulent motions in the flow. We investigate this in the following section. 4.4 Three-dimensional evolution of the billow As suggested by previous studies (e.g. Caulfield & Peltier, 2000), we divide the three-dimensional evolution of a KH billow into four phases. The first phase is the two-dimensional growth of the billow ending when a maximum height is reached. The second phase is the evolution of the billow without growing in two dimensions ending with the growth of secondary three-dimensional instabili- ties. The third phase is the growth of three-dimensional motions until reach- ing a maximum intensity. The fourth phase is the decay of three-dimensional motions and destruction of the billow until a laminar flow is reached. The quantitative definition of these phases is given by Rahmani et al. (2011a). In figure 4.2 we compare the density structure of these three-dimensional structures at their maximum intensity for Re0 = 300 and 400, for increas- ing Pr. At Re0 = 300 and Pr = 1, three-dimensional motions are mainly suppressed by viscosity and therefore the billow merely diffuses out without becoming turbulent. Because of this the overall amount of mixing for this case is considerably lower than billows with higher Reynolds numbers that strongly develop three-dimensional instabilities. However, as Pr increases for Re0 = 97 4.4. Three-dimensional evolution of the billow 300, the three-dimensional motions become stronger and more structured. So as suggested by the two-dimensional simulations in section 3, the more struc- tured density of the pre-turbulent billows at higher Pr and the small-scale distortions generated in the flow enhance the three-dimensionality. This effect is also pointed out by Cortesi et al. (1999); Staquet (2000), and attributed to higher baroclinically generated vorticity at sharper density gradients at higher Pr. This Prandtl number effect is strong enough to overcome viscous effects and sustain three-dimensional structures through the life-cycle of the billow even at very low Reynolds number of Re0 = 300. These motions increase the vertical extent of the billow, and also provide more interfacial area between the scalars for mixing. Hence, we expect higher overall amount of mixing at higher Pr for Re0 = 300. We will investigate this in section 5. At Re0 = 400, the Reynolds number is high enough for the flow to de- velop three-dimensional instabilities even at Pr = 1. By the time shown in figure 4.2, the counter-rotating streamwise vortices (Lin & Corcos, 1984) have rolled up the fluid into a mushroom-like structure and increased the vertical extent of the billow at all Pr. As Pr increases, the density gradients become sharper and smaller scale structures develop in the flow. However, this does not enhance the three-dimensionality and vertical extent of the flow as it did at Re0 = 300. This is because the less diffuse density structure at higher Pr is more strongly stratified and therefore requires greater turbulent energy to be vertically displaced. Hence, when three-dimensionality is strong, the increase in Pr does not act in favor of entrainment, although it locally enhances the intensity of three-dimensional motions. So, for Reynolds number within or above the mixing transition we expect less overall mixing for higher Pr. We will examine mixing and energy partitions quantitatively in the next section. 98 4.5. Mixing and energy partitions y/h0 z/ h 0 0 1 2 3 0 2 4 6 8 y/h0 0 1 2 3 y/h0 0 2 4 y/h0 0 2 4 y/h0 0 2 4 y/h0 z/ h 0 0 1 2 3 0 2 4 6 8 y/h0 0 2 4 y/h0 0 2 4 y/h0 0 2 4 Figure 4.2: Snapshots of the span-wise density structure in the middle of the domain at maximum intensity of turbulent motions for Re0 = 300 (top row) and Pr = 1, 9, 16, 25 and 64, from left to right, and Re0 = 400 (bottom row) and Pr = 1, 9, 16 and 25, from left to right. 4.5 Mixing and energy partitions We quantify mixing, M , as the amount of increase in the background poten- tial energy, Pb, where Pb is the minimum potential energy obtained when fluid particles are rearranged adiabatically to a stable state (Winters et al., 1995). This quantification of mixing and also our other energy partitions exclude dif- fusion of the initial background stratification. When studying mixing, it is also 99 4.5. Mixing and energy partitions useful to examine other energy partitions: the available potential energy, Pa, the two-dimensional turbulent kinetic energy, K2d, and the three-dimensional turbulent kinetic energy, K3d (Winters et al., 1995; Caulfield & Peltier, 2000). We discuss the time evolution of these energy partitions through the life-cycle of the billow for Re0 = 300 and 400. Re0 = 300: Figure 4.3 shows the time evolution of the energy partitions, normalized by the initial kinetic energy, K0, and added together so that the top line shows the total energy extracted from the mean flow, for Re0 = 300 and increasing Pr. As Pr increases, the total energy extracted by the instability during the two-dimensional phases increases and higher oscillations occur in K2d and P . This was predictable from the density structures discussed in section 3, which suggested higher shear production of turbulent kinetic energy and stronger oscillations of the core vortex at higher Pr. The significant influence of Pr on energy partitions of KH billows was previously pointed out by Klaassen & Peltier (1985a). At Pr = 700, the total energy of the billow at its peak is almost two times greater than that at Pr = 1. However, the mixing, i.e. the increase in Pb, during the two-dimensional phases is small, slow and limited by the diffusion of small scales at high Pr. During the first three phases, a large portion of the the total extracted energy is in the form of Pa for large Pr. An important effect of increasing Pr is enhancing the three-dimensional motions, seen as the increase in K3d for higher Pr in figure 4.3. This effect is more pronounced for Pr > 9, and gives rise to the total energy extracted by the instability in the three-dimensional phases. As the beginning of the mix- ing transition is marked by the sustainability of three-dimensional motions, this effect of Pr on the generation of three-dimensional motions at low Re0 suggests that the transition in mixing can occur at lower Reynolds numbers for higher Pr. However, the three-dimensional motions generated at high Pr have small influence on the entrainment and mixing. Firstly, this is because 100 4.5. Mixing and energy partitions 0 0.05 0.1 Pr = 1 Pr = 9 0 0.05 0.1 Pr = 16 Pr = 25 0 100 200 300 0 0.05 0.1 t* Pr = 64 0 100 200 300 t* Pr = 700 (2D) Figure 4.3: Time evolution of the energy partitions normalized by K0 for Re0 = 300. The lines from bottom to top show the background potential energy, Pb, the available potential, Pa, the two-dimensional turbulent kinetic energy, K2d, and the three-dimensional turbulent kinetic energy, K3d partitions added cumulatively together. Pa and K3d partitions are shaded. The vertical lines show the transition times between different phases. The Pr = 700 energy partitions are from a two-dimensional simulation. 101 4.5. Mixing and energy partitions the effect of Pr on increasing the energy of three-dimensional motions is much weaker than the effect of Re0 (as seen in Rahmani et al., 2011a). Secondly, due to the very slow mixing at high Pr a portion of the available potential energy returns to the turbulent kinetic energy and finally dissipated, without being converted to mixing. Hence, the small increase in mixing with increasing Pr is mainly due to the greater entrainment during the two-dimensional phases, and less significantly due to the entrainment generated by three-dimensional motions. This indicates that because of the small Re0, the two-dimensional roll-up is still the main mechanism for the majority of mixing, however this mixing is delayed until three-dimensional phases for high Pr. Re0 = 400: Time evolution of energy partitions for Re0 = 400 is shown in figure 4.4. This Reynolds number is sufficiently high to impose strong three-dimensional motions. The effects of Pr during the two-dimensional phases are similar to those at Re0 = 300: larger P and K2d, stronger oscillations and lower mixing as Pr rises. Similar to Re0 = 300, the increase in Pr prolongs the life-cycle of the billow, and delays the mixing caused by the two-dimensional roll-up until third and fourth phases. However, during the three-dimensional phases the increase in Pr does not enhance the three-dimensional motions. Espe- cially at Pr = 25, K3d becomes noticeably smaller, with less entrainment and mixing generated in the turbulent phase. As explained in section 4, the less diffuse density field at higher Pr suppresses buoyancy fluxes and therefore the entrainment in three-dimensional phases. This effect was also seen in experi- ments of Turner (1968). However, the stronger two-dimensional entrainment at higher Pr compensates for this effect and the total amount of entrainment and mixing only weakly drops as Pr rises. 102 4.5. Mixing and energy partitions 0 0.05 0.1 Pr = 1 Pr = 9 0 100 200 300 400 0 0.05 0.1 t* Pr = 16 0 100 200 300 400 t* Pr = 25 Figure 4.4: Time evolution of the energy partitions normalized by K0 for Re0 = 400. See the caption of figure 4.3 for the rest of the descriptions. Diffusion times scales: Figures 4.3 and 4.4 suggest that as Pr increases the flow remains unmixed for a longer period of time. We define a dimensionless time, t∗M , after which mixing starts to grow in a real sense, and quantitatively the first time M/K0 exceeds 10−3. Figure 4.5 shows variation of t∗M with Pr at different Re0. This figure shows that t∗M increases with increasing Pr and Re0. An appropriate measure of diffusion time scale for the background flow is δ20/κ. Therefore, the diffusion time scale of the backround flow is proportional to Pr. By increas- ing Pr from 1 to 700, t∗M becomes almost 5 times larger. This indicates that compared to the molecular diffusion time scale, t∗M grows much slower as the roll-up by the billow hastens mixing. 103 4.5. Mixing and energy partitions 1 10 100 100020 40 60 80 100 120 Pr t* M Figure 4.5: Variation of the dimensionless dissipation time, tD, with Pr for Re0 = 300, 400, and 600 from the lowest to highest line. Two-dimensional versus three-dimensional mixing: Figure 4.6 shows the rate of mixing, φM , from two and three-dimensional simulations at Re0 = 300 for different Pr. At Re0 = 300, KH billows simulated in three dimensions generated slightly higher mixing during phases 3 and 4 compared to KH billows simulated in only two dimensions. The difference was noticeable only at Pr = 16 and 25. As Re0 = 300 is below the mixing transition for Pr = 1 and 9, the effects of three-dimensionality on mixing are not significant. Therefore, two and three-dimensional mixing are relatively close throughout the entire life-cycle of the billow. At Re0 = 400, two and three-dimensional simulations induced significantly different mixing, as shown in figure 4.7. The two and three-dimensional φM start to diverge in phase 3. As Pr increases, the two-dimensional mixing occurs over longer time due to the slower mixing of the rolled-up billow. As will be seen in the next section, the total amount of the two-dimensional mixing, the 104 4.6. Overall mixing area below the φM curve, slightly rises for higher Pr. The mixing caused by three-dimensional motions diminishes as Pr rises. Figure 4.7 shows how significantly the turbulent mixing is suppressed by the increase in Pr. While the billow lives long enough for the rolled-up interface to completely mix, the duration of turbulence is not long enough for high Pr turbulent mixing. 4.6 Overall mixing Figure 4.8 presents the Pr dependence of the overall mixing induced over the lifetime of KH billows for different Re0. At Re0 = 100, the final amount of mixing does not depend on Pr as viscous effects are dominant. Just before the mixing transition, i.e. Re0 = 300, the final amount of mixing increases with increasing Pr. There is a transition-like behavior in mixing at Re0 = 300 for Pr > 9, for which small-scale features become significant (figures 4.1 and 4.2). The trend of mixing dependence on Pr at Re0 = 300 and the amount of mixing is to a great extent predictable from two-dimensional simulations, especially for Pr = 1 and 9. For Re0 > 300, while the two-dimensional amount of mixing slightly in- creases with Pr, the overall amount of mixing decreases as Pr increases. The two-dimensional mixing has a weak dependence on Pr as the entrained fluid by the two-dimensional roll-up gets completely mixed inside the region for these cases. However, Pr has strong effects on the entrainment and mixing rate of three-dimensional billows. Reduced overall mixing for higher Pr at Reynolds numbers above the mixing transition is in agreement with the re- sults of Breidenthal (1981); Konrad (1976); Smyth et al. (2005). Figure 4.8 suggests that both below and during the mixing transition the effect of Re0 on the final amount of mixing is more significant than the effect of Pr. This is be- cause of the dominant effect of Re0 on the two-dimensional roll-up before the transition, and on the three-dimensional turbulent motions after the transition. 105 4.6. Overall mixing 0 100 200 300 400 0 0.2 0.4 0.6 0.8 1 x 10−4 t* φ M Pr = 16 0 100 200 300 400 t* Pr = 25 Figure 4.6: Time variation of the mixing rate, φM , from two-dimensional sim- ulations (lower line) and three-dimensional simulations (upper line) at Re0 = 300. The filled areas show the three-dimensional mixing. The vertical dashed lines indicate the transition times between the phases. The results for Pr = 1 and 9 are not plotted since the two- and three-dimensional simulations are indistinguishable. 0 0.2 0.4 0.6 0.8 1 x 10−4 φ M Pr = 1 Pr = 9 0 100 200 300 400 0 0.2 0.4 0.6 0.8 1 x 10−4 t* φ M Pr = 16 0 100 200 300 400 t* Pr = 25 Figure 4.7: Time variation of the mixing rates at Re0 = 400. See the caption of figure 4.6 for details. 106 4.6. Overall mixing 100 101 102 0 0.02 0.04 0.06 0.08 100 300 400 600 Pr M /K 0 Figure 4.8: Variation of the total amount of mixing with the Prandtl number. The numbers below each solid line show the initial Reynolds number. The dashed line shows mixing from two-dimensional simulations at Re0 = 300 and the dashed-dotted line mixing from two-dimensional simulations at Re0 = 400. Effects of Pr on mixing transition: Previous mixing measurements of Konrad (1976) and Breidenthal (1981) suggest that the increase in Pr by almost three orders of magnitude does not affect the transitional Reynolds numbers, but reduces the amount of mixing for all Reynolds numbers. In figure 4.9 the mixing transition found by Rahmani et al. (2011a) at Pr = 1 is compared to the trend of mixing transition at Pr = 9. The transitional Re0 are the same for both Pr = 1 and 9. The amount of mixing is slightly lower at Pr = 9 for Re0 > 300. So by increasing Pr by almost one order of magnitude the trend of mixing transition does not change, for Reynolds number within and above the mixing transition the overall amount of mixing drops slightly, and for Reynolds numbers below the mixing transition the overall amount of mixing almost does not change. The first two results are in agreement with the previous laboratory observations, 107 4.7. Summary and conclusions 102 103 104 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Re0 M /K 0 Figure 4.9: Variation of the total amount of mixing gained at the end of the KH life-cycle, M , normalized by the initial kinetic energy, K0, for Pr = 1 (solid line) and Pr = 9 (dashed line). while the third result is different. 4.7 Summary and conclusions Two and three-dimensional DNS were performed to examine the effect of Pr, varying over a wider range compared to previous studies, on the evolution and overall mixing in KH billows. During the two-dimensional evolution of the billow, the increase in Pr generated a more structured density field, with small- scale features that triggered secondary three-dimensional instabilities even at a low Reynolds number of 300. This effect modified the three-dimensional structure of the flow and the amount of mixing, especially for Pr > 9. The time evolution of the billow and its energy budgets are highly depen- dent on Pr at all Re0. During the two-dimensional phases, a higher Pr gener- ates larger billows and energy partitions. During the three-dimensional phases, 108 4.7. Summary and conclusions a higher Pr locally enhances the turbulent motions and entrainment, but re- duces the entrainment by large eddies caused by strong three-dimensional in- stabilities. Mixing is always slower for higher Pr. The Pr dependence of the final amount of mixing shows two distinct be- haviors for Reynolds numbers below and above the beginning of the mixing transition. For Re0 = 300, below the mixing transition, the overall mixing increases as Pr rises. For Re0 > 300, the increase in Pr reduces the overall mixing. The latter trend is in agreement with the results of Konrad (1976); Breidenthal (1981) in mixing layers, and Smyth et al. (2005) in the double diffusion context. In general, the the results of this chapter suggest that, for the range of Reynolds numbers considered, the following effects of the Prandtl number are important for mixing: • Turbulent kinetic energy and entrainment: The increase in Pr has sig- nificant implications for the strength of turbulent three-dimensional mo- tions, buoyancy fluxes, and therefore for mixing. These effects are through the changes in the density structure of the flow. However, two-dimensional entrainment is not significantly dependent on Pr. • Mixing time scales: The increase in Pr changes the lifetime and diffusion time scales of the flow, however not in a predictable manner. The two- dimensional mixing occurs over longer time scales for slower molecular mixing at higher Pr, and therefore the overall two-dimensional mix- ing is weakly dependent on Pr. However, the three-dimensional mixing occurs over the limited time of existence of turbulence. Therefore, three- dimensional mixing is repressed by slower molecular mixing at high Pr. 109 Chapter 5 Conclusions 5.1 Summary Different aspects of the nonlinear evolution of Kelvin-Helmholtz (KH) insta- bilities in sheared density stratified flows have been studied using direct nu- merical simulations (DNS). The growth of KH instabilities in initially stably stratified sheared flow precedes the development of small-scale chaotic motions and turbulence. As described in chapter 3, the entire life-cycle of a KH in- stability driven flow consists of four phases: two-dimensional growth of the primary KH instability, pre-turbulent evolution of the KH billow, growth of turbulence, and decay of turbulence. In real geophysical flows at high Reynolds numbers (Re0 & 10 5) the growth and decay of turbulence in KH billows comprise a significant portion of the life- time of the billows (e.g. Seim & Gregg, 1994). Mixing generated by small-scale motions during these phases significantly enhances the background potential energy of the stratification at the cost of the kinetic energy of the background shear. Through this process also a portion of the energy extracted from the background shear converts to internal energy due to the viscous dissipation. While the ultimate goal of studies of shear instabilities is to characterize the mixing generated through turbulent phases, it is also essential to understand the properties of the primary instabilities in the pre-turbulent phases. This is because the transition to turbulence and characteristics of the turbulent flow is greatly under the influence of the primary KH instability (e.g. Dimotakis & Brown, 1976). In chapter 2, the structure of the primary KH instability was examined 110 5.2. Main contributions in detail. This study mainly focuses on the properties of the flow during the first phase of the evolution of KH billow. By the end of phase 1, the primary billow has reached its maximum height and also generated a very thin density interface on the braid. So the large and small scales of the two-dimensional flow field are well described at this stage. Chapters 3 and 4 of the thesis investigated the evolution of the KH billow through all four phases. In chapter 3, the effect of Reynolds number on mixing is assessed for a single Prandtl number of Pr = 1. In chapter 4, the effect of Prandtl number on mixing is studied at various Reynolds numbers. 5.2 Main contributions The two-dimensional evolution of the vorticity and density fields of KH in- stabilities were examined in chapter 2. This study was assisted by comparing the DNS results to the predictions of the model of Corcos & Sherman (1976). The effect of Prandtl number over a wide range of 1 < Pr < 700, on the evolution of the primary billow was studied for the first time. Previous in- vestigation of the effect of Reynolds number on the evolution of the primary billow by Patnaik et al. (1976); Staquet (1995); Smyth (2003) was extended to include Reynolds numbers as high as Re0 = 10 5. The results of this chapter showed that the evolution of the primary KH billow depends weakly on Re0 for sufficiently high Re0, but remains dependent on Pr for all Prandtl num- bers. The improvements in the predictions of the model of Corcos & Sherman (1976) at high Re0 and Pr makes this model an ideal tool for predicting some characteristics of KH billows such as their vertical extent, braid thickness and susceptibility to braid shear instabilities in real oceanic and atmospheric flows. Based on the predictions of the model of Corcos & Sherman (1976) for the braid thickness, and also suggestions by Lin & Corcos (1984) and Corcos & Lin (1984), a new Corcos-Sherman length scale was introduced in chapter 2. This length scale represents the smallest scale of motion in the two-dimensional 111 5.2. Main contributions KH flow field. The Corcos-Sherman length scale was compared to different important small length scales in the three-dimensional turbulent flow field. It was shown that a small fraction of the total dissipation of the density field is associated with the length scales smaller than the Corcos-Sherman scale. Hence, the Corcos-Sherman length scale is also a useful length scale in the turbulent three-dimensional KH flow field. The significance of this finding is that a length easily predictable from the model of Corcos & Sherman (1976) can be used as a guide for capturing small scales in numerical simulations, or laboratory and field measurements of KH billows. The effect of Reynolds number on mixing was investigated in chapter 3, for 100 < Re0 < 4800. The transition in mixing by increasing the Reynolds number was found for the first time using numerical simulations. This mixing transition manifested itself by an almost one order of magnitude increase in the overall amount of mixing generated through the entire life-cycle of a KH billow over a range of outer-scale Reynolds numbers of 900 . Reδ . 5000. This transitional range is within the transitional range of outer-scale Reynolds numbers of 103 . Reδ . 10 4 found by Breidenthal (1981) and Koochesfahani & Dimotakis (1986) in the laboratory. The mixing transition in DNS was ex- plained by studying the time variation of the rate mixing and energy partitions through the evolution of the billow during different phases. While it is known that the rate and characteristics of mixing vary significantly through the evo- lution of KH billow (e.g. Patterson et al., 2006; Smyth et al., 2007), numerical simulations of the time evolution of these quantities through the entire life- cycle of the billow have been very limited in the literature. The investigation in chapter 3 is one of the few mixing studies that covers the lifetime of KH billows (e.g. see Inoue & Smyth, 2009). Chapter 4 described the effect of Prandtl number on the evolution of KH billow, the amount of overall mixing and the time variation of energy parti- tions. The two-dimensional evolution of the primary KH billow was studied for a range of Prandtl numbers of 1 < Pr < 700. This study showed for the 112 5.3. Future research first time that the increase in the Prandtl number results in important changes in the density structure of the flow. These include small-scale features at high Prandtl numbers that significantly influence the mechanism of transition to turbulence. Understanding these effects is especially important because in numerical simulations usually the results of low Pr simulations are extended to high Pr flows. In three dimensions, the effects of Prandtl number in the range of 1 < Pr < 64 on the amount of mixing and time evolution of energy partitions were examined at different Reynolds numbers. This investigation is the most complete set of simulations performed with the goal of exploring the Re0-Pr plane for the overall amount of mixing. The results of this chapter demonstrated that depending on the Reynolds number the increase in Pr re- duces or increases the overall amount of mixing. The effect of Prandtl number on mixing is more pronounced for Pr > 9. Results of chapters 3 and 4 suggest that the dependence of mixing on Reynolds number is strong within the transitional Reynolds numbers and weak outside this range. Hence, for geophysical flows with sufficiently high Reynolds numbers it can be assumed that mixing does not change significantly with the Reynolds number. However, except for the Reynolds numbers below the mixing transition, the increase in Pr showed a monotonic, but gentle decrease in the overall amount of mixing. The effects of Pr are more important when the time evolution of mixing is of interest. The increase in Re0 and Pr both prolong the lifetime of the billow. 5.3 Future research Typical to all DNS studies, the numerical simulations in this thesis have been limited by the overwhelming cost of computations. As a result we did not have all the freedom we desired in choosing the values of Re0 and Pr. Considering the continuous progress in available computational resources, we recommend the following research for future work. While the variation of mixing in KH bil- 113 5.3. Future research lows with the Reynolds number exhibited a clear monotonic increase through a transition in chapter 3, the post-transition non-monotonic behaviour of mix- ing still needs further investigation for Re0 > 4800. Simulations at high Pr were restricted to low Re0 (e.g. Re0 6 600 at Pr = 25). Therefore, a com- prehensive mapping of the Re0-Pr plane for the overall amount of mixing was not possible. Vortex pairing is an important feature of KH billows, which was not inves- tigated in this thesis. Vortex pairing can have significant implications for the mixing transition (Moser & Rogers, 1991). Extra stretching of the braid and merging of two billow cores due to the vortex pairing influences the suscepti- bility of the braid and core to secondary instabilities and therefore transition to turbulence. In addition to changing the mixing behaviour, different types of secondary instabilities that appear in the flow because of the vortex pairing may introdue new small scales to the flow field. Under this condition, the usefulness of the Corcos-Sherman length scale, discussed in chapter 2, should be assessed. The results of this paper describe the evolution and mixing properties of KH billows for a bulk Richardson number of J = 0.03 (and 0.035 for some simulations). As J increases the instability and mixing are suppressed by the stronger stratification (e.g. Atsavapranee & Gharib, 1997). Geophysical flows exhibit different values of J depending on the strength of stratification. Understanding the effects of J on mixing transition, small scales, and complex- ities that high Prandtl numbers add to the flow is an important part of shear instabilities research. Especially, there is a transition from KH to Holmboe waves (Holmboe, 1962) at sufficiently high values of J . Studying the effects of Prandtl and Reynolds number on Holmboe wave field is an interesting topic for future research. 114 Bibliography Alexakis, A. 2005 On Holmboe’s instability for smooth shear and density profiles. Phys. Fluids 17, 84–103. Atsavapranee, P. & Gharib, M. 1997 Structures in stratied plane mix- ing layers and the efects of cross-shear small-scale variation of convected quantities like temperature in turbulent fluid. J. Fluid Mech. 342, 53–86. Baines, P. & Mitsudera, H. 1994 On the mechanism of shear flow insta- bilities. J. Fluid Mech. 276, 327–342. Batchelor, G. K. 1959 Small-scale variation of convected quantities like temperature in turbulent fluid. J. Fluid Mech. 5, 113–133. Batchelor, G. K. 1972 An introduction to fluid dynamics. Cambridge Uni- versity Press. Bernal, L. P. & Roshko, A. 1986 Streamwise vortex structure in plane mixing layers. J. Fluid Mech. 170, 499–525. Breidenthal, R. 1981 Structure in turbulent mixing layers and wakes using a chemical reaction. J. Fluid Mech. 109, 1–24. Broadwell, J. E. & Breidenthal, R. E. 1982 A simple model of mixing and chemical reaction in a turbulent shear layer. J. Fluid Mech. 125, 397– 410. Broadwell, J. E. & Mungal, M. G. 1990 Large-scale structures and molecular mixing. Phys. Fluids A 3(5), 1193–1206. Brown, G. L. & Roshko, A. 1974 On density effescts and large structure in turbulent mixing layers. J. Fluid Mech. 64, 775–816. Browning, K. A. 1971 Structure of the atmosphere in the vicinity of large- amplitude kelvin-helmholtz billows. Q. J. R. Meteorol. Soc. 97, 283–299. 115 Bibliography Buch, K. A. & Dahm, W. J. A. 1996 Experimental study of the fine-scale structure of conserved scalar mixing in turbulent shear flows. part 1. Sc > 1. J. Fluid Mech. 317, 21–71. Buch, K. A. & Dahm, W. J. A. 1998 Experimental study of the finne-scale structure of conserved scalar mixing in turbulent shear flows. part 2. Sc ≈ 1. J. Fluid Mech. 364, 1–29. Carpenter, J. R., Balmforth, N. J. & Lawrence, G. A. 2010a Iden- tifying unstable modes in stratied shear layers. Phys. Fluids 22, 054104–1– 054104–13. Carpenter, J. R., Lawrence, G. A. & Smyth, W. D. 2007 Evolution and mixing of asymmetric Holmboe instabilities. J. Fluid Mech. 582, 103– 132. Carpenter, J. R., Tedford, E. W., Rahmani, M. & Lawrence, G. A. 2010b Holmboe wave fields in simulation and experiment. J. Fluid Mech. 648, 205–223. del Castillo-Negrete, Diego 2000 Self-consistent chaotic transport in fluids and plasma. Chaos 10, 75–88. Caulfield, C. P. & Peltier, W. R. 1994 Three dimensionalization of the stratified mixing layer. Phys. Fluids 6A, 3803–3805. Caulfield, C. P. & Peltier, W. R. 2000 The anatomy of the mixing transition in homogeneous and stratified free shear layers. J. Fluid Mech. 413, 1–47. Corcos, G. M. & Lin, S. J. 1984 The mixing layer: deterministic models of a turbulent flow: Part 2. the origin of the three-dimensional motion. J. Fluid Mech. 139, 67–95. Corcos, G. M. & Sherman, F. S. 1976 Vorticity concentration and the dynamics of unstable free shear layers. J. Fluid Mech. 73, 241–264. Corcos, G. M. & Sherman, F. S. 1984 The mixing layer: deterministic models of a turbulent flow. part 1. introduction and the two-dimensional flow. J. Fluid Mech. 139, 29–65. 116 Bibliography Cortesi, A. B, Smith, B. L., Yadigaroglu, G., & Banerjee, S. 1999 Numerical investigation of the entrainment and mixing processes in neutral and stably-stratified mixing layers. Phys. Fluids 11, 162–184. Cortesi, A. B, Yadigaroglu, G., & Banerjee, S. 1998 Numerical in- vestigation of three-dimensional structures in stably-stratified mixing layers. Phys. Fluids 10, 1449–1473. Davidson, P. A. 2005 Turbulence, An introduction for scientists and engi- neers. Oxford University Press. DeSilva, I. P. D., Fernando, H. J. S., Eaton, F. & Hebert, D. 1996 Evolution of Kelvin-Helmholtz billows in nature and laboratory. Earth and Planetary Science Letters 143, 217–231. Dimotakis, P. E. 2000 The mixing transition in turbulent flows. J. Fluid Mech. 409, 69–98. Dimotakis, P. E. 2005 Turbulent mixing. Annu. Rev. Fluid Mech. 37, 329– 356. Dimotakis, P. E. & Brown, G. L. 1976 The mixing layer at high Reynolds number: the large-structure dynamics and entrainment. J. Fluid Mech. 78, 535–560. Drazin, P. G. & Reid, W. H. 1981 Hydrodynamic stability . Cambridge University Press. Farmer, D. & Armi, L. 1999 Stratified flow over topography: The role of small-scale entrainment and mixing in flow establishment. Proc. Roy. Soc. London 455A, 3221–3258. Gargett, A., Merryfield, W. & Holloway, G. 2003 Different numer- ical simulation of differential scalar diffusion in three-dimensional stratified turbulence. J. Phys. Oceanogr. 33, 1758–1782. Geyer, W. R. & Farmer, D. M. 1989 Tide induced variation of the dy- namcis of a salt wedge estuary. J. Phys. Oceanogr. 19, 1060–1072. Geyer, W. R., Lavery, A. C., Scully, M. E. & Trowbridge, J. H. 2010 Mixing by shear instability at high Reynolds number. Geophys. Res. Lett. 37, L22607,doi:10.1029/2010GL045272. 117 Bibliography Geyer, W. R. & Smith, J. D. 1987 Shear instability in a highly stratified estuary. J. Phys. Oceanogr. 17, 1668–1679. Goldstein, S. 1931 On the stability of superposed streams of fluids of dif- ferent densities. Proc. R. Soc. Lond. 132, 524–548. Gossard, E. E., Richter, J. H. & Atlas, D. 1970 Internal waves in the atmosphere from high-resolution radar measurements. J. Geophys. Res. 75, 3523.3536,doi:10.1029/JC075i018p03523. Haigh, S. P. 1995 Non-Symmetric Holmboe waves. PhD thesis, the University of British Columbia. Haigh, S. P. & Lawrence, G. A. 1999 Symmetric and non-symmetric Holmboe instabilities in an inviscid flow. Phys. Fluids A 11, 1459–1468. Haren, H. & Gostiaux, L. 2010 A deepocean Kelvin-Helmholtz billow train. Geophys. Res. Lett. 37, L03605,doi:10.1029/2009GL041890. Hazel, P. 1972 Numerical studies of the stability of inviscid parallel shear flows. J. Fluid Mech. 39, 39–61. Herbert, D., Moum, J. N., Paulson, C. A. & Cladwell, D. R. 1992 Turbulence and internal waves at the Equator. part ii: details of a single event. J. Phys. Oceanogr. 22, 1346–1356. Hogg, A. McC. & Ivey, G. N. 2003 The Kelvin-Helmholtz to Holmboe instability transition in stratified exchange flows. J. Fluid Mech. 477, 339– 362. Holmboe, J. 1962 On the behaviour of symmetric waves in stratified shear layers. Geophys. Publ. 24, 67–113. Howard, L. 1961 Note on a paper of John W. Miles. J. Fluid Mech. 10, 509–512. Inoue, R. & Smyth, W. D. 2009 Efficiency of mixing forced by unsteady shear flow. J. Phys. Oceangr. 39, 1150–1166. Ishiharaa, T., Kaneda, Y., Yokokawa, M., Itakura, K. & Uno, A. 2005 Energy spectrum in the near dissipation range of high resolution direct numerical simulation of turbulence. Journal of the Physical Society of Japan 74,5, 1464–1471. 118 Bibliography Ivey, G.N., Winters, K.B. & Koseff, J.R. 2008 Density stratication, turbulence, but how much mixing? Annu. Rev. Fluid Mech. 40, 169–184. Jackson, P. R. & Rehmann, C. R. 2003 Laboratory measurements of differential diffusion in a diffusively stable, turbulent flow. J. Phys. Oceangr. 33, 1592–1603. Jimenez, J. 1983 A spanwise structure in the plane shear layer. J. Fluid Mech. 132, 319–336. Klaassen, G. P. & Peltier, W. R. 1985a The effect of Prandtl number on the evolution and stability of Kelvin-Helmholtz billows. Geophys. Astrophys. Fluid Dynamics 32, 23–60. Klaassen, G. P. & Peltier, W. R. 1985b The onset of turbulence in finite-amplitude Kelvin-Helmholtz billows. J. Fluid Mech. 155, 1–35. Klaassen, G. P. & Peltier, W. R. 1991 The influence of stratification on secondary instability in free shear layers. J. Fluid Mech. 227, 71–106. Konrad, J. H. 1976 An experimental investigation of mixing in two- dimensional turbulent shear ows with applications to diffusion-limited chem- ical reactions. PhD thesis, California Institute of Technology. Koochesfahani, M. M. & Dimotakis, P. E. 1986 Mixing and chemical reaction in a turbulent liquid mixing layer. J. Fluid Mech. 170, 83–112. Koop, C. G. & Browand, F. K. 1979 Instability and turbulence in a stratified layer with shear. J. Fluid Mech. 93, 135–159. Kundu, P. K. & Cohen, I. M. 2004 Fluid Mechanics . Elsevier Academic Press. Kuthnur, P. S. & Clemens, T. 2005 Effects of unsteady strain rates on scalar dissipation structures in turbulent planar jets. Phys. Fluids 17, 125104–125104. Lamb, K. G. & Farmer, D. 2011 Instabilities in an internal solitary-like wave on the Oregon shelf. J. Phys. Oceangr. 41, 67–87. 119 Bibliography Lasheras, J. C . & Choi, H. 1988 Three-dimensional instability of a plane free shear layer: an experimental study of the formation and evolution of streamwise vortices. J. Fluid Mech. 189, 53–86. Lawrence, G. A., Browand, F. K. & Redekopp, L. B. 1991 The sta- bility of a sheared density interface. Phys. Fluids A3, 2360–2370. Lewis, G. S. & Swinney, H. L. 1999 Velocity structure functions, scal- ing, and transitions in high-Reynolds-number Couette-Taylor flow. Physical Review E 59, 5457–5467. Lin, S. J. & Corcos, G. M. 1984 The mixing layer: deterministic models of a turbulent flow. part 3. the effect of plane strain on the dynamics of streamwise vortices. J. Fluid Mech. 141, 139–178. Lorenz, E. N. 1955 Available potential energy and the maintenance of the general circulation. Tellus VII 2, 157–167. Luce, H., Mega, T., Yamamoto, M. K., Yamamoto, M., Hashiguchi, H., Fukao, S., Nishi, N., Tajiri, T., Nakazato, M., Lewis, G. S. & Swinney, H. L. 2010 Observations of Kelvin-Helmholtz instability at a cloud base with the middle and upper atmosphere (MU) and weather radars. J. Geophys. Res. 115, D19116, doi:10.1029/2009JD013519. Ludlam, F. H. 1967 Characteristics of billow clouds and their relation to clear-air turbulence. Q. J. R. Meteorol. Soc. 93, 419–435. Marmorino, G. O. 1987 Observations of small-scale mixing processes in the seasonal thermocline. part ii: wave breaking. J. Phy. Oceanogr. 17, 1348– 1355. Merryfield, W., Holloway, G. & Gargett, A. 1998 Differential verti- cal transport of heat and salt by weak stratified turbulence. Geophys. Res. Lett. 25, 2772–2776. Miles, J. 1961 On the stability of heterogeneous shear flows. J. Fluid Mech. 10, 496–508. Miura, A. & Sato, T. 1978 Theory of vortex nutation and amplitude os- cillation in an inviscid shear instability. J. Fluid Mech. 86, 33–47. 120 Bibliography Moin, P. & Mahesh, K. 1998 Direct numerical simulation: A tool in tur- bulence research. Annu. Rev. Fluid Mech. 30, 539–578. Moser, R. D. & Moin, P. 1987 The effects of curvature in wall-bounded turbulent flows. J. Fluid Mech. 175, 479–510. Moser, R. D. & Rogers, M. M. 1991 Mixing transition and the cascade to small scales in a plane mixing layer. Phys. Fluids A 3, 1128–1134. Moum, J. N., Farmer, D. M., Smyth, W. D., Armi, L. & Vagle, S. 2003 Structure and generation of turbulence at intefcaes strained by solitary waves propagating shoreward over the continental shelf. J. Phys. Oceangr. 33, 2093–2112. Patnaik, P.C., Sherman, F.S. & Corcos, G. M. 1976 A numerical sim- ulation of Kelvin-Helmholtz waves of finite amplitude. J. Fluid Mech. 73, 215–240. Patterson, M. D., Caulfield, C. P., McElwaine, J. N. & Dalziel, S. B. 2006 Time-dependent mixing in stratified Kelvin-Helmholtz billows: Experimental observations. Geophys. Res. Lett. 33, L15608, doi:10.1029/2006GL026949. Pawlak, G. & Armi, L. 1998 Vortex dynamics in a spatially accelerating shear layer. J. Fluid Mech. 376, 1–35. Peltier, W. R. & Caulfield, C. P. 2003 Mixing efficiency in stratified shear flows. Annu. Rev. Fluid Mech. . 35, 135–167. Pierrhumbertt, R. T. & Windall, S. E. 1982 The two- and three- dimensional instabilities of a spatially periodic shear layer. J. Fluid Mech. 114, 59–82. Rahmani, M., Lawrence, G.A. & Seymour, B.R. 2011a The effect of Reynolds number on mixing in Kelvin-Helmholtz instability. Submitted to J. Fluid Mech. -, –. Rahmani, M., Seymour, B.R. & Lawrence, G.A. 2011b Two- dimensional evolution of large and small-scale structures in Kelvin- Helmholtz. Submitted to J. Fluid Mech. -, –. 121 Bibliography Rosenhead, L. 1931 The formation of vortices from a surgace of discontinu- ity. Proc. Roy. Soc. A 134 (823), 170–192. Schowalter, D. G., Atta, C. W. Van & Lasheras, J. C. 1994 A study of streamewise vortex structure in a stratified shear layer. J. Fluid Mech. 281, 247–291. Scinocca, J. F. 1995 The mixing of mass and momentum by Kelvin- Helmholtz billows. Journal of the atmos. sci. 52, 2509–2530. Scotti, R. S. & Corcos, G. M. 1972 An experiment on the stability of small disturbances in a stratified shear layer. J. Fluid Mech. 52, 499–528. Seim, H. E. & Gregg, M. C. 1994 Detailed observations of a naturally occuring shear instability. Journal of Geophysical Research 99, 10049–10073. Seim, H. E. & Gregg, M. C. 1995 Energetics of a naturally occuring shear instability. Journal of Geophysical Research 1009, 4943–4958. Shih, L. H., koseff, J. R., Ivey, G. N. & Ferziger, J. H. 2005 Param- eterization of turbulent uxes and scales using homogeneous sheared stably stratied turbulence simulations. J. Fluid Mech. 525, 193–214. Smyth, W., Moum, J. & Caldwell, D. 2001 The efficiency of mixing in turbulent patches: inferences from direct simulations and microstructure observations. J. Phys. Oceanogr. 31, 1969–1992. Smyth, W. D. 1999 Dissipation range geometry and scalar spectra in sheared, stratified turbulence. J. Fluid Mech. 401, 209–242. Smyth, W. D. 2003 Secondary Kelvin-Helmholtz instability in weakly strat- ified shear flow. J. Fluid Mech. 497, 67–98. Smyth, W. D., Carpenter, J. R. & Lawrence, G. A. 2007 Mixing in symmetric Holmboe waves. J. Phys. Oceanogr. 37, 1566–1583. Smyth, W. D. & Moum, J. N. 2000 Length scales of turbulence in stably stratied mixing layers. Phys. Fluids 12, 1327–1342. Smyth, W. D., Nash, J. D. & Moum, J. N. 2005 Differential diffusion in breaking Kelvin-Helmholtz billows. J. Phys. Oceanogr. 35, 1004–1022. 122 Bibliography Smyth, W. D. & Peltier, W. R. 1991 Instability and transition in finite- amplitude Kelvin-Helmholtz and Holmboe waves. J. Fluid Mech. 228, 387– 415. Smyth, W. D. & Winters, K. B. 2003 Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr. 33, 694–711. Staquet, C. 1995 Two-dimensional secondary instabilities in a strongly strat- ified shear layer. J. Fluid Mech. 296, 73–126. Staquet, C. 2000 Mixing in a stably stratified shear alyer: two- and three- dimensional numerical experiments. Fluid. Dyn. Res. 27, 367–404. Stuart, J. T. 1967 On finite-amplitude oscillations in laminar mixing layers. J. Fluid Mech. 29, 417–440. Su, L. K. & Clemens, N. T. 1999 Planar measurements of the full three- dimensional scalar dissipation rate in gas-phase turbulent flows. Experiments in Fluids 27, 507–521. Su, L. K. & Clemens, N. T. 2003 The structure of ne-scale scalar mixing in gas-phase planar turbulent jets. J. Fluid Mech. 488, 1–29. Taylor, G. I. 1931 Effect of variation in density on the stability of superposed streams of fluid. Proc. R. Soc. Lond. 132, 499–523. Tedford, E. W., Carpenter, J. R., Pawlowicz, R., Pieters, R. & Lawrence1, G. A. 2009 Observation and analysis of shear instability in the fraser river estuary. J. Geophys. Res. 114, C11006, doi:10.1029/2009JC005313. Tennekes, H. & Lumley, J. L. 1972 A first course in turbulence. The MIT Press. Thorpe, S. 1987 Transition phenomena and development of turbulence in stratified fluids. J. Geophys. Res. 92, 5231–5245. Thorpe, S. A. 1968 A method of producing a shear flow in a stratified fluid. J. Fluid Mech. 32, 693–704. Thorpe, S. A. 1971 Experiments on the stability of stratified shear flows: miscible fluids. J. Fluid Mech. 46, 299–319. 123 Bibliography Thorpe, S. A. 1973 Experiments on instability and turbulence in stratified shear flows. J. Fluid Mech. 61, 731–751. Thorpe, S. A. 1978 The near-surfaceocean mixing layer in stableheating conditions. J. Geophys. Res. 61, 2875–2885. Thorpe, S. A. 1985 Labaratory observations of secondary structures in Kelvin-Helmholtz billows and consequences for ocean mixing. Geophys. As- trophys. Fluid Dyn. 34, 175–199. Thorpe, S. A. 2005 The turbulent ocean. Cambridge University Press. Turner, J. 1968 The influence of molecular diffusivity on turbulent entrain- ment across a density interface. J. Fluid Mech. 33, 639–656. Winant, C. D. & Browand, F. K. 1974 Vortex pairing: The mechanism of turbulent mixing layer growth at moderate reynolds number. J. Fluid Mech. 63, 237–255. Winters, K. B., Lombard, P. N., Riley, J. J. & D’asaro, E. A. 1995 Available potential energy and mixing in density-stratified fluids. J. Fluid Mech. 289, 115–128. Winters, K. B., MacKinnon, J. A. & Mills, B. 2004 A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Ocean. Tech. 21, 69–94. Woods, J. D. 1968 Wave-induced shear instability in the summer thermo- cline. J. Fluid Mech. 32, 791–800. Wunsch, C. & Ferrari, R. 2004 Vertical mixing, energy, and the general circulation of the oceans. Annu. Rev. Fluid Mech. 36, 281–314. 124
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Kelvin-Helmholtz instabilities in sheared density stratified...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Kelvin-Helmholtz instabilities in sheared density stratified flows Rahmani, Mona 2011
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Kelvin-Helmholtz instabilities in sheared density stratified flows |
Creator |
Rahmani, Mona |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | Kelvin-Helmholtz instabilities are the most commonly studied type of instability in sheared density stratified flows. Turbulence caused by these instabilities is an important mechanism for mixing in geophysical flows. The primary objectives of this study are the evolution of these instabilities and quantifying the mixing they generate using direct numerical simulations. The results are presented in three chapters. First, the evolution of primary Kelvin-Helmhlotz instabilities in two dimensions is studied for a wide range of Reynolds and Prandtl numbers, representing real oceanic and atmospheric flows. The results suggest that some properties of KH billows are predictable by a semi-analytical model. It is shown that a new Corcos-Sherman scale is a useful guide when simulating turbulent KH flow fields. The details of the mixing process generated by the evolution of Kelvin-Helmholtz instabilities as it goes through different stages, is analyzed. As the Reynolds number increases a transition in the overall amount of mixing is found, which is in agreement with previous experimental studies. This transition is explained quantitatively by the entrainment and mixing caused by three-dimensional motions, in addition to those resulted from the two-dimensional growth of the instability. The effect of Prandtl number on mixing is studied to understand the characteristics of high Prandtl number mixing events in the ocean; these cases have usually been approximated by low Prandtl number simulations. The increase in the Pradtl number has some significant implications for the evolution of the billow, the time variation of mixing properties, and the overall mixing. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-10-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0050710 |
URI | http://hdl.handle.net/2429/38098 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_rahmani_mona.pdf [ 3.91MB ]
- Metadata
- JSON: 24-1.0050710.json
- JSON-LD: 24-1.0050710-ld.json
- RDF/XML (Pretty): 24-1.0050710-rdf.xml
- RDF/JSON: 24-1.0050710-rdf.json
- Turtle: 24-1.0050710-turtle.txt
- N-Triples: 24-1.0050710-rdf-ntriples.txt
- Original Record: 24-1.0050710-source.json
- Full Text
- 24-1.0050710-fulltext.txt
- Citation
- 24-1.0050710.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0050710/manifest