Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Numerical simulations of linearly stratified flow around a sphere
(USC Thesis Other)
Numerical simulations of linearly stratified flow around a sphere
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
NUMERICAL SIMULATIONS OF LINEARLY STRATIFIED FLOW AROUND A SPHERE by Trevor Orr A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (AEROSPACE ENGINEERING) December 2014 Copyright 2014 Trevor Orr For Taylor Lily Orr & Duncan Thomas Orr ii Acknowledgments As allofthis comes toaclose, I think back onwhoI was atthe beginning ofallthis and what I would have written at that time. The words would have been similar, but there has been a fundamental shift in focus from the greatness of my own task onto the greatness of all the people who have helped me accomplish such a small task. Looking back, I can identify the important moments and times that have a direct impact on whatever path I walk today, and it is your actions that have made such adifference. Teachers, family, and, especially, friends, Iremember you all, and thank you. My parents get their own section. Without them, I would certainly not be where I am today. My dad’s love of flying, and my countless hours in the back of his Long-EZ and weekends atfly-ins arewhy I made the correct choice ofaerospace over computer science. My mom’s continuous focus on my scholastic aptitude, and countless arguments over why B’s were not acceptable, has resulted in this degree, thehighestlevelofacademicachievement. Ithankthembothfortheoverallpositive and nurturing environment that they were able to work toward and provide me. Professor Domaradzki has been instrumental in my completion of this degree. I have only begun to understand the level of responsibility involved in accepting, motivating,cultivating, andgraduatingastudent. Asabaseneed, Iamgratefulfor the financial support he secured as I pursued my studies. Professor Domaradzki is alsoresponsibleformyfirstformalintroductiontoCFDandguidingmysubsequent iii interest, understanding, and growth in this field. His depth of knowledge and understanding across multiple disciplines is so very much appreciated. My future self will undoubtedly continue thanking him for his interest and investment in me. It has been a pleasure and a privilege to work for and study under such a great man. Most importantly, I would like to thank my wife, Julie, for her endless support in my pursuit of this degree. This process has almost spanned our entire marriage. Her caring for our two children, maintaining our household, and looking after me has allowed the luxury of focus and late evenings that I required to succeed here. She is fantastic and a marvel. I could not have done this without you, Squish. Thank you. iv Table of Contents Dedication ii Acknowledgments iii List Of Tables vii List Of Figures viii Chapter 1: Introduction 1 1.1 The Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Non-Dimensional Equations of Motion and Their Similarity Parameters . . . . . . . . . . . . . . . . . . . 7 Chapter 2: Numerical Simulation of the Wake of a Towed Sphere 12 2.1 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.1 Solution of the Non-Dimensional Equations of Motion . . . . 13 2.1.2 Detached Eddy Simulation . . . . . . . . . . . . . . . . . . . 18 2.2 Computational Domain . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 24 2.2.2 Additional Domain Treatments . . . . . . . . . . . . . . . . 27 Chapter 3: Results 31 3.1 Wake of the Sphere at Re=200. . . . . . . . . . . . . . . . . . . . . 31 3.2 Wake of the Sphere at Re=1000 . . . . . . . . . . . . . . . . . . . . 37 3.2.1 Downstream Distance from the Sphere in x and Nt . . . . . 40 3.2.2 Averaging In The Near Wake . . . . . . . . . . . . . . . . . 41 3.2.3 Numerical Resolution at Re=1000, Fr=4 . . . . . . . . . . . 47 3.2.4 Length Scales of the Near Wake . . . . . . . . . . . . . . . . 51 3.2.5 The Near Wake Centerline . . . . . . . . . . . . . . . . . . . 57 3.2.6 The Near Wake in the Vertical and Horizontal Planes . . . . 72 3.2.7 Parameterizations vs. Simulation Data at Re=1000 . . . . . 77 3.3 Wake of the Sphere at Re=5000 . . . . . . . . . . . . . . . . . . . . 82 v Chapter 4: Conclusions 88 Bibliography 92 Appendix A Procedure to Initialize Sphere-less Simulations . . . . . . . . . . . . . . . 100 vi List Of Tables 2.1 Choice of simulation time step for Re/Fr parameter space . . . . . 17 3.1 Re=1000: Near wake vertical parameterization coefficients . . . . . 74 3.2 Re=1000: Near wake horizontal parameterization coefficients . . . . 74 vii List Of Figures 2.1 Sketch of physical problem relative to the sphere. . . . . . . . . . . 21 2.2 Sketch of the body-fitted coordinate system. . . . . . . . . . . . . . 21 2.3 Example of O-O grid topology.. . . . . . . . . . . . . . . . . . . . . 22 2.4 Low and high resolution grid comparison . . . . . . . . . . . . . . . 24 2.5 Sketch of the computational grid . . . . . . . . . . . . . . . . . . . 25 2.6 Example of viscous damping layer region . . . . . . . . . . . . . . . 30 3.1 Surface skin friction coefficient . . . . . . . . . . . . . . . . . . . . . 32 3.2 Re=200: Effects of stratification in sphere drag and wake internal wavelengths . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Re=200, Fr={0.4, 0.5}: Vertical vorticity in the horizontal plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4 Re=200, Fr=0.2: Drag coefficients and power spectra . . . . . . . . 35 3.5 Fr={1, 2, 4}: Instantaneous vertical velocity at Re=200 and Re=1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.6 Re=1000, Fr=∞: Statistics in the uniform density flow . . . . . . . 38 3.7 Sketches of stratified flow over the sphere . . . . . . . . . . . . . . . 42 3.8 Re=1000, Fr=4: Simulation statistics vs. DPIV methodology . . . . 43 3.9 Re=1000, Fr=4: Vertical velocity perturbation in the lee waves and the variance of the spatially averaged wake . . . . . . . . . . . 43 3.10 Re=1000, Fr=4: Effects of resolution on isosurface of λ 2 -criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 viii 3.11 Re=1000, Fr=4: Effects of resolution in perturbation quantities on wake centerline and edge . . . . . . . . . . . . . . . . . . . . . . 49 3.12 Re=1000, Fr=4: Effects of resolution in perturbation quantities on mean centerline defect velocity . . . . . . . . . . . . . . . . . . . 50 3.13 Re=1000, Fr=∞: Effects of resolution in computed half-widths of the near wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.14 Re=1000, Fr={1, 2, 4}: Computed half-widths of the near wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.15 Re=1000, Fr={1, 2, 4, 10}: The average body force located at the wake edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.16 Re=1000, various Fr: Computed half-widths of the near wake with parametric fitting . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.17 Re=1000, various Fr: Behavior and scaling of the centerline mean stream-wise velocity . . . . . . . . . . . . . . . . . . . . . . . 60 3.18 Re=1000, Fr=4: Anisotropy of centerline perturbation velocities in the near wake . . . . . . . . . . . . . . . . . . . . . . . 62 3.19 Re=1000, Fr={4, ∞}: Centerline stream-wise perturbation velocity components . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.20 Re=1000, various Fr: Centerline horizontal perturbation velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.21 Re=1000, various Fr: Scaling of vertical perturbation velocity by wake half-width . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.22 Re=1000, various Fr: Vertical perturbation velocity on the wake centerline, scaling, and fitted parameterizations . . . . . . . . 67 3.23 Re=1000, various Fr: Ratio of centerline density perturbations to vertical velocity perturbations . . . . . . . . . . . . . . . . . . . 68 3.24 Re=1000, various Fr: Centerline density perturbations . . . . . . . 71 3.25 Re=1000, Fr={4, 5, 16}: Parameterization of centerline density perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.26 Re=1000, various Fr: Numerical fit of planar time-averaged near wake data at various downstream distances. . . . . . . . . . . 78 3.27 Re=1000, Fr=∞: Parameteric fits versus simulation data . . . . . . 79 ix 3.28 Re=1000, Fr=6: Parameteric fits versus simulation data . . . . . . 81 3.29 Re=5000, Fr=4: Effect of azimuthal and spatial average at specific downstream location . . . . . . . . . . . . . . . . . . . . . . 84 3.30 Re=5000, Fr=4: Mean stream-wise velocity on wake centerline and horizontal and vertical slices at specific downstream location . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.31 Re=5000, Fr=4: Horizontal and vertical slices of perturbation velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.32 Re=5000, Fr=4: Centerline perturbation velocities . . . . . . . . . . 87 x Chapter 1 Introduction A large volume of work has been produced on the wake of submerged bodies trav- elling horizontally through stably stratified fluids of linear density gradients. The fluid can be categorized by its density profile under static conditions. A linearly stratified, stationary fluid has a density profile proportional to its depth, e.g.: ρ(z) = αz+ρ o (1.1) where α is the slope of the linear gradient profile, z is positive in the Cartesian ver- tical coordinate direction, and ρ o is a reference density. An isothermal atmosphere with an exponential density profile would fall into this category if the scale of the heightofvertical motionis smallcomparedtoitsscale height, see Turner 1973 [60]. A stratified fluid is classified as stable, neutral, or unstable. An unstable stratifica- tion produces a net buoyancy force on a displaced particle that encourages contin- ued advection away from its previous equilibrium position. A neutral stratification generates no net restorative force on the displaced particle. Uniform density fluids are categorized as neutrally stable. A stable stratification creates a net force that tends to return a displaced particle towards its original equilibrium position. If the 1 density gradient, dρ/dz, is negative with an increase in z, the fluid is stable. The gradient is linear and stable for all fluids under investigation. One of the distinctive features of the stratified flow, as opposed to the flow of a uniform density fluid, is that the stratification produces a gravity-induced restorative buoyancy force due to the displacement of particles away from their original hydrostatically balanced equilibrium positions. When the wake is oriented horizontally, this buoyancy force breaks the full axialsymmetry ofthe wake around the submerged body and may directly affect physical flow features ahead, around, andbehindthebody. Althoughmanyexperimentalconfigurationsforthetravelling body can be employed, the sphere is identified by Meunier & Spedding 2004 [37] and Meunier & Spedding 2006 [38] as a canonical configuration in the study of far wake characteristics. The wake of a sphere traveling horizontally through a linearly stratified fluid is categorized into three major regions: the near wake, the non-equilibrium (NEQ) regime by Spedding 1997 [52], and the far wake, also noted as the quasi-2D (Q2D) regime[52]. Eachoftheseregionshasbeenstudiedexperimentally, analytically,and numerically,butuntilnow,therehavebeennodetaileddescriptionsofself-contained multiple field components in the near wake of the sphere for any fluid of signifi- cant stratification. The parameter space in these studies is typically comprised of two variables. The motion of the sphere in the stratified wake is parameterized by the Reynolds number, Re, which is the momentum and viscous similarity param- eter of the sphere. The severity of stratification is parameterized by the internal Froude number, Fr, which characterizes the undisturbed density gradient of the fluid stratification. The current work is motivated by ongoing computational efforts to model the flowfield development ofthewake ofalinearly stratifiedfluid flowaroundasphere. 2 Thefirstnumericalstudytoincludeasphereexplicitlywithinthecomputationaldo- mainforalinearlystratifiedflowisthatofHanazaki1988[24]. Hanazakiperformed steady, Re = 200 simulations across a wide range of Fr = [0,200] to study the ef- fects of lee waves and stratification on the sphere itself. Recently, two quantitative studies have explicitly accounted forthe presence of the sphere inside the computa- tionaldomainforhighReynolds number simulations. Rottmanet al. 2010[43]per- formed simulations for Fr ={1,4} by use of an immersed boundary method. Sec- ond, Pasquetti 2011 [41] performed ”LES-like” simulations at Re = 10 4 , Fr = 50, accounting for the sphere with a ”psuedo-penalization” method. Other recent interest with respect to wake simulation in stratified flows is in the quantitative analysis of wake generated internal waves and the formation of far wake”pancake”eddies. Fung&Chang1996[19]performedturbulencemodellingof thefarwake with anisotropic closure schemes andstudied the formationofpancake eddies and the free surface reaction to their presence. Gourlay et al. 2001[21] con- tinued thenumerical investigations ontheformationofthepancake eddies byusing direct numerical simulation (DNS). Their Re =10 4 , Fr = 10 simulations showed that the far wake pancake eddies eventually form without the prior presence of coherent structures. The success of these simulations has encouraged successive simulations of the NEQ and far wake regions. Additional simulations investigating the NEQ and far wake using DNS include Brucker & Sarkar 2010 [9], the large eddy simulation (LES) of Dommermuth et al. 2002 [18], or multi-domain spectral element penalty methods from Diamessis et al. 2005 [16] and Diamessis et al. 2011 [17] who also initialize their stratified wake simulations without explicitly account- ing for the sphere. Eachnumericalstudythatinitializesitssimulationwithoutexplicitlyaccounting for the sphere utilizes an additional relaxation procedure [9, 16, 17, 18, 21] to allow 3 the turbulent velocity field to develop prior to analyzing the successive integration of the governing equations. Equipartition of turbulent energy is assumed during the initialization process, and an initial density field is not specifically prescribed. Pasquetti[41]initializedafarwakesimulationbyinterpolatingthewakeofasphere explicitly within the domain and found that the time/distance behind the sphere at which the pancake eddies form could change by a factor of five [41]. There is a continued interest in development of numerical stratified wake simu- lations, and related studies show a trend of considering the computational domain to represent a region that is physically close to the sphere; most commonly asso- ciated with the near wake region. The exact initial conditions and initialization procedures have varied: ad-hoc representation of the sphere [19], experimental fits of uniform density wakes (de Stadler et al. 2010 [15])[9, 18, 21], and extrapolation of post-NEQ experimental data[16, 17]. Experimentalstudiesofthewakedensityhavetypicallyusedconductivityprobes and/or rakes, but it has historically been a challenge, or not a focus, to captur- ing the full density field within the near wake. A full-field density description of 13 tablulated runs at 99 entries each is performed by Stockhausen 1966 [56] using a conductivity probe behind a self-propelled disk mounted to an ellipsoid body. Sarpkaya & Massidda 1997 [45] used conductivity probes mounted along the ex- perimental tank wall in a lightly stratified fluid at Fr∼O(100)to detect ensemble averages of the density fluctuations seen after the passage of a submarine-shaped self propelled body. Bonneton et al. 1996 [4] used a conductivity probe located at various points behind a towed sphere and showed that the vertical and per- turbation components in the near wake were highly correlated. The stratified far wake was studied by Bonnier et al. 2000 [6] by using a translating conductivity probe to study the internal structure of quasi-2D pancake vortices. Only recently, 4 Brandt & Schemm 2011 [7] have used an array of conductivity probes to measure the density perturbation field behind the sphere, albeit at Re& 10,000. The com- monality in these studies is that all tracked the density perturbation field in some respect but do not collect or reproduce the full velocity components of the flow field. Computational investigations into the near wake that explicitly account for the sphere body are in an advantaged position to generate the entire near wake velocity and density field across a wide, uncoupled (which Spedding 1996[55] high- lights as a useful feature) parameter space in an efficient manner. This work is a successor to thatofHanazaki [24]because the sphere is explicitly accounted for in the simulation domain, and the wake features are unsteady. The simulations at Re = 1000 do not rely on turbulence modelling or immersed bound- ary treatments, and a self-consistent data field of velocity and density is collected. The investigation is then concluded with a Re = 5000 turbulent near wake case that uses a detached eddy simulation (DES) formulation. 1.1 The Equations of Motion The motion of the stratified fluid in a gravity field is governed by the equation of continuity and the Navier-Stokes equations under the Boussinesq approximation: ∂u ∗ i ∂x ∗ i =0 (1.2) ∂u ∗ i ∂t ∗ + ∂u ∗ i u ∗ j ∂x ∗ j = − 1 ρ ∗ o ∂p ∗ ∂x ∗ i +ν ∗ o ∂ ∗2 u ∗ i ∂x ∗ j 2 −ρ ′ ∗ g ∗ δ i3 (1.3) andanequationofstatedescribingtheconvectionofthedensityperturbationvalue: ∂ρ ′ ∗ ∂t ∗ + ∂ρ ′ ∗ u ∗ i ∂x ∗ i = −w ∗ dρ ∗ (z ∗ ) dz ∗ +κ ∗ o ∂ 2 ρ ′ ∗ ∂x ∗ i 2 (1.4) 5 where time is t ∗ , the position vector is x ∗ i ~ e i with ~ e i as the standard Cartesian basis set, the velocity vector field is u ∗ i ~ e i , p ∗ is the perturbation pressure field, ρ ′ ∗ are the density perturbations, ρ ∗ o is the fluid reference density, ν ∗ o is the reference kinematic viscosity of the fluid, κ ∗ o is the reference value describing the thermal diffusivity of the fluid, and −ρ ′ ∗ g ∗ δ i3 accounts for the gravity force acting on the fluid in the vertical direction. Note that~ e i is implicit in Eqs. (1.3) & (1.4) and will remain implicit in subsequent discussions. Because of the familiarity of Cartesian coordinates,resultsarecommonlyreferencedintheCartesianset,andforreference, the following is explicitly stated: x ∗ i ~ e i ≡~ x ∗ =x ∗ ~ e 1 +y ∗ ~ e 2 +z ∗ ~ e 3 =x ∗ ~ i+y ∗ ~ j +z ∗ ~ k = (x ∗ ,y ∗ ,z ∗ ) (1.5) u ∗ i ~ e i ≡~ u ∗ =u ∗ ~ e 1 +v ∗ ~ e 2 +w ∗ ~ e 3 = u ∗ ~ i+v ∗ ~ j +w ∗ ~ k =(u ∗ ,v ∗ ,w ∗ ) (1.6) where ~ i is directed opposite to the steady motion of the sphere and (x ∗ ,y ∗ ,z ∗ ) =0 is located at the sphere center. By convention, all quantities denoted by an asterik are dimensional, and this convention is carried throughout the remainder of the work. To obtain Eq. (1.3), the fluid density is decomposed and hydrostatic balance is used to equate a background pressure gradient with a background stratification: ρ ∗ (x ∗ i ,t ∗ ) =ρ ∗ 0 +ρ ∗ (z ∗ )+ρ ′ ∗ (x ∗ i ,t ∗ ) (1.7) P ∗ (x ∗ i ,t ∗ ) =P ∗ (z ∗ ,t ∗ )+p ∗ (x ∗ i ,t ∗ ) (1.8) ∂P ∗ (z ∗ ,t ∗ ) ∂z ∗ =−g ∗ ρ ∗ 0 − g ∗ ρ ∗ (z ∗ ) (1.9) 6 where P ∗ (x ∗ i ,t ∗ ) is the total pressure and P ∗ (z ∗ ,t ∗ ) is the background pressure. The total density field is ρ ∗ (x ∗ i ,t ∗ ), ρ ′ ∗ (z ∗ ) is the undisturbed vertical density dis- tribution, and g ∗ is the sole component of gravity in the − ~ k direction. The left hand side of Eq. (1.3) is referred to as the inertial term because it represents the Eulerian advection of fluid particles, and the right hand side terms are referred to as the body force terms as they represent the mechanical forces influencing the fluid motion. 1.2 Non-Dimensional Equations of Motion and Their Similarity Parameters Equations (1.2), (1.3), and (1.4) are non-dimensionalized in preparation for nu- merical solution. The non-dimensional equations of motion for a stable, linearly stratified fluid under the Boussinesq approximation are: ∂u i ∂x i =0 (1.10) ∂u i ∂t + ∂u i u j ∂x j = − ∂p ∂x i + 1 Re ∂ 2 u i ∂x 2 j − ρ ′ Fr 2 D = 0 (1.11) and the equation of state describing the convection of the density perturbations becomes: ∂ρ ′ ∂t + ∂ρ ′ u i ∂x i = w+ 1 PrRe ∂ 2 ρ ′ ∂x 2 i (1.12) 7 Non-dimensionalization is achieved by creating the following non-dimensional variables: t = t ∗ U ∗ s D ∗ , x i = x ∗ i D ∗ , u i = u ∗ i U ∗ s , p= p ∗ ρ ∗ o U ∗ s 2 , ρ ′ = ρ ′ ∗ −D ∗ ∂ρ ∗ (z ∗ ) ∂z ∗ (1.13) Non-dimensional time is t, and U ∗ s and D ∗ are the sphere speed and diameter, re- spectively. The perturbation pressure is normalized by a representative dynamic head. The density perturbation is normalized by ∂ρ ∗ (z ∗ )/∂z ∗ , which is the up- stream, undisturbed fluid density gradient and is constant in a linearly stratified fluid. Field quantities (i.e., ~ u, p, and ρ ′ ) vary with space and time. The Brunt-V¨ ais¨ al¨ a frequency for an incompressible, stratified fluid is denoted by N ∗ (z ∗ ) and is defined by: N ∗ (z ∗ ) = −g ∗ ρ ∗ o ∂ρ ∗ (z ∗ ) ∂z ∗ 1/2 (1.14) In a linearly stratified fluid, N ∗ (z ∗ ) is constant with respect to z ∗ and is denoted simply as N ∗ . The Brunt-V¨ ais¨ al¨ a frequency is the maximum buoyancy frequency of the fluid and is in units of rads −1 . The non-dimensional parameters of Eqs. (1.11) & (1.12) are: Re = U ∗ s D ∗ ν ∗ o , Pr = ν ∗ o κ ∗ o , Fr = U ∗ s N ∗ R ∗ = 2 U ∗ s N ∗ D ∗ =2Fr D (1.15) The three non-dimensional parameters are Re, Pr, and Fr which are called the Reynolds number, the Prandtl number, and the Froude number, respectively. The Reynolds number is commonly interpreted as the ratio of inertial forces to viscous forces. Here, it is based on the sphere diameter and speed, and it is an indicator of the boundary layer behavior on the sphere. For a sphere moving 8 in a steady motion through a neutrally stratified fluid, a Reynolds number below 100,000 contains a sub-critical boundary layer. In the sub-critical boundary layer, the boundary layer flow is laminar from the forward stagnation point until the point of separation. When the Reynolds number increases past 400,000 then the flow is considered critical. In a critical boundary layer, the flow transitions into a turbulent boundary layer between the stagnation point and the separation point. The Reynolds number is also an indicator of the size of the Kolmogorov length scale. The Kolmogorov length scale represents the smallest overturning length scale in a turbulent flow and is proportional to Re −3/4 . Ideally, the simulation of a turbulent wake should contain sufficient numerical resolution to resolve both the largest and smallest scale structures within the flow. The Re defined within Eq. (1.11) is based on the sphere diameter and speed but also remains a primary indicator of the length scales present in the sphere wake. The Prandtl number, defined by Pr = ν ∗ o /κ ∗ o , is the measure of momentum diffusivity to thermal diffusivity. If the Pr > 1, Eq. (1.12) requires increased nu- merical resolution over the required resolution for Eq. (1.11) to resolve all dynamic scales in the flow. If Eq. (1.12) is parameterized with a Pr < 1, the range of dynamic scales of the density perturbations in Eq. (1.12) is expected to be smaller than those based on the Re in Eq. (1.11) alone. Salt-stratified water has an equivalent Pr ≃ 700 (although it is formerly the Schmidtnumber,Sc),thermallystratifiedwaterhasPr≃ 7,andairhasaPr≃ 0.7. Prior Prandtl number sensitivity studies for stratified non-wake configurations are reviewed by deStadler & Sarkar 2010 [15]. Their direct numerical simulation of a representative stratified wake at a parameter of Re = 10,000 investigated the sensitivity of the wake to Pr. They found that the mean flow and total wake energycharacteristicsareonlyweaklyinfluencedbythePratdownstreamdistances 9 beyond the near wake, and that increasing the Prandtl number to 7 has less overall effects than decreasing to Pr = 0.2 [15]. The behavior of the near wake for Nt≤ 3 is presumed unaffected by choice of Pr such that even simulations at Pr = 1 statistically represent the wake characteristics of a salt-stratified medium. Many experimental configurations of stratified flow behind a sphere are salt- based. It is shown in Section 3.2.2 that even at Re = 1000, the statistical descrip- tions of the Pr = 1 simulations compare well to the experiments performed in a salt-stratified tank at Re= 5000. Although this indicates success in some respects, Pr = 1 is not Pr = 700 and the behavior of the energy cascade for the density perturbations at the smallest scales and very long wake times will not necessarily correlate with experiments. This work is focused on the near wakes of submerged bodies in stably stratified fluids, regardless of the stratification agent, as this is the origin of the motivation for the current work. Performing simulations of the far wake is not a priority here and the choice of Pr = 1 is justified. Throughout this work, the internal Froude number, Fr, is based on the sphere radius, R ∗ , and Fr = U ∗ s /(N ∗ R ∗ )= 2Fr D . The square of the Froude number may be viewed as a representation of the mean kinetic energy of a particle to the energy required for that particle to move over the top of the sphere when initially located withinthez ∗ = 0plane. Fralsorepresentsotherbuoyantcharacteristicsoftheflow. IfaparticleisfreetooscillateverticallyatafrequencyN ∗ ,thenitsbuoyancyperiod must be T b ∗ = 2π/N ∗ . Because N ∗ is constant, non-dimensionalization consistent with Eq. (1.11) results in T b =πFr, and likewise, ω b = Fr −1 . Thus, the Fr is also the characteristic angular frequency of the oscillation. Inalinearized system foralinearlystratifiedfluid, thenon-dimensionalinternal wave frequency of an oscillating particle is ω b . If the observer’s frame of reference is attached to the steadily, horizontally translating sphere, and an observed wave 10 is now stationary with respect to the observer, then the wave vector aligned with the wake centerline is k = ω b . This is due to the non-dimensional Doppler relation (see Lighthill 1979[33], Hopfinger et al. 1991 [25], or Bonneton et al. 1993 [3]): ω obs. =ω b −k i (1.16) and to the internal wave dispersion relation in a linearly stratified fluid [33, 25, 3]: ω b = Fr −1 cos(θ) (1.17) where θ is the defined as the angle ofthe wave-vector with respect to the horizontal planedefinedbyz =0. Inthissteady,linearapproximation,anobservedstationary wave would be related to buoyant wavelength, λ, by ω b = k = 2π/λ, or λ = 2πFr. Thus, Fr is also the metric of the characteristic buoyant wavelength in the wake of the sphere. 11 Chapter 2 Numerical Simulation of the Wake of a Towed Sphere 2.1 Numerical Method The numerical method solves the discretized form of Eqs. (1.11) - (1.12) using a finite-difference, alternating direction implicit (ADI), pseudo-transient generalized curvilinear partial transformation formulation on a body-fitted, non-staggered grid (e.g., Constantinescu 1998 [12]) in the parallelized framework seen in Koken & Constantinescu 2009 [30]. Inafullcoordinatetransformationformulation, allcomponents aretransformed into the body fitted basis set. In a partial coordinate transformation, only the spatial derivatives of the governing equations are transformed, i.e.: ∂ ∂x i ~ e i = ∂ξ j ∂x i ∂ ∂ξ j ~ e i (2.1) The spatial derivatives are taken along the body-fitted coordinate system which is aligned with the body-fitted grid, but the governing equations remain in the Cartesian basis set. 12 2.1.1 SolutionoftheNon-DimensionalEquationsofMotion For the linearly stratified fluid solved in an orthogonal body-fitted coordinate sys- tem, Equations (1.10)-(1.12) are partially transformed as: J ∂ ∂ξ i V i J = 0 (2.2) ∂u i ∂t +V j ∂u i ∂ξ j + ∂ξ j ∂x i ∂p ∂ξ j − 1 Re J ∂ ∂ξ j 1 J g jk ∂u i ∂ξ k + ρ ′ Fr 2 D δ i3 = 0 (2.3) ∂ρ ′ ∂t +V j ∂ρ ′ ∂ξ j −u 3 − 1 PrRe J ∂ ∂ξ j 1 J g jk ∂ρ ′ ∂ξ k =0 (2.4) where J is the Jacobian of transformation, the contravariant metric tensor g jk is defined: g jk = ∂ξ j ∂x i ∂ξ k ∂x i (2.5) and the contravariant components of velocity are defined: V j = u i ∂ξ j ∂x i (2.6) The Prandtl number is a fluid property and represents the ratio of momentum to thermal diffusivity of the fluid. If Pr > 1 then the range of all dynamic scales in Eq. (1.12) will increase over what is indicated by the presence of the Re in the Navier-Stokes equations, Eq. (2.3). As discussed in Section 1.2, setting Pr = 1 in the current simulations is expected to physically represent the large scale flows in the near wake region but are not expected to resolve the finest scales of higher Pr fluids. 13 The numerical implementation of Arnone et al. [1] is followed. Equations (2.3) & (2.4) are numerically integrated by adding a pseudo-transient term: ∂Q ∂τ +GE(Q) = 0 (2.7) where Q represents any primitive variable, GE() is its governing equation, and τ is the additional psuedo-time variable. A formulation containing both the psudeo- time and temporal derivative is the so-called dual-time stepping scheme. When Eq. (2.7) reaches a steady state in pseudo-time, the governing equations are iden- tically satisfied regardless of whether GE() contains a temporal derivative, ∂Q/∂t, or not. The following notation for specific operators within the governing equations are introduced: UN( )= ∂ ∂t (2.8) NL( )= V j ∂ ∂ξ j (2.9) PR() = ∂ξ j ∂x i ∂ ∂ξ j (2.10) VS( )= J ∂ ∂ξ j 1 J g jk ∂ ∂ξ k (2.11) Both continuous and discrete forms the operators will use this notation. The momentum equation, Eq. (2.3), is written in discretized form: u i m+1 −u i ∼ Δτ − u i ∼ −u i m Δτ +UN (u m i )+NL m (u i m+1 ) = −PR(p m+1 ) − ρ ′m Fr 2 D δ i3 + 1 Re VS(u i m+1 ) (2.12) 14 and the density perturbation equation, Eq. (2.4), is written: ρ ′ m+1 −ρ ′ m Δτ +UN(ρ ′ m+1 ) +NL m+1 (ρ ′ m+1 ) = u m+1 3 + 1 Re VS(ρ ′ m+1 ) (2.13) where ”m” indicates the explicit evaluation of a quantity at the m-th discrete level in the pseudo-time integration and ”∼” indicates an intermediate quantity state. The temporal derivative is discretized using a second order backwards time difference: UN(q m )= 3q m −4q n +q n−1 2Δt (2.14) where ”n”indicatestheprevious timelevel inthediscretized integration, and”n-1” the previous level to that. The convective velocity in the non-linear momentum equation convective term is phase-lagged in pseudo-time for the momentum equa- tions which is denoted by NL m (). The non-linear operator is spatially discretized using either a 2 nd -orderor 5 th -order upwind scheme. Both PR() and VS()contain spatial derivatives on stencils comprised of central differences. The discretizations of the operators NL( ), PS(), and VS( ) appear in Constantinescu 1998 [12]. Flow quantities are linearized for implicit treatment in the pseudo-time integra- tion: u ∼ i =u m i +Δu i (2.15) p m+1 =p m +Δp (2.16) ρ ′ m+1 =ρ ′ m +Δρ ′ (2.17) 15 The linearized form of Equation (2.12) is split to form a convection-diffusion step followed by a continuity step [12]: 1 Δτ 1+ΔτNL m ()+ Δτ Re VS() Δu i = NL m (u i m )−PR(p m )+ 1 Re VS(u i m ) − ρ ′m Fr 2 D δ i3 −UN (u m i ) (2.18) u i m+1 −u i ∼ Δτ =−PR(Δp) (2.19) FollowingthemethodologyofSortiropolous&Patel[48],anartificialdissipation term is added to the discrete continuity equation to avoid an unphysical pressure field because the pressure is collocated with the velocity at each node. The dis- cretized continuity equation, Eq. (2.2), is modified and written in implicit form as: DIV u i m+1 =ǫ(L− ¯ L)p m+1 (2.20) where DIV () is the divergence operator of Eq. (2.2). The pressure operators L and ¯ L are stencils of separate size that introduce a stabilizing dissipation term to the pressure equation. The specific discretizations of the operators L and ¯ L are found in Constantinescu 1998 [12]. By incorporating Eqs. (2.16) & (2.19) into Equation (2.20): λ ¯ L−L +DIV [−ΔτPR()] Δp =ǫ L− ¯ L p m +DIV [u ∼ i ] (2.21) which provides the updated pressure based on the intermediate velocity field. The coefficient of the pressure dissipation term, ǫ, may vary between 0 <ǫ< 1 [12, 48]. A dissipation term of ǫ = 0.05 is selected based on the previous direct numerical simulations of flow around a sphere at Re =300 by Johnson & Patel 1999 [28]. 16 Re = 200 Re = 1000 Fr 0.10 0.15 0.20 0.275 0.35 [0.40,∞] [1.0,∞] Δt 0.005 0.005 0.01 0.0125 0.0125 0.025 0.01 Table 2.1: Choice of Δt for each Re/Fr pairing. The density perturbation equation, Eq. (2.4), is linearized as: Δρ ′ Δτ +NL m+1 (Δρ ′ )+ 1 PrRe VS(Δρ ′ ) = NL m+1 (ρ ′ m )+ 1 PrRe VS(ρ ′ m ) +u m+1 3 −UN ρ ′ m (2.22) and also implicitly formulated: 1 Δτ 1+ΔτNL m+1 ()+ Δτ PrRe VS() Δρ ′ = NL m+1 (ρ ′ m )+ 1 PrRe V(ρ ′ m ) −u m+1 3 −UN ρ ′ m (2.23) In summary of the algorithm, once the momentum and pressure field are inte- grated to level ”m+1”, the density perturbation equation is then advanced to level ”m+1”. The pseudo-transient process is repeated for each unsteady time step until any flow field quantity is ΔQ ≃ 0. As lim ΔQ→0 Q m+1 = Q m ; indicating a steady stateinthepsuedo-transient termandasatisfaction ofthegoverning equationsuch that when Q m+1 =Q m , then Q m+1 =Q n+1 in Eq. (2.18) and (2.22). The pseudo-transient residuals are based on Δ ~ Q max where ~ Q represents the primitive variables in Eqs. (1.11) & (1.12)and arerequired tobe less thanorequal to 10 −4 to obtain the unsteady converged solution at each physical timestep. The required iterations to converge varies between each Fr case but are typically less than 100 iterations in the Re = 1000 cases. 17 2.1.2 Detached Eddy Simulation Reynolds Averaged Navier Stokes (RANS) formulations are used extensively to simulate turbulent flows. Because RANS is formulated to average flow field prop- erties under the Reynolds decomposition, they are suited for flows where only the mean flow properties are desired (e.g. Wilcox pg. 30 [62]). RANS formu- lations are useful in modelling attached near-wall boundary layers but are not suitable if the time-dependent, spatially developing structures within the wake are of interest, especially in the axisymmetric case of the sphere as is shown by Constantinescu et al. [11]. An alternative approach that attempts to represent the large-scale structures of the wake is large eddy simulation (LES). LES resolves the largest scale motions in space and time but models the interactions of scales of motion that are spatially unresolvedbythecomputationalgrid. Thespatiallyunresolvedmotionsarereferred to as the sub-grid scale motions. A major drawback of LES occurs in modelling flows in the presence of a wall and, consequently, the boundary layer of the sphere. This is because the dynamic scales within the boundary layer are small compared to the large-scale motions of the wake. The spatial resolution required to include near-wall LES is computationally prohibitive because of this difference in scale. Detached eddy simulation (DES) is a blend between RANS and LES that uses RANS as a wall-model in the boundary layer and an LES sub-grid scale model to resolvetheunsteadyfeaturesinthewake. DESisacategoryofturbulencemodelling with a variety of formulations (e.g. Spalart et al. 1997 [50] or Strelets 2001 [57]), and it is not specific to a particular turbulence model; only that it behaves in the prescribed RANS/LES manner. 18 Only the Re = 5000 investigation uses the Spalart-Allmaras DES (SA-DES) formulation ofConstantinescu 2004[13], previously applied to uniform density flow over a sphere at Re = 10,000, on a grid resolution of 101×101×201. Near the sphere boundary, SA-DES uses the one-equation Spalart-Allmaras RANS formu- lation of Spalart 1994 [49] to predict boundary layer and separation behavior in the boundary layer region. Outside of the boundary layer region, the model au- tomatically switches to an eddy resolving LES-type mode. The LES-type model behaveslikeaSmagorinskymodel(Smagorinsky1963[47]),whichinitselfattempts to represent the Kolmogorov spectrum for sub-grid scales as determined by Lilly in 1967 [34]. The SA-DES eddy viscosity is also used for the density perturbation eddy viscosity since Prandtl number is set formally to one. The addition of the DES turbulence model to the equations of motion, Eq. (2.3), modifies the viscous operator to be time dependent; i.e. VS( )→ VS m () in the following way: VS m ( )= 1 Re + 1 Re m T J ∂ ∂ξ j 1 J g jk ∂ ∂ξ k where Re m T = (U ∗ s D ∗ s )/ν ∗m t , where ν ∗m t is the eddy viscosity at iteration ”m”. Consequently, Eq. (2.3) stays effectively unmodified to the user, and the tur- bulence model is turned off as the term 1 Re T → 0. The SA-DES formulation for the eddy viscosity is reproduced from Nikitin et al. 2000 [40] as follows: ∂˜ ν ∂t +V j ∂˜ ν ∂ξ j =c b1 ˜ S˜ ν + 1 σ (▽·(ν + ˜ ν)▽ ˜ ν)+c b2 (▽˜ ν) 2 −c w1 f w ˜ ν ˜ d 2 (2.24) 19 where: ν t = ˜ νf ν1 , f ν1 = χ 3 χ 3 +c 3 ν1 , χ = ˜ ν ν ˜ S = S + ˜ ν/κ ˜ d 2 , f ν2 =1−χ/(1+χf ν1 ) f w =g 1+c 6 w3 g 6 +c 6 w3 1 6 , g =r+c w2 r 6 −r , ˜ r = ˜ ν ˜ Sκ 2˜ d ˜ d = min(d w ,C DES △) and S is the magnitude of the vorticity; d w is the physical distance to the wall; △ is the local grid spacing; c b1 = .1355, σ = 2/3, c b2 = 0.622, κ = .41, c w1 = c b1 /κ 2 +(1+c B 2)/σ, c w2 = 0.3, c w3 = 2, and c ν1 = 7.1aremodeling constants. The value of C DES is generally set at 0.65 following the work done by Shur et al. [46] for homogeneous turbulence. 2.2 Computational Domain The simulation grid is constructed in the naturally orthogonal spherical coordinate system in an O-O grid topology, where the analytic relations between r, θ, and φ (radial, polar, and azimuthal, resp.) and x, y, z are known. A sketch of the body-fitted coordinate system is found in Figure 2.2. The low resolution grid has a uniform distribution of points in the polar direc- tion, andtherelationshipbetween theCartesian coordinates,~ x, andthecurvilinear coordinates, ~ ξ, (e.g., [28]) are: r =R min +(R max −R min ) tanhk 1 ζ−1 ζmax−1 −1 tanhk 1 +1 (2.25) 20 Figure 2.1: Sketch of physical problem relative to the sphere. Figure 2.2: Sketch of the body-fitted coordinate system. 21 Y X Z Figure 2.3: Example of O-O grid topology. 22 θ = π ξ−1 ξ max −1 (2.26) φ =2π η−1 η max −2 (2.27) subsequently: x = rcosθ, y =r sinθcosφ, z =rsinθsinφ (2.28) where r =r ∗ /D ∗ , R max =R ∗ max /D ∗ , and θ and φ are in radians. Figure 2.3 gives a graphical representation of a typical physical grid used. These simulations are per- formed on a grid size of N ξ ×N η ×N ζ = 120×120×201. The computational grid is constructed from the known relations, i.e. Eqs. (2.25)-(2.28),and a visualization ofthecomputationalgridisinFigure2.5. Convective termsforthegoverningequa- tions are discretized using a 2 nd order upwind scheme. The equations of motion, Eqs. (1.11) & (1.12), are solved directly without any additional explicit turbulence modelling. R max is chosen to approximate a far-field condition in the simulations at R max = 60, and R min =0.5 is the surface of the sphere. The simulations at this resolution run across a fairly wide parameter space of Fr for Re ={200,1000}. Two 5th-order upwind, high fidelity simulations are performed at a higher res- olution of N ξ ×N η ×N ζ = 180×180×551 for Re = 1000 and Fr ={4,∞} using an hyperbolic stretching function similar to the radial one in Eq. (2.25) and for the polar grid point distribution as similar to Mittal 1999 [39]: θ =π tanhk 2 ξ−1 ξmax−1 −1 tanhk 2 +1 (2.29) The radialstretching parameterinthehigh resolutioncases is modifiedtok 1 =3.3, and the polar stretching parameter is k 2 = 1.5. R max is set to 20 for the radial stretching function in these cases for ζ ≤ 501, and additional 50 grid points are 23 Figure 2.4: Low and high resolution grid comparison for 6.25 . x . 7.5 and 0≤ y. 0.6. Dashed line is low resolution grid. added to the domain up through R max = 60 under a different radial stretching function. The grid size in the radial direction does not see much transition in grid size where these two stretching functions meet. This distribution concentrates two- thirds of the grid points behind the sphere, increases the uniformity of the grid point distribution along the wake centerline, and increases the grid resolution in the region of interest over the low resolution cases by about a factor of ten. A qualitative comparison for a radial and aziumthal slice is shown in Figure 2.4 for a downstream centerline near wake location of x≃ 7. 2.2.1 Boundary Conditions There is a periodic overlap of one grid point in the azimuthal direction such that η 1 = η max−1 and η 2 = η max . Because of this periodicity, there are no explicitly stated boundary conditions in this direction. An analytic singularity is located at 24 Figure 2.5: Sketch of the computational grid θ = {0 = ξ min ,π = ξ max } is created by the choice of grid construction coordinate system in the physical domain because of the Jacobians in Eqs. (2.2) - (2.4). On the surface of the sphere, ζ min ≡R min , a no-slip condition is enforced: u i =0 (2.30) If the inflow boundary is far upstream from the sphere, it is sufficient to use DirichletconditionsforlowReynoldsnumber, linearystratifiedcases, inaccordance with Hanzaki 1988 [24] who used an outer domain distance of 20 sphere diameters [24] at Re = 200. This distance was also found sufficient for low Reynolds number non-stratified simulations in [28] and also for turbulent, high Reynolds number 25 (Re > 10 4 ) simulations [13]. For the inflow (ζ max ≡ R max , ξ ≡ θ = [0,0.55π)), Dirichlet boundary conditions are enforced for the momentum equations: u i =δ i1 (2.31) The second condition on the domain boundary away from the sphere is referred to as the outflow boundary condition and will always extrapolate information from insidethecomputationaldomain. Theoutflowboundaryiscontainedontheoutflow surface defined on ζ max ≡R max , ξ ≡θ =[0.55π,π]: Du i Dt + ρ ′ Fr 2 D δ i3 = 0, (2.32) and represents theadvection ofa fluid particle andits tendency toaccelerate in the vertical due to any buoyancy forcing. All variables for momentum, the turbulence model,anddensityareextrapolatedfromtheinterioronthesingularitylinespresent at θ ={0,π} and averaged over the repeated nodes. The pressure boundary condition is controlled by using the normal component of the momentum equation on the computational boundaries located at R min and R max to enforce conservation of mass as required by Gresho & Sani 1987 [23]: ∂P ∂x i = − ∂u i ∂t −u j ∂u i ∂u j + 1 Re ∂ 2 u i ∂x 2 j − ρ ′ Fd 2 ~ e i · ~ ζ (2.33) This provides an explicit equation for the pressure variable when put into dis- cretized form along both the sphere surface and the far field domain boundary. The pressure is extrapolated from the interior along the singularity lines present at θ ={0,π}; subject to the conservation of mass constraints and similar to what is 26 recommended in Johnson & Patel 1999 [28]. Along the periodic overlap in the η direction, periodicity is applied to the pressure field as well. The density perturbations have a zero gradient normal to the surface of the sphere: ∂ρ ′ ∂~ n = ∂ρ ′ ∂ζ = 0 (2.34) wheretheorthogonalnatureofthegridisexploitedtoachieveasimpleexpressionin ζ. Physically, thisconditionrequires azeroheatfluxconditiononthesurfaceofthe sphere. Pasquetti 2011[41]allows theimmersed sphere tobe thermally conductive, but Hanzaki 1988 [24] considers the non-diffusive equivalent of a salt-stratified medium such that the initial conditions for the density distribution on the sphere surface are unchanged. Because the domain is body-fitted in a diffusive system, a control, zero-flux condition for the density perturbations is the appropriate choice. A zero value is enforced at the inflow for the density perturbations which rep- resents the undisturbed upstream condition. The density perturbation boundary conditions defined on the outflow are: Dρ ′ Dt −u 3 = 0 (2.35) TheformofthisconditionissimilartoEquation (2.32),buttheverticalcomponent ofvelocity represents the advection ofthe background density, i.e. Kundu &Cohen 2004 [31]. 2.2.2 Additional Domain Treatments The physical problem being modelled exists within a theoretically infinite domain and will contain lee and internal waves that propagate from around the sphere 27 and within the turbulent wake structure. An infinite computational domain would allow any internal waves topropagateout towards infinity. Introducing acomputa- tionaldomainandanartificialdomainboundarycomplicatessimulationscontaining waves because waves may reflect from the numerical boundary. A turbulent wake will produce internal waves adjacent to the computational boundary. The exact correction required to ensure the radiation condition to be completely satisfied in a fully three dimensional domain is an area of active investigation, and a review is available from Givoli 2004 [20]. A possible method would be to use the eddy viscosity of the turbulence model itself to damp turbulent oscillations in the wake, which will damp turbulent oscil- lations towards or near the boundary. The idea is to use the ”switch” inherent to the DES formulation Eq. (2.24) to explicitly turn off the LES mode towards the outflow boundary and return the turbulence modeling to RANS; explicitly: ˜ d = d w for~ x∈ Ω(~ x damp ) (2.36) where d w is the physical distance to the wall as in the original SA-RANS model [49] and Ω(~ x damp ) is the domain in which damping is desired. This would create a RANS simulation at the outflow and would maintain the DES behavior desired in the near wake region. An additional possibility of using the DES formulation as a damper is to increase the value of C DES near the computational outflow boundary. Constantinescu & Squires 2004 [14] show that by increasing the value of C DES the smallest scales ofenergy aredamped. WithintheDESframework, increasing C DES may create a damping characteristic similar to the RANS switch approach. InpreliminarysimulationsusingtheRANS-as-dampercondition, theoscillatory nature ofthe wake decreased, but the presence of lee waves is still detected towards 28 thedomainboundary. Forleewavelengths thatcouldcontainsignificantamplitude, such as those at Fr = 4, one characteristic lee wavelength is 12 sphere diameters, and this length increases with increasing Fr. Within the range 0 . Fr . 100, the characteristic lee wavelength is directly proportional to Fr, which means that reflections of the lee wave could occur for any practical domain size. Itisdesirabletofindamethodcapableofreducingallpotentiallyreflectedwaves while at the same time maintaining a certain simplicity in its implementation. A tempting methodology is to introduce some artificial viscous damping to the equations of motion, Eqs. (2.3) & (2.4), near the region of the boundary, but far enough from the region of interest so as not to adversely affect the investigation’s observations: NS(q) =−ν(~ x)q (2.37) NS(q) = µ (~ x) ∂ 2 q ∂x 2 (2.38) where NS(q) represents the original Navier-Stokes bracketed term in Eq. (2.3) operating on any desired quantity q that is to be damped. ν(x) and µ (x) are the dampingcoefficientswhichvaryspatially. IsraeliandOrszag1981[26]showthatthe damping term must vary smoothly in space to avoid reflections off of the damper interface itself, and they show that Eq. (2.37) is superior in its wave damping properties. Along with the gradually varying damping coefficient, a long length of damping layer in conjunction with a sufficient placement of nodes within the layer mitigates numerical reflection of waves from the computational boundary and their interaction within the computational region of interest, i.e. Bodony 2006 [2]. A damping layer is used in the current simulations due to the potential of numerical boundary reflections of lee and internal waves that are generated by the 29 sphere and its wake. The layer is implemented by adding a forcing term to the governing equations, Eqs. (2.3) & (2.4): GE( ~ Q)+ν(~ x)( ~ Q− ~ Q ave )= 0 (2.39) whereGE isthegoverningequationforanycomponentof ~ Q = (u,v,w,ρ ′ )andν(~ x) containsthespatiallyvaryingdampingcoefficient. Bychoice: ~ Q ave = (U(~ x),0,0,0), whereU(~ x)isthespatiallyvarying, time-averagedstream-wise velocity component. Thedampingcoefficientiszeroforx≤ 15,|y|≤ 8,|z|≤ 8andreachesitsmaximum just before the outer domain limits. Figure 2.6: Example of viscous damping layer region in the vertical plane (z ≥ 0); 0 indicates no damping, 1 indicates full damping. Mesh density visible in vertical plane (z ≤ 0) for qualitative comparison. 30 Chapter 3 Results 3.1 Wake of the Sphere at Re=200 ThewakeofaspheretravellingsteadilythroughauniformdensityfluidatRe = 200 is commonly studied, such as in Taneda 1956 [58], Johnson & Patel 1999 [28] and Lee 2000 [32]. The wake is considered steady and generates an axisymmetric recirculationregionattachedtotheleeofthesphere. Thelowresolutionsimulations indicate the non-dimensional length of this laminar recirculation region, measured from the center ofthe sphere, isL b =1.95. The skin friction coefficient can be used insteadyflowsasanindicatorforseparation. Instantaneousskinfrictioncoefficients for the steady homogeneous case is shown in Figure 3.1a. The separation occurs at approximately 120 degrees in both the horizontal and vertical, whose lines overlap in that Figure. The stream-wise drag coefficient is computed to be C Dx = 0.757. These quantities are consistent with prior findings [58, 28, 32]. To ensure the simulations are producing viable results for stratified flow cases, it is important to compare the code behavior against existing data sets. A natural startingpointistocomparesimulationresultswithHanazaki[24]atRe = 200. The flow field is steady, except in the case of very stratified flows which should exhibit 31 C f 0 50 100 150 0 0.05 0.1 0.15 0.2 0.25 0.3 Figure 3.1: Re=200: Surface skin friction coefficient, C f vs. angle, θ, on the (−) vertical(z ≥ 0) and (−·−) horizontal (y ≤ 0) plane. θ is taken from the negative x-axis. unsteady vortex shedding as noted by Chomaz et al. 1993 [10]. Results obtained for time-averaged ΔC Dx as a function of Fr are presented in Figure 3.2a where: ΔC Dx =C Dx (Re,Fr)−C Dx (Re,∞) (3.1) where C Dx (Re,∞) is the stream-wise drag coefficient of the sphere in a uniform density fluid at a particular Re. The overbar of ΔC Dx denotes a time-averaging for the unsteady cases. By Eq. (3.1), ΔC Dx represents the change in drag coefficient of the sphere due to the stratification ofthe fluid. The simulation agrees with prior investigation [24, 36,43] values of ΔC Dx for Fr& 0.5. When Fr. 0.5, the current simulations disagree completely with Hanazaki [24], corroborate those of Rottman et al. [43], and reproduce the experimental results of Lofquist & Purtell 1984 [36] until Fr < 0.2. When Fr < 0.5, the flow has a horizontal orientation, and vertically oriented, coherent vortices are shed, as seen in Figure 3.3a; whereas when Fr ≥ 0.5 in 32 (a) (b) 0.1 1 10 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Fr ΔC D x 1 2 0 2 4 6 8 Fr λ Figure 3.2: (a) Change in stream-wise drag coefficient due to stratification. (−) Hanazaki[24], (−−) Lofquist & Purtell [36], {× (low-res), + (med-res)} Rottman et al. [43], low resolution simulations at (•) Re=200, △ Re=1000. (b) Calculated λ: (−) Hanazaki [24], (−−) linear theory, low resolution simulations {+,△,•} at z ={0.5,1.0,1.5} for Re=200, at z = 1.0 for Re=1000. 33 Figure 3.3b there is no vortex shedding. At Fr = 0.2, the simulations produce a St y = 0.19 at Re = 200. Chomaz et al. [10] also noted a general Strouhal number of approximately 0.2 at Re = 2000 and Fr = 0.3. In this regime, there are small- amplitude oscillations ofthe stream-wise dragcoefficient where St x ≃ 2St y =0.37; as shown by Figure 3.4. Simulation results are in agreement with experimental values [36] of ΔC Dx until Fr <0.2. When Fr < 0.2, the current simulations likely have an insufficient resolution for this severely stratified fluid. The current simulations and the low-resolution simulations [43] at Fr = 0.1 concur with a value of ΔC Dx ≃ 0.80. Between the low- and mid-resolution cases [43] at Fr = 0.1, the discrepancies in ΔC Dx and experiment [36] tend to disappear; indicating that increased resolution in future simulations would allow further investigation into these very low Fr flows if they were of interest. Hanazaki 1988 [24], Greenslade 2000 [22], and Voisin 2007 [61] suggest that a Re effect accounts for the discrepancy between experimental data [36] and the computations of Hanazaki [24]. For Re = 200, current results remain consistent with experiments until low Fr values. Review of the experimental data [36] in- dicates that for Fr . 0.25, the highest Re case is Re ≃ 550, and the lowest Re case is Re ≃ 175 at Fr = 0.125. Turbulent fluctuations in the wake typically occur around Re = 1000, and a Re dependence does not sufficiently explain the discrepancies between the ΔC Dx curves. The current simulations show, via Figure 3.2b, comparable lee wavelength, λ, values to both Hanazaki[24] and linear theory. Lines of constant phase are lo- cated where the vertical component of velocity is zero according to linear theory as ω = N∗sin(θ) [33], and θ is the angle of the motion of the fluid particle with respect to the horizontal; making lines of w = 0 equivalent to lines of constant 34 (a) (b) Figure3.3: Re=200: Instantaneousvertical vorticity, ω z ,onplanedefinedbyz =0: (a) Fr = 0.4 (b) Fr =0.5. (a) (b) 0 10 20 30 40 50 0 0.5 1 1.5 t C D 0 0.2 0.4 0.6 10 −1 10 0 10 1 10 2 10 3 St Power Figure 3.4: Re=200, Fr=0.2: (−) C Dx , (−·−) C Dy , (−−) C Dz . (a) Components of drag coefficient. (b) Power spectra (C Dz omitted). 35 (a) (b) (c) (d) (e) (f) Figure 3.5: Contours of instantaneous vertical velocity, w, on plane at y = 0. The sphereiscenteredatx = y =z =0. Re=200: (a) Fr =1 (b) Fr =2 (c) Fr =4. Re=1000: (d) Fr =1 (e) Fr =2 (f) Fr = 4. 36 phase. Measurements of λ are calculated by the distance between inflections of w(~ x,t)= W(~ x)= 0 (i.e, Figure 3.5) at different vertical heights on the plane y =0 at z ={0.5,1.0,1.5}, where the overbar indicates time-averaging for the unsteady flow cases. Considering the differences between ΔC Dx evident in Figure 3.2a and the agreement on λ in Figure 3.2b, sufficient resolution of the lee waves is also not a cause for discrepancy with Hanazaki [24] when Fr < 0.5. By Figure 3.2b, there is no significant departure of calculated lee wave length from linear theory due to the vortex shedding. Brighton 1978 [8] postulated that there is a lack ofdirect interaction between lee waves and vortex shedding, andthis is partly confirmed as the Hanazaki low-resolution (32×32×62) [24] simulations areabletoreproducetheexpected λvaluesinFigure 3.2abutnottheexperimental ΔC Dx seen in Figure 3.2b. In the range of highly stratified, Fr ≃ 0.2, to uniform density fluids, Fr =∞, thecurrent simulations accurately reproduce thecharacteristics ofastratifiedfluid. Results for Re= 1000 are focused on Fr ≥ 1. Severely stratified cases, Fr < 0.2, require higher numerical resolution and are not of current interest. 3.2 Wake of the Sphere at Re=1000 At Re = 1000, the uniform density flow is unsteady and very close to the near- transition region as stated by Sakamoto & Haniu 1990 [44]. Samples of C Dx pro- ducedbythelowandhighresolutionsimulationsareshowninFigure 3.6a. Thetwo curvesarealignedforcomparativepurposes. Thedifferencesareinpeakmagnitudes and the higher frequency oscillations between the two cases. Overall, these two curves remain similar in quantitative and qualitative behavior. The time-averaged stream-wise drag coefficient in the standard resolution simulation is C Dx = 0.466 37 (a) (b) 0 50 100 0.42 0.44 0.46 0.48 0.5 0.52 C D 0 50 100 0.42 0.44 0.46 0.48 0.5 0.52 t C D 0 0.5 1 0 100 200 Power 0 0.5 1 0 50 100 150 St Power (c) (d) 5 10 15 10 -1 10 0 x U o 0 5 10 15 10 -2 10 -1 x u ′ o Figure 3.6: Re=1000, Fr=∞: (a) Sampled stream-wise drag coefficient. Top - low resolution,bottom-highresolution. (b)Powerspectraoftimehistoryofuatx≃ 5. Top-lowresolution,bottom-highresolution. (c)U o isthecenterline,meanstream- wise defect velocity. (d) u ′ o is the centerline r.m.s. stream-wise perturbation veloc- ity. (−) low resolution, (−−) high resolution, (−·−) Tomboulides & Orszag[59], (X) Wu & Faeth[63] at Re=930. 38 and the high resolution simulation produces a C Dx = 0.460, which is comparable to experimentally reported values of Roos & Willmarth 1971 [42]. Fourprobesareplacedinthewakeatx = 5.75andintheverticalandhorizontal planes about the wake centerline at a distance of r p = 0.6 in the high resolution simulations and x≃ 5,r p ≃ 0.675 in the low resolution ones. The power spectrum of u(~ x,t) is averaged between the probes and produces the spectrum shown in Figure 3.6b. The strongest frequency for both resolutions is at St= 0.191. This is the Strouhal number associated with the primary vortex shedding and the overall structure of the wake. Both simulations capture a shear instability frequency as reportedbyKim&Durbin1988[29]andSakamoto&Haniu1990[44]atSt =0.36, shown in Figure 3.6b. The high resolution case peak is more prominent, perhaps due to the position of the probe placements. Values of the mean stream-wise centerline velocity defect, U o , and root-mean- square stream-wise velocity perturbation component are reported in Figure 3.6c,d. There is some downstream variation between the high resolution and standard res- olution simulations in these quantities that is attributed to grid coarsening. The unsteady recirculation region length based on the inflection of U o about 1 is mea- sured to be L rr ≃ 2.25, which agrees with the prior experimental investigation of Linet al. 1992[35] andnumerical investigation ofTomboulides &Orszag2000[59]. There is some measure of error in defect velocity between the experiments of Wu & Faeth 1993 [63], the spectral results[59], and the high and standard resolution simulations along the centerline. Because there is no turbulence modelling used in the current simulations, both simulation resolutions exhibit behavior similar to an implicit large eddy simulation. 39 3.2.1 Downstream Distance from the Sphere in x and Nt For linearly stratified fluids, the relationship between the commonly used non- dimensional time, Nt≡N ∗ t ∗ , and downstream distance, x≡x ∗ /D ∗ , is given by: Nt = x(Fr/2) −1 (3.2) whereoriginsforxandtarelocatedatthespherecenter. AnalysisofthewakeinNt isarelative distance behindthesphere interms ofthefractionoftheFr-associated characteristic buoyant wavelength and portion of completed buoyant period. Consequently, comparisons between the flow of a uniform density fluid and stratified fluid cannot be properly expressed in terms of Nt since the wake of the uniform density fluid (i.e., N ∗ ≡ 0) cannot have an analytic description in Nt. In- tuitionsuggeststhatdescriptionoftheflowinthenearwakeshouldasymptotewith increasing Fr to the behavior of a uniform density fluid. Thus, x is the preferable choice for descriptions of the wake because the simulations explicitly contain the sphere, and physically relevant descriptions in the limit of Fr →∞ can be made. Notwithstanding, Nt remains useful in illustrating buoyant commonalities between the stratified unsteady wakes. Both x and Nt are used throughout remaining dis- cussions, and care is taken to ensure that descriptions in x and Nt correlate to the same physical location behind the sphere. For convenience, the equivalent distance and time behind the sphere are commonly given as a set {x,Nt}. 40 3.2.2 Averaging In The Near Wake Experimental studies of stratified wakes frequently use a digital particle image velocimetry (DPIV) setup, such as in Spedding 2001 [53] and Spedding 2002 [54], where a realization of the flowfield within a viewing window is spatially integrated: hqi xwc (y,z,t)= 1 ΔL Z xwc+ ΔL 2 xwc− ΔL 2 q(~ x,t)dx (3.3) where q is any flow field quantity, x wc is the mid-window downstream position, ΔL is the stream-wise width of the window. A temporal average over the course of several vortex shedding cycles will elim- inate variance of collected wake statistics from the simulation data, i.e.: q(x,y,z)= 1 ΔT Z t 1 t 0 q(~ x,t)dt (3.4) where ΔT is the representative time scale of averaging. The period used for com- puting q(x,y,z) is longer than that which is involved in the assumption that a spatial and time average are equivalent [18]. To highlight the differences between q(x,y,z) and hqi xwc (y,z,t), and the vari- ance of hqi xwc (y,z,t), a computationally generated unsteady flow field is analyzed using a DPIV setup [53, 54] for a Re = 1000, Fr = 4 flow. Figure 3.7a serves as a physical reference to the planes referred to in the DPIV arrangement, and Figure 3.8a serves as an illustrative reference to the DPIV setup and its subse- quent affect on statistical data when applied within the near wake. The visualiza- tion/analysis window in Figure 3.8a is centered at {x,Nt} = {6,3} and z = 0, of 41 (a) (b) Figure 3.7: Sketches (not to scale). (a) Vertical and horizontal center planes and wake centerline. (b) Galilean transformation of stream-wise velocity profiles in upper vertical plane and negative horizontal plane and wake length scales L σv and L σh . width ΔX = 3.5, and is coplanar with the y = 0 vertical center plane. The vi- sualization plane is integrated spatially to obtain an averaged stream-wise velocity profile [53]: U xwc (0,z,t) =hu(~ x,t)i xwc = 1 ΔX Z xwc+ ΔX 2 xwc− ΔX 2 u(~ x,t)dx (3.5) This average velocity is then assumed to be the average stream-wise velocity every- where within the visualization window. The root mean square perturbation values are then calculated spatially along the stream-wise direction[53]: hu ′ i xwc (0,z,t) = " 1 ΔX Z xwc+ ΔX 2 xwc− ΔX 2 [u(x,0,z,t)−U xwc (0,z,t)] 2 dx #1 2 (3.6) hv ′ i xwc (0,z,t) = " 1 ΔX Z xwc+ ΔX 2 xwc− ΔX 2 v(x,0,z,t) 2 dx #1 2 (3.7) hw ′ i xwc (0,z,t) = " 1 ΔX Z xwc+ ΔX 2 xwc− ΔX 2 w(x,0,z,t) 2 dx #1 2 (3.8) 42 (a) (b) x z Wake 0 5 10 15 −10 −5 0 5 10 −3 −2 −1 0 1 2 3 0 0.01 0.02 0.03 0.04 0.05 z w ′ Figure 3.8: Re = 1000, Fr = 4: (a) Shaded area is a DPIV win- dow centered about x wc = 6. The sphere at (x,z) =(0,0) is moving to the left. (−−) damping layer. (+) detected locations for |z| > 2 and 3.75≤x≤ 14.75 where w ′ (x,0,z,t). 10 −3 ≃ 0 over a ΔT = 40, sampled every 0.4 time units. (b) (−·−) Simulation comparison of hw ′ i xwc (0,z) to (N) Spedding [53] at Re = 5000. (a) (b) −8 −6 −4 −2 0 2 4 6 8 0 0.01 0.02 0.03 z w ′ −3 −2 −1 0 1 2 3 0 0.01 0.02 0.03 0.04 0.05 z w ′ Figure 3.9: Re = 1000, Fr = 4: (a) w ′ (~ x,t = 20) in vertical plane for (−•−) x =4.45, (−−) x =5.25, (−− w/N) x =6.05, (−− w/ △) x =6.85, (−− w/ ⋄) x = 7.65. (b) hw ′ i xwc (0,z,t) in vertical plane at (◦) t =0, (−) t =20, and (−·− w/ ⋄) t = 40 at x wc = 6;N Spedding [53] at Re = 5000 43 where x wc =6 is the x location about which the example analysis window (shaded area of Figure 3.8a) is centered in the stream-wise direction. Comparisons of the low resolution simulation are made against experimentally collected DPIV field data taken in the near wake at Re = 5000, Fr = 4 (Fig. 8 of Spedding [53]; reproduced in Figure 3.8b). In the stratified near wake, this is possibly the only in-plane stream-wise and vertical perturbation velocity data set available. Although the current simulation is performed at Re =1000, the uniform density case is a highly unsteady, near-transition flow, and comparisons between the statistically collected data sets remain of interest as wake behavior between Re = 1000 and Re =5000 are expected to be similar. At Fr = 4, the computationally obtained quantities in Eqs. (3.6)-(3.8) are spatially averaged as indicated and an additional average in time is applied to a series of windows of the same dimension centered about{x,Nt} ={6,3} such that at y = 0: hqi xwc (0,z) = 1 ΔT Z t 1 t 0 hqi xwc (0,z,t)dt (3.9) The interval of the time average is ΔT = 40 and is intended to remove any vari- ance in the collected statistics at that location. If time and spatial averaging are equivalent, or roughly equivalent, then there should be no significant penalty in comparison from the double averaging technique. Equation (3.8)assumesthatW(x wc ,0,z,t) =0,andallfluctuatingquantitiesin the vertical direction contribute tohw ′ i xwc (0,z,t) only. The computational results forhw ′ i xwc (0,z) are compared against the Re =5000, Fr =4 experimental data in Figure 3.8b. In the double-averaging technique, regions of non-negligible vertical perturbation velocities appear in Figure 3.8b, and hw ′ i xwc (0,z) extends out of the 44 wake region for |z|≥ 2, where the value of |z| = 2 is the current choice in demar- cating the wake and free stream regions. In that same Figure, the computed profile ofhw ′ i xwc (0,z) is also wider in z than the experimentally collected hw ′ i xwc (0,z,t). This collective difference is referred to as the ”tails” of the perturbation quantity data. They could suggest evidence of internal waves being generated outward from the unsteady wake motion, wave reflections from the finite numerical boundary, or some combination thereof. However, the origin of both the non-negligible regions of hw ′ i xwc (0,z) in Figure 3.8b for |z|≥ 2 and the width of the hw ′ i xwc (0,z) curve for|z| <2 are results of the averaging techniques being applied in the near wake. If W(x wc ,0,z,t) ≡ 0, the z location where w ′ (x,0,z,t) ≃ 0 (locations of lee wave crests and troughs as per linear theory [33]) will certainly change in x. This behavior is qualitatively shown in Figure 3.8a, and the lines of constant phase for lee waves are curved but not necessarily in accordance with linear theory [10]. Us- ing a temporally sampled subset of the data used to create Figure 3.8a, the results in Figure 3.9a for values of w ′ (x,0,z,t) show that there are indeed x locations where w ′ (x,0,z,t) ≃ 0 at a specific z location; occurring within Figure 3.9a at (z,x) = (±5.6,4.45) and (z,x) =(±3.2,5.25). In Figure in 3.8a), there are subse- quent downstream locations between the lee wave crests and troughs in x where w ′ (x,0,z,t) ≇ 0 anywhere in the far field between 6.05≤x≤ 7.65, denoted by the double dashed lines. These trends correlate with Figure 3.8a where between 8. x. 12 no w ′ (x,0,z,t)≃ 0 is detected. Because the width of the wake profile for hw ′ i xwc (0,z) is also wider than the experimental results presented in Figure 3.8b, the variation of the profile width of hw ′ i xwc (0,z,t) in time is also addressed. Three spatially averaged windows for hw ′ i xwc (0,z,t) at arbitrary times t = {0,20,40} are compared against the exper- imental data in Figure 3.9b. From the three discrete-time profiles, it is apparent 45 that the hw ′ i xwc (0,z,t) profile significantly varies in time, and at some times (e.g. t = 20), the computational results produce a profile width forhw ′ i xwc (0,z,t) com- parable to the experimental results. The fluctuations in the width of the profile forhw ′ i xwc (0,z,t) are related tothe overall unsteadiness ofthe near wake structure itself. If the DPIV window is spatially integrated and W(x wc ,0,z,t) =0, the majority of x locations contribute some non-zero quantity to hw ′ i xwc (0,z,t) in the free- stream region, thereby creating the tails’ non-zero asymptote with increasing|z| as seen in Figure 3.8b. The free-stream values for hw ′ i xwc (0,z,t) will be dependent on the width and the location of the DPIV averaging window and the time of the spatial integration. A similar argument can be made for profiles of v ′ in the NEQ regime (e.g., Fig. 9 of Spedding [53] at Nt = 9). When the assumption V(~ x,t)= W(~ x,t) =0 is removed, which is the case for all data presented after this Section, the asymptotic free-stream values ofv ′ (~ x,t)andw ′ (~ x,t)will approach zero by the onset of the damping layer. The computational results are presented using only temporal averaging. An inferred benefit of the analysis of the free stream perturbation quantities is that when w(~ x,t) is only averaged in time then there will be an expected free stream z locationwhereW(~ x)≃ 0fordownstream locationsofxthatcoincide withthecrest or trough of a lee wave, and this feature can be utilized (i.e. Figure 3.2b, Figure 3.5, Section 3.2.4) when identifying the crests or troughs of lee waves at the edge of the unsteady wake. The flow field is described using a Reynolds decomposition: q(~ x,t) =Q(~ x)+q ′ (~ x,t) (3.10) 46 where q =(u,v,w,ρ ′ ) represents the full-field ofthe primitive variables. Their time averaged values are Q = (U,V,W,ρ ′ ), and q ′ = (u ′ ,v ′ ,w ′ ,ρ ′ ) are the perturbation quantities. Notethatρ ′ (~ x,t)ofEqs. (1.12) and (1.7)isnowcomprisedofanaverage and perturbation component. The root mean square of a perturbation component is then defined: q ′ (~ x) = 1 ΔT Z t 1 to [q(~ x,t)−Q(~ x)] 2 dt 1/2 (3.11) where ΔT = 200 for the low resolution Re = 1000 cases. t o begins from a statis- tically steady state based on the time history of the stream-wise drag coefficient, and t 1 = t 0 +ΔT. For brevity, the root mean square perturbation components cal- culated by Eq. (3.11) are denoted by the prime only. Although the simulations are performed by considering the inertial frameas attached tothe sphere, the flow field is presented in the laboratory reference frame. This is accomplished by Galilean transformation in the stream-wise direction because the sphere travels in a steady, horizontal motion. A sketch of the transformation is in Figure 3.7b. 3.2.3 Numerical Resolution at Re=1000, Fr=4 In Section 3.2, it is shown that the high resolution and standard resolution sim- ulations exhibit similar characteristics at Re = 1000 and Fr = ∞. That case represents the upper limit of the stratified regime of interest. Above Fr = 4, the near wake behavior maintains its highly unsteady structure, and a Fr = 4 rep- resents the lowest limit of this flow regime of interest. Understanding the effects of the difference in numerical resolution on the range Fr = [4,∞], and how the statistical trends are affected by this difference are addressed. 47 (a) (b) Figure 3.10: Re=1000, Fr=4: Instantaneous isosurface of λ 2 -criterion, defined by λ 2 =−0.001. (a) Low resolution, (b) high resolution. The grey region is the z =0 plane. In these stratified simulations, it is challenging to quantitatively separate the internal waves at the egde of the wake from the wake structure itself using the λ 2 -criterion of Jeong & Hussain 1995 [27] as vorticity is generated by both inter- nal waves and vortical wake structure. Regardless, based on the isosurface of the λ 2 -criterion in Figure 3.10a, the low resolution simulations are able to capture the largest scale structures present in the stratified near wake. In qualitative compar- ison with Figure 3.10a, the higher resolution simulation in Figure 3.10b captures higher detail of the instantaneous near wake structure. Comparison of the resultant centerline perturbation quantities, denoted by an ”o”, between the low and high resolution simulations are shown in Figure 3.11a, andtheperturbationquantitiesattheverticaledgeofthewake(definedbyL σv )are presentedinFigure 3.13b. Theperturbationquantitiesareconsistentlylargerinthe high resolution simulation when compared against the low resolution simulation. The perturbations in u ′ are larger in the high resolution simulation, but there is also an increased centerline velocity defect, U o , seen in Figure 3.12 over the low resolution values. The magnitudes for quantities v ′ , w ′ , and ρ ′ also trend similarly 48 (a) (b) 0 5 10 15 10 -2 10 -1 x u ′ o 0 5 10 15 10 −2 10 −1 u ′ L σv x 0 5 10 15 10 -2 10 -1 x v ′ o 0 5 10 15 10 −2 10 −1 v ′ L σv x 0 5 10 15 10 −2 10 −1 x w ′ o 0 5 10 15 10 −2 10 −1 w ′ L σv x 0 5 10 15 10 −2 10 −1 x ρ ′ o 0 5 10 15 10 −2 10 −1 ρ ′ L σv x Figure 3.11: Re=1000, Fr=4: Perturbation quantity comparison (a) on centerline, (b) at z = L σv . (−) Standard resolution, (−−) High resolution. 49 5 10 15 10 -1 10 0 x U o Figure 3.12: Re=1000, Fr=4: Centerline velocity defect.(−) standard resolution, (−−) high resolution. in the near wake along both the centerline and wake edge. It is of note that the peak magnitude of w ′ o at w ′ at L σv for x = 6 in the high resolution simulations is similar totheexperimental values obtainedby Spedding[53]inFigure 3.8battheir respectivelocations. Significantdeparturesinqualitativebehaviorbetweendifferent resolution simulations occur around or after the near wake value of Nt ≃ 2−3, where the grid resolution for the low resolution simulation starts to coarsen away from the near wake region. Later sections will also focus on a data collapse description of the near wake as a function of Fr and L σv . As such, a comparison between the length scales cal- culated between the high resolution simulation and standard resolution simulation is provided in Figure 3.13. Critically, the vertical length scales are nearly identi- cal between the high and standard resolution case until a downstream position of x ≃ 8 or Nt ≃ 2−3. This range fully encompasses the region considered to be the near wake of the sphere. Thus, the low resolution simulations are expected to have sufficient resolution to capture the relevant near wake physics for a range of 50 0 5 10 15 0 0.5 1 1.5 2 x L σ Figure 3.13: Re=1000, Fr=4: Computed half-widths of the near wake at (−) low resolution and (−−) high resolution. No symbol is vertical length scale, (N) is horizontal length scale. stratified fluids between Fr = 4 and Fr = ∞ based on the parity of calculated vertical length scales and the agreement of quantitative and qualitative behavior both at the wake centerline and its edges. The low resolution simulations are used in the remainder of the discussion of the near wake behavior. 3.2.4 Length Scales of the Near Wake The horizontal and vertical half-widths of the wake are the definitional choice of wake length scales. The average stream-wise velocity component along the cen- terline U(x,0,0) = U o (x) = U o , called the defect velocity, is used to define these length scales. The edges of the wake are defined by locations where U(~ x) =0.15U o as sketched in Figure 3.7b. The horizontal half-width of the wake is denoted by L σh ≡ L ∗ σh /D ∗ , and the vertical half height is denoted as L σv ≡ L ∗ σv /D ∗ . Physical locations where |y|<L σh and|z| < L σv are referred to as inside the wake region, 51 0 5 10 15 0 0.5 1 1.5 x L σ Figure 3.14: Re=1000. Half-widths of the near wake. Open symbols correlate to L σh , closed symbols to L σv . (◦,•) Fr=1, (,) Fr=2, (N) Fr=4 for comparison at x = 0.6, (−) Fr =∞. and outside of these bounds is referred to as the free stream, which may or may not contain lee waves of significant amplitude. When Fr <4, the simulations indicate that the stratification significantly sup- presses unsteady motion; observed experimentally [10] at Re = 2000. At Fr = 2, the flow is essentially steady, and when Fr = 1, the flow is certainly steady, with stratification having completely suppressed the statistical fluctuations in the flow. Significant lee waves are also present in the free stream for Fr < 4 as inferred through Figure 3.5d,e. The lee waves have a pronounced effect on the length scale behavior of the near wake forthe Fr ={1,2}cases as indicated by Figure 3.14. Initial values of L σv are significantly reduced for Fr ={1,2} when compared with the Fr = 4 case. After an initial reduction of L σv for both Fr ={1,2}, the vertical length scale continues to increase in some average sense with x. The vertical flow over the sphere in the Fr = 1case islikely suppressed duetoenergyrestrictions, butL σv between thetwo 52 cases scales similarly asdownstream locationincreases in x. The oscillations in L σh for the Fr =1 case (Figure 3.14) do not match the Fr = 2 case when x& 1.5 even though the widths begin at similar values. This is attributed to the significance of the lee wave influence, and inspection of Figure 3.14 shows local maxima and minima for Fr = 1 length scales that are out of phase directly behind the sphere. The oscillatory behavior of L σv is correlated with the influence of the average perturbation density, ρ ′ . From Eq. (1.11), 4ρ ′ /Fr 2 is representative of the average vertical bodyforce onthe fluid. The correlation between 4ρ ′ /Fr 2 andL σv is shown through Figure 3.15a,b, collected at positive z = L σv , the vertical edge of the wake, where minima and maxima are in phase with the minima and maxima of L σv in the near wake. As 4ρ ′ /Fr 2 decreases at the edge of the wake with increasing Fr (e.g. Figure 3.15c,d), the influence of the average body forcing on the vertical length scales in the near wake decreases. The lower Fr limit to obtain a fully three dimensional near wake region [10] occurs at Fr≃ 4, and this claim is supported by the near wake behavior of L σv in Figure3.16aandtheunsteadywakecontoursofw(~ x,t)inFigure 3.5f. Byinspection ofFigure3.15cfor2.x. 7,theaveragedensityisinphasewithL σv ,butthebody force is both reduced in magnitude and appears to slip out of phase with L σv , as downstream distance increases. In Figure 3.16a, each case of Fr ≥ 4, L σv grows in x after x = 2.25 until {x,Nt}≃{3Fr/2,3}. ForFr& 4,thelocationofx≃ 2.25isthecommonlocation where it is possible to parameterize the wake widths. The near wake region begins at this location, and the region where x. 2.25 is referred to as the recirculation region, which is also unsteady. The length of the region is denoted by L rr and is the location at which U o (L rr )≃ 1. The length of this region is only weakly dependent on Fr as L rr ≃ 2.25 for Fr =∞ and L rr ≃ 2.1 for Fr = 4. No attempt 53 (a) 0 5 10 15 −1 −0.5 0 0.5 1 4 F r 2 ρ ′ 0 5 10 15 0 0.25 0.5 0.75 1 L σ v (b) 0 5 10 15 −0.5 −0.25 0 0.25 0.5 4 F r 2 ρ ′ 0 5 10 15 0 0.25 0.5 0.75 1 L σ v (c) 0 5 10 15 −1 −0.5 0 0.5 1 4 F r 2 ρ ′ 0 5 10 15 0 0.25 0.5 0.75 1 L σ v (d) 0 5 10 15 −0.02 −0.01 0 0.01 0.02 4 F r 2 ρ ′ x 0 5 10 15 0 0.375 0.75 1.125 1.5 L σ v Figure 3.15: Re=1000. Comparison of average body force located at z = L σv , y = 0. (−−) 4ρ ′ /Fr 2 for each case and (−) is the zero axis line of the ρ ′ axis. Symbols indicate lines of L σv . (a) Fr = 1, (b) Fr =2, (c) Fr = 4, (d) Fr =10. 54 to parameterize this region is made because downstream characteristics in length scale behavior and flow field statistics are not compatible with the scaling form of parameterizations (Section 3.2.5) available for x& 2.25. At Fr =∞, a fitted parameterization of L σv and L σh in the near wake between 2.25≤x≤ 15 is: L σv∞ = L σh∞ = 0.37x 0.55 2.25≤x≤ 15 (3.12) where the subscript ∞ denotes the case Fr = ∞. Parameterizations of L σv and L σh in stratified near wakes should asymptote to Eq. (3.12) as Fr→∞. After the recirculation region, L σv of each lower-Fr case grows at a reduced rate from the Fr =∞ case. Based on the decreased growth rate of L σv in Figure 3.16a, for Fr ≥ 4, accounting for stratification will affect the exponent in Eq. (3.12), where L σv ∼x 0.55f(Fr) such that f(Fr) = 1 for Fr → ∞ and f(Fr) ≃ 0 when Fr = 4. Stratification will affect the amplitude of L σv at x≃ 2.25 due to suppression of the flow over the sphere as Fr decreases. A numerically fitted parameterization for L σv is given by: L σvFr = A(Fr)x 0.55 0.94+3.2/(Fr−4) 2.25.x. 3 2 Fr (3.13) valid for Fr≥ 4 where the subscript Fr of L σv indicates the Fr-specific case. It is worth noting that the description of A(Fr) is sensitive to the function f(Fr), and is empirically found to be: A(Fr)= 0.11sech[0.24(Fr−7.42)]+0.37 (3.14) 55 (a) (b) 0 5 10 15 0.6 0.8 1 1.2 1.4 1.6 1.8 x L σ v 0 5 10 15 0.6 0.8 1 1.2 1.4 1.6 1.8 x L σ h Figure 3.16: Re=1000: Half-widths of the near wake. (a) Vertical center plane at y = 0. (b) Horizontal center plane at z = 0 (every other Fr case omitted for clarity). (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞. (−−) Plotted through Nt = 3.5, Eq. (3.13) in (a), Eq. (3.15) in (b). Eq. (3.13) displays an asymptotic behavior towards Eq. (3.12) as Fr → ∞. The lower limit of Fr =4 is chosen because if Fr. 4, L σv does not scale similarly to the Fr≥ 4 cases by comparison between Figures 3.14 and 3.16a. The change in growth of L σv for{x,Nt}≥{3Fr/2, 3} is the beginning of the NEQ regime and is outside ofthe intended parameterization range. The agreement between Eq. (3.13) and the simulation data is qualitatively shown in Figure 3.16a. Figure 3.16b suggests that growth rates of L σh only weakly depend on Fr in the near wake. Again, x≃ 2.25 is a common starting point for significant growth in L σh . When Fr≥ 4, and using Fr =4 as the illustrative case, L σh does show an initial growth and changes its growth rate around{x,Nt}≃{3Fr/2, 3}. Opposite the effects of Fr on L σv , L σh is initially wider at x ≃ 2.25 with decreasing Fr. All L σh for Fr ≥ 4 begin at similar values directly behind the 56 sphere as, even for strongly stratified flows, horizontal motion around the sphere is not inhibited by potential energy requirements. After the recirculation region, at x≃ 2.25, initial values of L σh become dependent on Fr. As with L σv , L σh should approach Eq. (3.12) in the Fr → ∞ limit of the uniform density case. If growth rate of L σh in the near wake is only weakly dependent on Fr, then Fr affects only the starting width of the wake at x& 2.25, and a numerically fitted description is given of the form: L σhFr = B(Fr)x 0.55 2.25. x. 3 2 Fr (3.15) where: B(Fr)= 0.46exp(−0.31Fr)+0.37 (3.16) when Fr ≥ 4 and {x,Nt}≤{3Fr/2,3}. The parameterization for L σh is plotted against the simulation results in Figure 3.16b. The parameterization generally describes values of L σhFr , and deviations are related to the choice of the constant exponent in Eq. (3.15). Errors in estimation from the data because of this choice are typically less than 10%. 3.2.5 The Near Wake Centerline The centerline defect velocity, U o , in the near wake develops with downstream distance as U o ∼x −c in a uniform density fluid wake. The exponent c is commonly constant, O(1), and is sensitive to the range and behavior of the data downstream of the sphere. The exponent is sensitive to the range and the methodology of the fit, where directly behind the recirculation region, U o may behave as U o ∼x −2 , and fartherdownstream, self-similardevelopment isexpectedasU o ∼x −2/3 . Inthenear 57 wake, expectations of truly self-similar scaling should be abandoned because the wake is not self-similar. Nevertheless, scaling descriptions are commonly used and remain useful in describing the downstream behaviors of flow quantities. At Fr = ∞, the simulation reproduces an expected decay in U o , available in Figure 3.17a, and a simple scaling description is used to parameterize the defect velocity of the uniform density fluid wake: U o =3.58x −1.64 2.25≤x≤ 15 (3.17) where Eq. (3.17) is numerically fitted from the simulation data available from x≃ 2.25 to x =15. In Figure 3.17a, Fr = 4 is representative of stratification effects because U o departs from the uniform density case immediately downstream of x ≃ 2.25. Up- stream of x≃ 2.25, which is within the recirculation region, U o is indiscernible be- tween the Fr cases. Between {x,Nt}≃{2.25,1.1} and {x,Nt}≃{7,3.5}, depar- tures of the stratified wakes from the uniform density case are not severe compared to the departure further downstream. After {x,Nt}≃{7,3.5} U o increases very slightly, and after{x,Nt}≃{10,5}, it returns to decay. Similar behavior has been reportedbyBonnier&Eiff2002[5]forFr ={3,6,10}atRe ={3400,6900,11500}, respectively. At the constant Re = 1000, current results only indicate downstream increase of U o in the Fr = 4 case, and the severity of this departure from the Fr =∞ case lessens with increasing Fr as neither Fr = 5 or Fr = 6 exhibit this increase for the same time in Nt. Spedding et al. [55] notes that mean stream-wise velocities in the NEQ and far wake regions can be greater than that of high Fr by nearly an order of magnitude, and it is clear by maintaining a constant Re and 58 changing Fr, as in Figure 3.17a, that this is caused by stratification effects in the near wake. It is shown in Section 3.2.4 that stratification influences the growth of the ver- tical wake width in the near wake. A stream-wise, mean local Froude number is defined as Fr Lσv = U ∗ o /(N ∗ L ∗ σv ). Like the internal Fr, Fr Lσv is representative of the mean energy available to a fluid particle on the wake centerline as compared to the potential energy required to vertically displace that same particle to the edge of the wake. If U o is normalized by L σv , then that ratio may be manipulated: U o L σv = U ∗ o U ∗ s D ∗ L ∗ σv = U ∗ o N ∗ L ∗ σv N ∗ D ∗ U ∗ s = 2 Fr Lσv Fr (3.18) By Figure 3.17b, regardless of Fr, U o /L σv collapses completely onto the uniform density case when {x,Nt}.{5,10/Fr}. Downstream of this point, the stratified casesbegintoslightlydivergefromtheuniformdensitycaseuntil{x,Nt}≃{Fr,2} whereafter U o /L σv severely departs from the uniform density case. The data collapse in Figure 3.17b indicates that in the stratified fluid, the localratioofavailable meankinetic topotentialenergy requirements directly scales with internal Froude number, Fr. When 4 ≤ Fr ≤ ∞ and Nt . 2 this ratio is independent of the stratification level. Additional interpretation of U o /L σv in the near wake is not immediately clear, although U ∗ o /L ∗ σv is proportional to some wake vorticaltimescale,andU o /L σv isthattimescalenormalizedbytheinverseadvective timescale U ∗ s /D ∗ ; previously utilized in studies of stratified wakes [53]. Values of U o /L σv in cases of Fr = {4,5,6,8} significantly transition away from the uniform density case limit at Nt ≃ 2 and appear to reach a new, sim- ilar rate of decay at corresponding values x≃{8.4,9.0,9.9,12} or, equivalently, Nt≃{4.2,3.6,3.3,3}; indicating a Fr-dependence on the length of this transition. 59 (a) (b) 0 5 10 15 10 -1 10 0 x U o 0 5 10 15 10 -1 10 0 x U o L σ v Figure 3.17: Re=1000. (a) Centerline defect velocity, U o . (b) Normalized U o /L σv . (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞. (−−) Eq. (3.17) in (a), (−−) Eq. (3.19) in (b). The common departure at Nt =2 of U o /L σv indicates that the mean flow starts to transitionintotheNEQregimeearlierthanL σv alonemightindicate(Section 3.2.4). A straightforward modification to Eq. (3.17) that can parameterize U o in the near-wake for Fr ≥ 4 is found by using U o /L σv , applied downstream of the recir- culation region. The general collapse of the Fr ≥ 4 data in Figure 3.17b onto the Fr =∞ data for Nt. 2 is used to generalize the curve fit for any Fr ≥ 4. After Nt≃ 2, a virtual origin description in lieu of Eq. (3.17) for the NEQ region would require a shift in origin on a per-Fr basis. A modified parameterization, for the near wake, consistent with Eq. (3.17) is: U oFr L σvFr =9.68x −2.19 2.25.x.Fr (3.19) 60 where the upper limit in x is equivalent to Nt =2; valid for Fr≥ 4 and x≤ 15. A largerdomaincouldincrease theaccuracy ofthisstatement forthehigherFr cases, where there is notable departure from the description provided by Eq. (3.19). Descriptions of the perturbation components are presented in two ways. First, the unmanipulated data is available for all perturbation components for direct comparison. Second, the collapsed parameterized data is presented in terms of physical downstream distance x. In the cases of w ′ o and ρ ′ o , the data in Nt is used as a guide in an attempt to collapse the data of the near wake. Centerlineperturbationvelocitiesfortheuniformdensitycasebegininaslightly anisotropic state for the Re = 1000 flow. Peak perturbation velocity values for ~ u ′ o =(u ′ o , v ′ o , w ′ o )= (0.159, 0.166, 0.156) at x = (2.2,2.5,2.5) from Figure 3.18. While not beginning with the same energy distribution, each perturbation compo- nent in the uniform density case decreases at a similar rate. This behavior is in contrast to the behavior of u ′ o , v ′ o , and w ′ o for the stratified Fr = 4 case. Quali- tatively, it is also clear that the stratified perturbation velocities in the near wake do not scale in a similar fashion to each other whatsoever. Parameterizations of velocity perturbations in uniform density case are: u ′ o = 0.32x −0.92 2.25≤x≤ 15 (3.20) v ′ o = 0.38x −0.86 2.25≤x≤ 15 (3.21) w ′ o = 0.38x −0.94 2.25≤x≤ 15 (3.22) As evident in Figure 3.19a, stratification immediately affects the magnitude of u ′ o in the recirculation region and subsequent values in the near wake, but it does not appear to generate a pronounced effect on the decay rate of u ′ o . For Fr = 4, the magnitude of u ′ o at x = 2.2 is roughly two-thirds that of the uniform density 61 0 5 10 15 10 -2 10 -1 x ~ u ′ o Figure3.18: Re=1000. Anisotropyof(−)u ′ o ,(−·−)v ′ o ,(−−)w ′ o inthenearwake. Unmarked lines are uniform density case. The Fr =4 case: (N) u ′ ,(•) v ′ ,() w ′ . case. Only the Fr = 4 case hints at any significant departure from the decay rate of u ′ o at{x,Nt}≃{7,3.5}. If an Ozmidov length scale in the near wake is defined as: l ∗ o = ǫ ∗ N ∗ 3 1/2 = u ′ o ∗3 L ∗ σv N ∗ 3 ! 1/2 (3.23) A manipulation of Eq. (3.23) proceeds in the following fashion: l ∗ o L ∗ σv = l o L σv = (u ′ o U ∗ s ) 3 (L σv D ∗ ) 3 N ∗ 3 ! 1/2 = u ′ o L σv Fr 2 3/2 (3.24) and after rearrangement of terms provides: u ′ o L σv = 2 Fr l o L σv 2/3 (3.25) Eq. (3.25)indicates thatthenon-dimensionalquantity u ′ o /L σv isrepresentative of the ratio of the the buoyant overturn length scale to turbulent integral length 62 (a) (b) 0 5 10 15 10 -2 10 -1 x u ′ o 0 5 10 15 10 -2 10 -1 x u ′ o L σ v Figure 3.19: Re=1000. Centerline stream-wise perturbation velocity (a) unmod- ified, (b) normalized. (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞. (−−) Eq. (3.20) in (a), Eq. (3.27) in (b). scale [54]. If Fr → ∞ and l o /L σv ∼ Fr 3/2 , Eq. (3.25) remains consistent in the unstratified limit. In a uniform density fluid, u ′ o /L σv = (u ′ o ∗ /L σv ∗ )(D ∗ /U ∗ ) is representative of an eddy turnover time normalized by the advective time scale of theflow. Ifalocalstream-wiseFroudenumberisdefined[54]alongthecenterline as Fr l1 =u ′ o ∗ /(N ∗ L σv ∗ ), thennormalizingu ′ o byL σv is alsoequivalent tonormalizing the local Froude number by the internal Froude number, u ′ o /L σv =u ′ o D ∗ /(U ∗ s L σv ) and subsequently u ′ o /L σv =2Fr l1 /Fr, which is also consistent in the Fr → ∞ limit, and the localperturbation kinetic energy available in the near wake will scale with Fr. Once the normalization u ′ o /L σv is made, in Figure3.19b, the near wake stream- wise perturbation velocities collapse onto the data of the uniform density case, confirming that that the local stream-wise kinetic energy scales with the potential energy requirements of particles passing over the sphere. The kink in the decay rate of u ′ o /L σv remains perceptible, yet under dramatic, as the flow enters the 63 NEQ region. Generally, the higher the Fr, the more closely this collapse appears to follow the uniform density case for downstream values in x. There are some notable departures for the Fr. 8 cases of this behavior further downstream, but this may be attributed to the entry of the wake into the NEQ regime as well as edge-of-domain proximity of the damping layer. At sufficiently high Re, if u ′ o (or v ′ o or w ′ o ) in the uniform density case, and both L σv∞ and L σvFr are known (i.e., by shadowgraph), then the value of the stream-wise perturbation velocity in the near wake can be predicted by: u ′ oFr =u ′ o∞ L σvFr L σv∞ 2.25.x. 3 2 Fr (3.26) valid for Fr ≥ 4, where subscripts Fr and ∞ indicate the stratified and uniform density case values, respectively. If this data is not available then a consistent parameterization for u ′ o in the near wake can be obtained by a fit of the uniform density case for u ′ o /L σv : u ′ oFr L σvFr = u ′ o∞ L σv∞ = 0.87x −1.49 2.25.x. 3 2 Fr (3.27) and is valid for Fr ≥ 4. Because of the u ′ o /L σv collapse in Figure 3.19b and Eq. (3.26), the subscripts of ∞ and Fr are interchangeable, and Eq. (3.13) may be used in place of L σv for a full parameterization in terms of x or Nt only. Values of v ′ o in Figure 3.20a indicate that stratification can affect the decay rateofcenterline perturbationvalues inthenearwakewithin shortdistances down- stream of the sphere. In the cases Fr ={4,5}, there is a notable change in decay where v ′ o begins to level off further downstream from the sphere. This occurs at x ≃ {6,7.5} or equivalently Nt ≃ 3 for both Fr = 4 and Fr = 5. Following the 64 (a) (b) 0 5 10 15 10 -2 10 -1 x v ′ o 0 5 10 15 10 -2 10 -1 x v ′ o L σ v Figure 3.20: Re=1000. Centerline horizontal perturbation velocity (a) unmod- ified, (b) normalized. (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞. procedureappliedtou ′ o ,v ′ o isrescaledbyL σv inFigure3.20bsuchthatatx≃ 2.25 the values of the uniform density case and stratified cases for v ′ o /L σv reach a closer agreement. Similar to u ′ o , this is equivalent to normalizing a lateral, local Froude number Fr l2 =v ′ o ∗ /(N ∗ L ∗ σv ) by the internal Froude number. This rescaling, seen in Figure 3.20b, is insufficient to collapse the data in the near wake region for v ′ o . An attempt at correction to the collapse is made by the use of a fitted asymptotic function such that parameterizations of v ′ o , for x& 2.25, scale with the uniform density case. The fitted form of the asymptotic function, f v (Fr) is: f v (Fr)=0.26tanh 0.39 Fr 2 (3.28) and asymptotes to 1 as Fr → ∞. It is reiterated that the parameterizations are intended tobeapplied forFr≥ 4, which isalso inherent tothe parameterizationof Eq. (3.27). The behavior of Eq. (3.28) as Fr→ 0 is not necessarily representative 65 0 5 10 15 10 -1 10 0 x 1 f v v ′ o L σ v Figure 3.21: Re=1000. Collapse of v ′ o /L σv by Eq. (3.28). (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞, (−−) Eq. (3.29). of the wake behavior. The parameterization that is consistent with Eq. (3.21) becomes: v ′ oFr L σvFr =3.95f v (Fr)x −1.43 2.25.x. 3 2 Fr (3.29) and Fr has replaced the ∞ subscript in Eq. (3.29) through the same reasoning applied to u ′ o . This parameterization is shown in Figure 3.21 where the near wake values and Eq. (3.29) become indistinguishable. Out of all perturbation velocity values along the centerline, w ′ o is the most dramaticallyaffectedbystratification. Figure3.22aindicatessignificant,oscillatory behavior in w ′ o within short downstream distances from the sphere, even at a moderate values of Fr, such as Fr = 8. Furthermore, the curvature of the initial decay for lower Fr is inflected from that of the u ′ o or v ′ o quantities. It appears from Figure 3.22a that w ′ o will always depart from the uniform density case for low-to-moderate Fr values at some x and Fr pairing. The first visible valley of w ′ o occurs at locations x≃{6,7.5,9,12} for Fr ={2,2.5,3,4}, respectively, or at Nt≃ 3 for all these Fr cases. 66 (a) (b) 0 5 10 15 10 −2 10 −1 x w ′ o 0 2 4 6 8 10 −2 10 −1 10 0 (Nt) 0 .9 4 w ′ o ! F r 2 " 0 .9 4 (c) (d) 0 5 10 15 10 −1 10 0 x 1 f w w ′ o L σ v 0 5 10 15 10 −2 10 −1 x w ′ o Figure 3.22: Re = 1000; (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16, (−) Fr=∞. (a) Highly oscillatory nature of w ′ o , (−−) Eq. (3.22), (−·−) Eq. (3.30). (b)Collapseofw ′ o inNt, (−−) ∼Nt −0.94 , (−·−) ∼f w Nt −0.94 . (c) Reduced oscillations in perturbation velocity. (d) Pa- rameterized vertical perturbation velocity, (−·−) Eq. (3.33). 67 0 2 4 6 8 10 −1 10 0 10 1 Nt 2 F r ρ ′ o w ′ o Figure 3.23: Re=1000,Fr=[4,16]: Ratio of centerline density perturbation ρ ′ o to vertical velocity perturbation w ′ o . An estimate of Fr value where this oscillatory behavior could be considered negligible is possible. The valleys are described by: w ′ o | Nt=3 = .033x −0.54 (3.30) forFr≥ 4,wherew ′ o | Nt=3 indicatestherecordedvalueofw ′ o inthefirstvalley. The Fr = 10 case is not included in the parameterization due to the Nt = 3 condition beingcoincidentwiththeonsetofthedampinglayer. TheintersectionofEq. (3.22) and Eq. (3.30) would occur at x ≃ 450. Based on the visible trend, the case at which stratification becomes negligible in the near wake would also occur within a negligible valley at Nt≃ 3. This negligible reduction correlates with a Fr ≃ 150. Naturally, this prediction is sensitive to the accuracy of the description in both Eq. (3.22) and Eq. (3.30). With a larger domain, it may be possible to obtain a better estimate, but perhaps it is reasonable to suggest that when Fr ∼ O(100), the differentiation between astratified flow andauniform density flow is negligible. 68 The valleys of w ′ o at Nt = 3 are located at the same location of the local maximum of the collapsed ρ ′ o /w ′ o at Nt = 3 in Figure 3.23. Potential energy in a linearly stratified fluid is associated with ρ ′ o because ρ ′ o is representative of the vertical displacement of fluid particles away from their equilibrium position. The vertical turbulent kinetic energy is associated with w ′ o . Thus, Figure 3.23 can be seenasrepresentativeofthecenterlineratioofturbulentpotentialenergytovertical turbulent kinetic energy. The valleys in w ′ o of Figure 3.22a are representing the process of transferring vertical turbulent kinetic energy into turbulent potential energy. As the values of w ′ o increase out of the valleys, some of the turbulent potential energy is being converted back into vertical turbulent kinetic energy and vice versa for the descent of w ′ o into a valley. Both u ′ o and v ′ o behave as the uniform density case when Fr → ∞, and Eqs. (3.27) and (3.29) reflect this condition. Any parameterization of w ′ o should asymptote in a similar fashion. Oscillations are not likely removed by normalizing w ′ o by L σv alone because, as seen in Figure 3.16a, there are no notable oscilla- tions of L σv for Fr ≥ 4 that would account for such large amplitude oscillations in w ′ o . Oscillations could be smoothed if the quantity w ′ o is normalized by the ”correct” choice of oscillatory function. Because of the form chosen to describe w ′ o in Eq. (3.22), a parameterized equation is sought in the form: w ′ o Fr 2 0.94 =f w (Fr,Nt)(Nt) −0.94 (3.31) where f w (Fr,Nt) is some oscillatory function dependent on Fr and Nt that com- pliments thedecayrateintimeandshould asymptotetounitywhen Fr→∞. The final form of Eq. (3.31) is expected to translate consistently into x, and because f w asymptotes to unity as Fr → ∞, w ′ o must be multiplied by (Fr/2) −0.94 , as done 69 in Figure 3.22b. Like f v , there is no particular restriction on f w as Fr→ 0 because Fr = 4 is intended to be the limiting case. When viewed in Nt through Figure 3.22b, w ′ o (Fr/2) 0.94 is collapsed and the decay appears independent of Fr for 4.5/Fr . Nt ≃ 3, at which point there is growth in Nt until Nt& 4 when there is a return to decay. The proposed form for f w is then: f w (Fr,Nt) =0.15 h 2+sin 1.16 π 2 (Nt) 0.94 i (3.32) A description of w ′ o /L σv in the Fr =∞ case provides an estimate of the strat- ified behavior of w ′ o /(f w L σv ) by Figure 3.22c. A parameterization that describes w ′ o /L σv at Fr =∞ is then: w ′ oFr L σvFr = w ′ o∞ L σv∞ = 2.9f w (Fr,x)x −1.49 2.25≤x≤ 3 2 Fr (3.33) which is a parameterization of the vertical velocity coefficient and when combined with Eq. (3.13) may be used for a parameterization of w ′ o purely in x. This pa- rameterization is shown in Figure 3.22d where Eq. (3.33) tends to over or under predict the Fr cases controlled by the generated collapse seen in Figure 3.22c. Centerline values of ρ ′ o for the near wake are presented in Figure 3.24, where views of ρ ′ o show a relative independence from Fr for x . 2.25, which is within the recirculation region. The development of ρ ′ o is physically described by a quick growthwithin and behind therecirculation region, followed by anoscillatory decay. ThisisrelatedtothetransferofenergyaspreviouslydescribedthroughFigure 3.23. When viewed in x, there is no obvious simple correlation between the value of the magnitudes and location ofthe peaks between the various Fr cases. Analysis of ρ ′ o in Nt, as in Figure 3.24b may aid in the parameterization of ρ ′ o because it appears 70 (a) (b) 0 5 10 15 10 −2 10 −1 x ρ o ′ 0 2 4 6 8 10 −2 10 −1 Nt ρ ′ o Figure 3.24: Re = 1000. Centerline density perturbation. (N) Fr=4, () Fr=5, (+) Fr=6, (×) Fr=8, () Fr=10, (⋆) Fr=16. that Nt ≃ 5 may be the time at which ρ ′ o is at a local minimum for all cases of Fr≥ 4. Unfortunately, domainlimitsprohibitutilizingthetrends ofthehigherFr cases to generate an oscillatory non-linear parameterization of ρ ′ o in Nt alone, as at- tempted for w ′ o . Results for the Fr = 4 case in the near wake leads to use of an empirical nonlinear parameterization in Nt of the form: ρ ′ o = 0.116 1+2.4 sin 2 (0.32Nt) (Nt) 1.26 exp(−0.64 (Nt)) 0.3≤Nt. 4 (3.34) which isincluded inFigure3.25againstthesimulation resultofρ ′ o forFr = 4. The result of this parameterization is quite representative of the behavior of ρ ′ o in the near wake and the form of Eq. (3.34) is chosen based on the behaviors of ρ ′ o across the various cases in Fr in Figure 3.24b. The severity of the oscillations appears to decrease based on the Fr ={4,5,6} cases at Nt≃ 5 in Figure 3.24b, although the Nt = 5 condition for the Fr = 6 case is located at the onset of the damping layer. 71 0 2 4 6 8 10 −2 10 −1 Nt ρ ′ o Figure 3.25: Re = 1000. Comparison of simulation result to parameterization of ρ ′ o . (N) Fr=4, () Fr=5, (⋆) Fr=16, (−−) Eq. (3.34), (−·−) Eq. (3.35). The tendency towards a decrease in severity in the oscillations of ρ ′ o suggests that, in the near wake, the sinusoidal term of Eq. (3.34) may be removed for a higher Fr fluid. For Fr≥ 5, the description in x can be empirically fitted as: ρ ′ o = 0.056tanh Fr 6 x 17.5Fr −1 exp −15Fr −2 2.25≤ x. 3 2 Fr (3.35) and this is also shown, translated into Nt, in Figure 3.25 for the Fr = 5 and Fr = 16 cases. For Fr ≥ 5, Eq. (3.35) appears as a good, simple alternative to the form of Eq. (3.34) for the higher Fr cases, although for the Fr → ∞ case, Eq. (3.35) results in a negligible value. This is not of concern as the buoyant body force in Eq. (1.11) is identically zero as Fr→∞. 3.2.6 TheNearWakeintheVerticalandHorizontalPlanes The near wake of the sphere at Re = 1000 for moderate to high Fr is fully three- dimensional and anisotropic. A description of the near wake in the vertical and 72 horizontal planes is useful in completing descriptions ofthenear wake. Because the gravityvectorisimplicitlypresentas−g ˆ k inEq. (1.11),andthesphereisanaxially symmetric body, there are two planes about which the statistical descriptions are expected to be symmetric, which are the planes shown in Figure 3.7a. Theverticalcenterplaneisaplaneofsymmetrybecausegravitydoesnotcontain a horizontal bias in Eq. (1.11), and the horizontal centerplane is a plane of symme- try because particles should flow equally over and under the sphere because, in a linearlystratifiedfluid,thereisnoverticalbiasinthethebuoyantbodyforcesabout the horizontal centerplane. Unless otherwise stated, all references to the vertical plane and horizontal plane refer to the respective centerplanes of symmetry. At Re =1000, the uniform density fluid flow around the sphere and the wake is assumed axisymmetric such that variance in the statistical description of the flow field is also invariant to the rotation of the y or z axes about the wake centerline. Any directional bias in the wake caused by choice of Re in the uniform density fluid would likely manifest itself in asymmetric vertical and horizontal velocity field distributions between the vertical or horizontal planes about the centerline. The same vertical and horizontal planes used to describe the stratified fluids in Figure 3.7a are also used to describe the near wake for the uniform density fluid. The general form used to parameterize the wake in the vertical plane is of the form: q(x,z) q o (x) Fr = exp " − (z/L σv −B z ) 2 C 2 z + B 2 z C 2 z # +exp " − (z/L σv +B z ) 2 C 2 z + B 2 z C 2 z # (3.36) 73 Fr U/U o u ′ /u ′ o v ′ /v ′ o w ′ /w ′ o ρ ′ /ρ ′ o B z C z B z C z B z C z B z C z B z C z 4 0.00 0.69 0.54 0.50 0.00 1.18 0.00 1.02 0.00 0.95 6 0.00 0.72 0.62 0.54 0.00 1.16 0.00 1.07 0.00 1.24 10 0.00 0.74 0.58 0.52 0.00 1.08 0.00 1.12 0.00 1.43 ∞ 0.00 0.73 0.62 0.60 0.00 0.98 0.00 1.31 n/a n/a Table 3.1: Coefficients for vertical parameterizations in Figure 3.26a,c,d,e,f. n/a: not applicable. Fr U/U o u ′ /u ′ o v ′ /v ′ o w ′ /w ′ o ρ ′ /ρ ′ o B y C y B y C y B y C y B y C y B y C y 4 0.41 0.44 0.63 0.53 0.55 0.47 0.47 0.56 0.00 0.58 6 0.37 0.46 0.63 0.56 0.57 0.62 0.46 0.65 0.00 0.64 10 0.32 0.55 0.61 0.54 0.54 0.76 0.44 0.76 0.00 0.79 ∞ 0.29 0.59 0.61 0.58 0.47 0.91 0.01 1.08 n/a n/a Table 3.2: Coefficients for horizontal parameterizations in Figure 3.26b,d,f,h,j. n/a: not applicable. and, likewise, a description in the horizontal plane is of the form: q(x,y) q o (x) Fr = exp " − (y/L σh −B y ) 2 C 2 y + B 2 y C 2 y # +exp " − (y/L σh +B y ) 2 C 2 y + B 2 y C 2 y # (3.37) where q(x,z) = q(x,y = 0,z), or q(x,y) = q(x,y,z = 0), is a flow field variable, q o represents its centerline value at q(x,y = 0,z =0), and the subscript Fr indicates thespecific internalFroudenumber. AtRe = 1000,thecoefficients B z , C z , B y , and C y aretunableforEq. (3.36) and (3.37)andaredependent onFr andthevariable, q, of interest. By allowing separate descriptions in the vertical and horizontal planes, it is possible to discuss asymmetry in the near wake that exists within the vertical and horizontal planes. In obtaining the coefficients for Eqs. (3.36) and (3.37), the physical coordinates of the horizontal and vertical planes are also normalized by the appropriate half- widthlengthscale; L σv forz andL σh fory. Thevarianceofthedatafortheuniform 74 densityandstratifiednearwakecasesreasonablycollapsesfor2.25≤x.Fr,where the upper limit is equivalent to Nt = 2. Because of the damping layer, Fr = 10 is the highest stratified simulation that reaches Nt = 2 prior to the onset of the damping layer at x = 15. The Fr = ∞ case is only limited by the onset of the damping layer. After the data collapse, the data is numerically fitted to the descriptions of Eq. (3.36) or (3.37), depending on the plane. Following this procedure, Figure 3.26 contains the best-fit curves of the near wake forselected values ofparameter Fr. Fromthe numerically fitted curves of the Fr =∞ case, it appears that assumptions of the axisymmetric wake are satisfied, with only minor asymmetry in the fitted profiles compared between Figure 3.26e,f forv ′ /v ′ o andFigure3.26g,hforw ′ /w ′ o . WithhigherRe,thevariancebetweenfitted profile widths in the uniform density case should decrease. While these arguments are based on best-fit curves, Section 3.2.7 confirms that the Fr =∞ simulation data can be considered axially symmetric. In comparison between the vertical and horizontal planes, the profiles are wake- width normalized, and in the uniform density case, the length scales are identical (Eq. (3.12))suchthatprofilewidthsdirectlycomparetosymmetryofthenearwake width. This is not true for the stratified cases, and althoughprofiles may appear to have similar profile widths in the vertical and horizontal planes, the quantities are physically wider in width than height (Figure 3.7b and Figure 3.16). Additionally, directcomparisonsacrossFrcasesthroughFigure3.26isdifficultbecauseeachflow fieldvariableisnormalizedbyitscenterlinequantity,whichbySection3.2.5isshown to vary with parameter Fr. Quantitative commentary on the physical processes involved in creating the near wake behavior is desirable. This would require an additionalformalanalysisusingeachindividualtermsofEqs. (1.11) and (1.12)and is reserved for future investigations. Qualitative commentary remains possible for 75 someofthedefiningfeaturesofthewakeintheverticalandhorizontalcenterplanes, and the distribution of these quantities at various Fr can be compared. In the vertical centerplane of Figure 3.26a, stratification does not appear to have any pronounced effect on the mean stream-wise velocity distribution within the near wake. In the horizontal centerplane, Figure 3.26b, the mean stream-wise velocity tends towards a wider distribution throughout the wake as Fr decreases. At Fr = 4, the mean stream-wise velocity reaches its maximum value off-centerline which is quite different from the higher Fr cases. In the vertical plane, the stream- wise velocity perturbation quantities, u ′ /u ′ o in Figure 3.26c are higher, relative to their centerline values, than in the uniform density case, but the vertical plane velocity v ′ /v ′ o distributions are unaffected by the presence of stratification (Fig- ure 3.26e). The variance of the stratified case values of w ′ /w ′ o in Figure 3.26g is slightly less than the uniform density case. It appears that once the fluid is strati- fied, thedistributionofw ′ /w ′ o between Fr cases inthenearwake appearsgenerally unaffected, relative to the centerline values, by the severity of stratification. All horizontal plane velocity perturbation quantities (Figure 3.26d,f,h) contain maxi- mum normalized values atoff-centerline locations, and, in relationtotherespective centerline values, the value of this maximum increases as Fr decreases. For the stratified cases, the centerline-normalized perturbation density profiles inFigure3.26i,jdisplaysanaxialasymmetrybetweenverticalandhorizontalplanes whencomparingtheirprofilewidths. Intheverticalplane,asFrincreases, thevari- ance of the profile in the vertical plane increases and appears to create a ”front” of ρ ′ /ρ ′ o . Thisislikelyduetotheeaseofdisplacingaparticleinalightlystratifiedfluid because as Fr increases the ratio of available kinetic energy to the energy required for a vertical fluid displacement also increases. The profile of ρ ′ /ρ ′ o maintains its horizontal distribution within the near wake fairly regularly across the stratified 76 cases, similar to what is seen in w ′ /w ′ o in Figure 3.26h. If ρ ′ is generated primarily through vertical motions, it is not surprising that the horizontal distribution of ρ ′ /ρ ′ o in the wake is relatively unaffected by stratification. The coefficients of Eqs. (3.36) and (3.37) that are obtained from numerically fitting the data to each curve in Figure 3.26 are available in Tables 3.1 and 3.2, respectively. For cases of Fr that are not explicitly listed in Tables 3.1 and 3.2, linear interpolation between entries may be used. If the stratified case of interest requires Fr > 10, substituting a value of Fr ≃ 30 for the Fr = ∞ value is recommended in determining high-Fr interpolated values. For future simulations interested in initialization of the near wake without accounting for the sphere, fully three dimensional profiles can be generated by interpolation between Eqs. (3.36) and (3.37). 3.2.7 Parameterizations vs. Simulation Data at Re=1000 Two examples of the effectiveness of the Re = 1000 near wake parameterizations are given; one for Fr =∞ and the other for Fr = 6. In Figures 3.27 and 3.28 flow field quantities described by Eqs. (3.19), (3.27), (3.29), and (3.33), in addition to (3.35) for the Fr = 6 case, are compared against the time-averaged quantities obtained directly from the simulations. 77 (a) (b) −2 −1 0 1 2 0 0.5 1 U U o −2 −1 0 1 2 0 0.5 1 U U o (c) (d) −2 −1 0 1 2 0 1 2 u ′ u ′ o −2 −1 0 1 2 0 1 2 u ′ u ′ o (e) (f) −2 −1 0 1 2 0 1 2 v ′ v ′ o −2 −1 0 1 2 0 1 2 v ′ v ′ o (g) (h) −2 −1 0 1 2 0 0.5 1 w ′ w ′ o −2 −1 0 1 2 0 0.5 1 w ′ w ′ o (i) (j) −2 −1 0 1 2 0 0.5 1 z/L σ v ρ ′ ρ ′ o −2 −1 0 1 2 0 0.5 1 y/L σ h ρ ′ ρ ′ o Figure3.26: Re=1000. Numericalfitoftime-averageddatafor2.25≤x≤ Fr. Left column is on y = 0 plane. Right column is on z = 0 plane. (△) Fr=4, (+) Fr=6, (♦) Fr=10, (−) Fr=∞. 78 (a) (b) −2 −1 0 1 2 0 0.5 1 U −2 −1 0 1 2 0 0.5 1 U (c) (d) −2 −1 0 1 2 0 0.1 0.2 0.3 u ′ −2 −1 0 1 2 0 0.1 0.2 0.3 u ′ (e) (f) −2 −1 0 1 2 0 0.1 0.2 v ′ −2 −1 0 1 2 0 0.1 0.2 v ′ (g) (h) −2 −1 0 1 2 0 0.1 0.2 z/L σ v w ′ −2 −1 0 1 2 0 0.1 0.2 y/L σ h w ′ Figure 3.27: Re=1000, Fr=∞. Parameterizations versus simulation data. (−) w/ filled symbols are simulation data, (−−) w/ open symbols are parame- terizations. (•) x = 2.25, () x = 3, (N) x = 6, () x =9. 79 By Figure 3.27a-h, at Fr = ∞, the claim of axisymmetry (Section 3.2.6) produced in the simulation is reasonable. Despite some over prediction of v ′ at x ≃ 2.25, the parameterizations appear to replicate the simulation flow field well. Figure 3.28a-j also indicate that the parameterizations of the stratified Fr = 6 case generally reproduce the vertical and horizontal flow fields when {x,Nt} > {2.25,2.25(Fr/2) −1 }. As expected, parameterizations of w ′ /w ′ o in Figure 3.28g,j are not quite as accurate as those for other velocity components, and this is related to the relatively coarse collapse of the centerline quantities seen in Figure 3.22c. The parameterized density perturbation quantities, ρ ′ , in Figure 3.27i,j reproduce the simulation results fairly well. 80 (a) (b) −2 −1 0 1 2 0 0.5 1 U −2 −1 0 1 2 0 0.5 1 U (c) (d) −2 −1 0 1 2 0 0.1 0.2 0.3 u ′ −2 −1 0 1 2 0 0.1 0.2 0.3 u ′ (e) (f) −2 −1 0 1 2 0 0.1 0.2 v ′ −2 −1 0 1 2 0 0.1 0.2 v ′ (g) (h) −2 −1 0 1 2 0 0.1 0.2 w ′ −2 −1 0 1 2 0 0.1 0.2 w ′ (i) (j) −2 −1 0 1 2 0 0.2 z/L σ v ρ ′ −2 −1 0 1 2 0 0.2 y/L σ h ρ ′ Figure 3.28: Re=1000, Fr=6. Parameterizations versus simulation data. (−) w/ filled symbolsaresimulationdata, (−−) w/ open symbols areparameter- izations. (•) {x,Nt} ={2.25,0.75}, () {x,Nt} ={3,1}, (N){x,Nt} ={6,2}, (){x,Nt} ={9,3}. 81 3.3 Wake of the Sphere at Re=5000 Additional investigations are performed using the SA-DES formulation described in Section 2.1.2 at Re = 5000, Fr = 4. The results presented in Section 3.2 at Re = 1000 provide insight on how a stratified near wake simulation should qualita- tively behave without the influence of a turbulent dissipation model. As mentioned inSection2.1.2,theSA-DESapproachhasbeensuccessfullyappliedinuniformden- sity fluids, in particular for Re = 10,000 by Constantinescu & Squires 2003 [14]. This allows an additional investigation into the higher Reynolds number ranges be- tween Re = [10 3 ,10 4 ]. A Reynolds number of 5000 is chosen for exploration of this regime because of available experimental data [53, 51]; also that previous stratified wake simulations [16] have initialized at this condition. Thereisaparticularinterest isinwhetherLESsimulations usingspectralmeth- odscanbeinitializedwithmoreappropriate,physicalinitialconditions. Simulation dataatRe = 5000,Fr = 4iscomparedagainstthestandardinitializationtechnique of Diamessis et al. 2005 [16] and Diamessis et al. 2011 [17]. Their LES domain is initializedbyassumingaspatiallyaveragedx-averagedfieldwhereV xwc =W xwc =0, an axisymmetric Gaussian distribution of perturbation values about the wake cen- terline, and an equipartition of energy amongst the perturbation velocities. This assumption is made forthe velocity perturbation initialcondition by Dommermuth et al. [18], both Diamessis et al. [16, 17], and Bruker & Sarkar [9]. Effectively, this is an azimuthal average of perturbation quantities about the wake centerline in addition to the spatial average of the DPIV arrangement in Section 3.2.2, i.e.: 82 U xwc (r) = 1 2πΔX Z xwc+ ΔX 2 xwc− ΔX 2 dx Z 2π 0 u(x,y,z,t)dθ (3.38) hu ′ i xwc (r) = 1 2πΔX Z xwc+ ΔX 2 xwc− ΔX 2 dx Z 2π 0 [u(x,y,z,t)−U xwc ] 2 dθ !1 2 (3.39) hv ′ i xwc (r) = 1 2πΔX Z xwc+ ΔX 2 xwc− ΔX 2 dx Z 2π 0 [v(x,y,z,t)] 2 dθ !1 2 (3.40) hw ′ i xwc (r) = 1 2πΔX Z xwc+ ΔX 2 xwc− ΔX 2 dx Z 2π 0 [w(x,y,z,t)] 2 dθ !1 2 (3.41) subject to the constraint y 2 +z 2 = r 2 and θ = tan(y/z); where x wc represents the spatial center of the initialized computational domain and r represents the radial distance of the velocity component normal the the wake axis centerline at x wc . Intended as an initial condition, the Gaussian profile does not take into account thetime-variationofthevelocity field. Therefore, theexplicittemporaldependence in Eqs. (3.38)-(3.41) is subsequently ignored. The DES simulation data is collected at the same location that Diamessis et al. [16] initialize their simulations at {x = 6, Nt = 3} The mean stream- wise velocity profile U xwc (r) is compared against the assumed Gaussian profile in Figure 3.29a, and the perturbation quantities, assumed hu ′ i xwc (r) = hv ′ i xwc (r) = hw ′ i xwc (r) =hu irms i ′ , are compared in the same fashion in Figure 3.29b. The DES perturbation velocities hu ′ i xwc (r), hv ′ i xwc (r), hw ′ i xwc (r) in the Figure have been explicitly averaged together from the simulation data. Neither of these profiles accurately describe state of the near wake behind the sphere at Re =5000. The time-averaged centerline mean stream-wise velocity is compared against the experimental data collected by Spedding 2001 [51] in Figure 3.30a. At x = 6, the value of the defect velocity computed by the DES, U o , is comparable to the 83 spatially averaged DPIV data. Outside of the near wake region of Nt≃ 2−3, the resolution of the simulation reduces towards x = 15 and is likely the cause of the discrepancy at that point. Through at least {x = 6,Nt = 3}, the DES is capable of providing the correct centerline defect velocities. As described in Section 3.2.6 for Re = 1000, the mean velocity profile is not symmetric about the wake centerline. The mean flow quantities of the DES simu- lation follow the same trend in the near wake region. By inspection of the vertical andhorizontaltime-averagedvelocity profilesinFigure 3.30b,itisclearthatanax- isymmetric assumption as an initial condition is not correct because the horizontal and vertical length scales of the wake, L σh and L σv , respectively, are not equal. (a) (b) Figure 3.29: Re=5000, Fr=4: Azimuthally averaged at x wc = 6 (a) stream-wise velocity profile (b) hu ′ i i velocity profile at x = 6 (−•−) DES, (−) Diamessis ’05. (i.c.) Inaddition,theequipartitionofenergy, andaxialsymmetry, cannotbeassumed for the perturbation velocities. DES simulation results are compared against the Gaussian fitted profile of Diamessis et al. [16] in Figure 3.29b. In Figure 3.29b, the non-zero asymptotic value with increasing r captured by the averaged DES perturbation quantities, has been fully accounted for in Section 3.2.2, and it is a 84 (a) (b) 10 0 10 1 10 2 10 −2 10 −1 10 0 x U 0 −2 −1 0 1 2 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 y | z U Figure 3.30: Re=5000, Fr=4: (a) U o . (−) DES, (N) Spedding[51]. (b) U x at x =6. ΔT = 36 sampled at Δt = 0.5. (−) vertical plane along y = 0, (−−) horizontal plane along z =0. result of the spatial averaging of lee waves of significant amplitude being generated by the sphere motion through the stratified fluid. Coupled with the Gaussian distribution for the mean centerline stream-wise velocity with a zero asymptote in Figure 3.29a, it immediately follows that the LES initialization technique does not account for the sphere lee waves of significant amplitude. The near wake DES simulation perturbation quantities are shown on the hori- zontal and vertical planes in Figure 3.31a-c. The u ′ o values exhibit an asymmetry on the vertical and horizontal planes, owing to the horizontal and vertical asym- metry of the mean stream-wise velocity profile. The w ′ also captures internal lee waves and matches the data of Spedding [53] reasonably well. It is interesting to note that the DES centerline perturbation quantities for u ′ o , v ′ o , and w ′ o in Figure 3.31 are reasonably equivalent. However, the off-centerline peak magnitudes of w ′ o are not contributing equally when compared to u ′ o and v ′ o . When the DESsimulation perturbation quantities alongthe wake centerline are plotted, as in Figure 3.32, the perturbation quantities are closest at {x = 6,Nt = 85 (a) (b) (c) −2 −1 0 1 2 0 0.02 0.04 0.06 0.08 0.1 y|z u ′ −2 −1 0 1 2 0 0.02 0.04 0.06 0.08 0.1 y v ′ −2 −1 0 1 2 0 0.02 0.04 0.06 0.08 0.1 z w ′ Figure 3.31: Re=5000, Fr=4 at {x = 6,Nt = 3} (a) u ′ (b) v ′ (c) u ′ . ΔT = 36 sampled at Δt = 0.5. (−) vertical plane along y = 0, (−−) horizontal plane along z = 0. 3}. Furtherupstreamfromthisposition,thequantitiesshowsimilartrendsindecay but do not begin with an equipartition of energy behind the recirculation region of the sphere. Downstream of x = 6, w ′ o takes a radical departure from its counter- parts. This correlates with the same Nt = 3 position at Re = 1000 in Figure 3.23 where there is a maximum of potential energy, exhibited by a maximum of density perturbations. At Nt ≃ 4.5, there is a local maximum of w ′ o , corresponding to the point at which there is greater vertical turbulent kinetic energy than potential energy, which correlates to that same position in Figure 3.23. There are two main observations from this comparison. First, the DES simula- tionsarecapableofaccuratelyreproducingthenearwaketurbulentbuoyantbehav- ior and indicates its ability to be used as a predictive tool. Second, the end of the near-wake region, i.e.Nt = 3, is the only specific location at which w ′ o ≃ v ′ o ≃ w ′ o is true, and including upstream and downstream locations under this assumption is inaccurate. 86 10 0 10 1 10 2 10 −3 10 −2 10 −1 10 0 x u i Figure 3.32: Re=5000, Fr=4: DES centerline perturbation velocities. Sampled over ΔT = 36, every Δt = 0.5. (−) u ′ o , (−−) v ′ o , (−·−) w ′ o , (N) Spedding [51] hu ′ i xwc (0,z,t), and (+) Spedding [51] hv ′ i xwc (0,z,t). 87 Chapter 4 Conclusions Simulations havebeenperformedacrossawide parameterspace ofFr atRe = 200, Re = 1000, and Re =5000. At Re = 200, the code was validated against previous experimental and nu- merical investigations. Improvements over the prior numerical investigation [24] at Re = 200 of stratified fluid are shown. Additional analysis is offered over what was available through that study by correctly simulating the shedding of coherent vortices at low Fr. The internal lee wave lengths are measured at multiple heights and found tobe unaffected by the coherent vortex shedding. The current Re = 200 investigation shows that the cause of the discrepancy with experimental results is related to numerical resolution. At Re =1000, the simulations provide quantitative descriptions of the velocity anddensity field inthe nearwake ofthe sphere, andindicate thattheinternal wave field (Figure 3.5) for Fr ≥ 1 is not qualitatively affected by differences between Re = 200 and Re = 1000. The internal lee wave lengths and the separate wake regime behavior for Fr <4 is correlated with the average density field. 88 A detailed investigation explained the inherent differences in results between a DPIVexperimentalsetupandanumericalinvestigationduetothedifferentmethod- ologies in statistical collection techniques. The near wake varies significantly such that time and spatial averaging are not equivalent in the stratified near wake. The ”tails” of the data at low Fr are artifacts of the averaging techniques being applied within the unsteady near wake region. Removing the steady lee wave field velocity components (i.e., w(~ x) in Figure 3.5) from the velocity field before the spatial in- tegration would likely remove the non-zero asymptote, and the temporal variation couldbesmoothedthroughanensembleaveragingofwindoweddatasubsets, aswas done for the simulation results in Figure 3.8b. Internal waves were not detected by inflections ofW(~ x,t)anywhere inthe freestream forFr≥ 5. Regardless, some lee- way should be given in directly comparing between experimentally and numerically collected field data in the near wake of stratified fluids when different averaging techniques are employed. The beginning of the unsteady near wake region starts at x =L rr ≃ 2.25 when Fr ≥ 4. Based on the analysis of L σv , L σh , mean stream wise velocity, and the perturbation quantities for Fr < O(100), there is no region in the wake of the sphere where buoyancy will not affect the near wake at Re = 1000. The near wake region in stratified fluids should be given a better descriptor than the ”3D near wake” lest the statistical behavior of the stratified near wake becomes erroneously synonymous with the behavior of the near wake of a uniform density fluid. A full description of the major field components is described for 4≤Fr≤∞ at a single Re = 1000, and appears to be the first full component field data set available in the stratified near wake. The simulations used DES to run a set of numerical experiments at Re = 5000, Fr = 4 to compare against the DPIV data of Spedding [51] but also comment on 89 the applicability of Gaussian axisymmetric, equipartition initial conditions of wake simulations that do not account for the presence of the sphere explicitly within the domain. The Gaussian approach is explained as an azimuthal and spatial average of the wake. The downstream position Nt = 3 is shown to be a special position at which the distribution of perturbation velocities are close and is correlated to the same behavior seen at the Re =1000 case. The Gaussian distribution matches the DES data when the data is post-processed in that fashion, but the averaging is ultimately misleading as an initial condition because it does not account for the spatial variation or the anisotropic distribution of perturbation velocities in the near wake of the sphere. Additional investigations into the DES simulations show that matches experimental data reasonably well. It is able to reasonably predict the centerline defect velocity, the antisymmetry of the mean velocity profiles in the vertical and the horizontal planes, predicts the perturbation velocities, and most importantly, reproduces the appropriatebuoyant vertical velocity andperturbation density behavior in the near wake. Because of this, DES is expected to perform well as a predictive tool for future stratified wake simulations. Future simulations may attempt to account for the sphere within the computa- tionaldomain intheir own particularway, andaself-contained pointofcomparison is now available for unsteady stratified flows. Other numerical simulations can ini- tialize physically close to the sphere without having to actually account for the sphere explicitly within the computational domain by using the given parameter- izations. The benefits of doing this may result in achieving a higher correlation of simulation-to-physical downstream distances from the sphere, spatially varying initial conditions that account for the density field, and shorter initial simulation relaxation times [9]. 90 Making this claim, it is known that there are two competing requirements in ”proper” initialization of simulations within the near wake that do not explicitly account for the sphere. The first is maintaining a total coherence of the near wake, which would likely require full-interpolations between near wake domains [41]. This is not necessarily a trivial task between various grid geometries espe- cially with respect to the sensitivities to interpolation errors of highly-accurate non-dissipative numerical methods. Second, there is likely a practical requirement toquickly implement the initial conditions into existing numerical schemes without major modification to the previously initialization procedures. The parameterizations offer a relative optimum by spatially varying along the wakecenterlinewhilevaryinginthehorizontalandverticaldirectionstoaccountfor downstream development of the near wake. They can be easily implemented into a numerical simulation while leaving the existing initialization scheme relatively intact, and an implementation guide is in the Appendix for convenient reference. 91 Bibliography [1] A. Arnone, M.S. Liou, and L.A. Povinelli. Integration of Navier-Stokes equa- tionsusingdualtimesteppingandamultigridmethod. AIAA Journal,33:985– 990, 1995. [2] D.J. Bodony. Analysis of sponge zones for computational fluid mechanics. Journal of Computational Physics, 212:681–702, 2006. [3] P. Bonneton, J.M. Chomaz, and E.J. Hopfinger. Internal waves produced by theturbulentwakeofaspheremovinghorizontallyinastratifiedfluid. Journal of Fluid Mechanics, 254:23–40, 1993. [4] P. Bonneton, J.M. Chomaz, E.J. Hopfinger, and M. Perrier. The structure of the turbulent wake and the random internal wave field generated by a moving sphere in a stratified fluid. Dynamics of Atmospheres and Oceans, 23:299–308, 1996. [5] M. Bonnier and O. Eiff. Experimental investigation of the collapse of a turbu- lent wake in a stably stratified fluid. Physics of Fluids, 14(2):791–801, 2002. [6] M. Bonnier, O. Eiff, and P. Bonneton. On the density structure of far-wake vorticesinastratifiedfluid. DynamicsofAtmospheresandOceans,31:117–137, 2000. 92 [7] A. Brandt and C.E. Schemm. Small-scale structure in the near-field of a strat- ified wake. 7th Int. Symp. on Stratified Flows, August 2011. [8] P.W.M. Brighton. Strongly stratified flow past three-dimensional obstacles. Quarterly Journal of the Royal Meteorological Society, 104:289–307, 1978. [9] K.A. Brucker and S. Sarkar. A comparative study of self-propelled and towed wakes in a stratified fluid. Journal of Fluid Mechanics, 652:373–404, 2010. [10] J.M. Chomaz, P. Bonneton, and E.J. Hopfinger. The structure of the near wake of a sphere moving horizontally in a stratified fluid. Journal of Fluid Mechanics, 254:1–21, 1993. [11] G.S. Constantinescu, M. Chapelet, and K. Squires. Turbulence modeling ap- plied to flow over a sphere. AIAA Journal, 41:1733–1742, 2003. [12] G.S.Constantinescu andV.C.Patel. Numericalmodelforsimulationofpump- intake flow and vortices. Journal of Hydraulic Engineering, 124:123–134,1998. [13] G.S. Constantinescu and K. Squires. Numerical investigations of flow over a sphere in the subcritical and supercritical regimes. Physics of Fluids, 16:1449– 1466, 2004. [14] G.S. Constantinescu and K.D. Squires. LES and DES investigations of tur- bulent flow over a sphere at Re = 10,000. Flow, Turbulence and Combustion, 70:267–298, 2003. [15] M.B. de Stadler, S. Sarkar, and K.A. Brucker. Effect of the Prandtl number on a stratified turbulent wake. Physics of Fluids, 22:095102, 2010. 93 [16] P.J. Diamessis, J.A. Domaradzki, and J.S. Hesthaven. A spectral multido- mainpenalty method model forthe simulation ofhigh Reynolds number local- ized incompressible stratified turbulence. Journal of Computational Physics, 202:298–322, 2005. [17] P.J. Diamessis, G.R. Spedding, and J.A. Domaradzki. Similarity scaling and vorticity structure in high-Reynolds-number stably stratified turbulent wakes. Journal of Fluid Mechanics, 671:52–95, 2011. [18] D.G. Dommermuth, J.W. Rottman, G.E. Innis, and E.A. Novikov. Numerical simulation of the wake of a towed sphere in a weakly stratified fluid. Journal of Fluid Mechanics, 473:83–101, 2002. [19] Y.T.FungandS.W.Chang. Surfaceandinternalsignaturesoforganizedvortex motions in stratified fluids. Physics of Fluids, 8:3023, 1996. [20] D.Givoli. High-orderlocalnon-reflectingboundaryconditions: areview. Wave Motion, 39:319–326, 2004. [21] M.J. Gourlay, S.C. Arendt, D.C. Fritts, and J. Werne. Numerical modeling of initially turbulent wakes with net momentum. Physics of Fluids, 13(12):3783– 3802, 2001. [22] M.D. Greeslade. Drag on a sphere moving horizontally in a stratified liquid. Journal of Fluid Mechanics, 418:339–350, 2000. [23] M.P. Gresho and R.L. Sani. On pressure boundary conditions for the in- compressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 7:1111–1145, 1987. 94 [24] H. Hanazaki. A numerical study of three-dimensional stratified flow past a sphere. Journal of Fluid Mechanics, 192(1):393–419, 1988. [25] E.J. Hopfinger, J.-B. Flor, J.-M. Comaz, and P. Bonneton. Internal waves generated by a moving sphere and its wake in a stratified fluid. Experiments in Fluids, 11:255–261, 1991. [26] M. Israeli and S. A. Orszag. Approximation of radiation boundary conditions. Journal of Computational Physics, 41:115–135, 1981. [27] J. Jeong and F. Hussain. On the identification of a vortex. Journal of Fluid Mechanics, 285:69–94, 1995. [28] T.A. Johnson and V.C. Patel. Flow past a sphere up to a Reynolds number of 300. Journal of Fluid Mechanics, 378:19–70, 1999. [29] H.J. Kim and P.A. Durbin. Observations of the frequencies in a sphere wake and of drag increase by acoustic excitation. Physics of Fluids, 31(11):3260– 3265, 1988. [30] M.KokenandG.Constantinescu. Aninvestigationofthedynamicsofcoherent structures in a turbulent channel flow with a vertical sidewall obstruction. Physics of Fluids, 21:085104, 2009. [31] P. Kundu and I. Cohen. Fluid Mechanics, Third Edition. Elsevier Inc., 525 B Street, Suite 1900, San Diego, California, 92101-4495, USA, 2004. [32] S. Lee. A numerical study of the unsteady wake behind a sphere in a uniform flow at moderate Reynolds numbers. Computers & Fluids, 29:639–667, 2000. [33] J. Lighthill. Waves in Fluids. Cambridge University Press, New York, NY, 1978. 95 [34] D.K. Lilly. The representation of small-scale turbulence in numerical simu- lation experiments. In Proc. IBM Sci. Comput. Symp. Environ. Sci., pages 195–210, White Plains, N.Y., 1967. IBM Data Process. Div. [35] Q.Lin,D.L.Boyer, andH.J.S.Fernando. Turbulentwakesoflinearlystratified flow past a sphere. Physics of Fluids A, 4:1687–1696, 1992. [36] K.E. Lofquist and L.P. Purtell. Drag on a sphere moving horizontally through a stratified liquid. Journal of Fluid Mechanics, 148:271–284, 1984. [37] P. Meunier and G.R. Spedding. A loss of memory in stratified momentum wakes. Physics of Fluids, 16(2):298–305, 2004. [38] P. Meunier and G.R. Spedding. Stratified propelled wakes. Journal of Fluid Mechanics, 552:229–256, 2006. [39] R. Mittal. A Fourier-Chebyshev spectral collocation method for simulating flow past spheres and spheroids. International Journal for Numerical Methods in Fluids, 30(7):921–937, 1999. [40] N.V. Nikitin, F. Nicoud, B. Wasistho, K.D. Squires, and P.R. Spalart. An approach to wall modeling in large-eddy simulations. Physics of Fluids, 12(7):1629–1632, July 2000. [41] R.Pasquetti. Temporal/spatialsimulationofthestratifiedfarwakeofasphere. Computers & Fluids, 40:179–187, 2011. [42] F. Roos and W. Willmarth. Some experimental results on sphere and disk drag. AIAA Journal, 9(2):285–291, 1971. 96 [43] J.W. Rottman, K.A. Brucker, D.G. Dommermuth, and D. Broutman. Param- eterization of the Near-Field Internal Wave Field Generated by a Submarine. Symposium On Naval Hydrodynamics, September 2010. [44] H. Sakamoto and H. Haniu. A study on vortex shedding from spheres in a uniform flow. Transactions of the ASME, 112:386–392, 1990. [45] T. Sarpkaya and T. Massidda. Conductivity measurements in the wake of submerged bodies in density-stratified media. In Twenty-First Symposium on Naval Hydrodynamics, pages 266–277. National Academies Press, 1997. [46] M.Shur,P.R.Spalart,M.Strelets,andA.Travin. Detached-eddysimulationof anairfoilathighangleofattack. InProceedingsof4thinternationalsymposium on engineering turbulence modelling and measurements. Elsevier, May 1999. [47] J. Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiment*. Monthly Weather Review, 91(3):99–164, 1963. [48] F.SortiropoulosandS.Abdallah. Thediscretecontinuityequationinprimitive variable solutions of incompressible flow. Journal of Computational Physics, 95:212–227, 1991. [49] P.R. Spalart and S.R. Allmaras. A one-equation model for aerodynamic flows. La Recherch Aerospatiale, 1:1, 1994. [50] P.R. Spalart, W-H. Jou, M. Strelets, and S.R. Allmaras. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. In 1st AFOSR Int. Conf. on DNS/LES. 1st AFOSR Int. Conf. on DNS/LES, Aug. 4-8 1997. 97 [51] G.R. Spedding. Regarding “Anisotropy in turbulence profiles of stratified wakes”. Personal communication. [52] G.R. Spedding. The evolution of initially turbulent bluff-body wakes at high internal Froude number. Journal of Fluid Mechanics, 337:283–301, 1997. [53] G.R. Spedding. Anisotropy in turbulence profiles of stratified wakes. Physics of Fluids, 13(8):2361–2372, Aug. 2001. [54] G.R. Spedding. Vertical structure in stratified wakes with high initial Froude number. Journal of Fluid Mechanics, 454:71–112, 2002. [55] G.R. Spedding, F.K. Browand, and A.M. Fincham. Turbulence, similarity scalingandvortexgeometryinthewake ofatowedsphere inastablystratified fluid. Journal of Fluid Mechanics, 314:53–103, 1996. [56] P.J. Stockhausen, C.B. Clark, and J.F. Kennedy. Three-dimensional momen- tumless wakes in density-stratified liquids. Technical Report 93, Department of Civil Engineering, Massachusetts Institute of Technology: Hydrodynamics Laboratory, June 1966. [57] M. Strelets. Detached eddy simulation of massively separated flows. In 39th AIAA Aerospace Sciences Meeting and Exhibit.39thAIAAAerospace Sciences Meeting and Exhibit, Jan. 8-11 2001. [58] S. Taneda. Experimental investigation of the wake behind a sphere at low Reynolds numbers. Journal of the Physical Society of Japan, 11:1104–1108, 1956. 98 [59] A. G. Tomboulides and S. A. Orszag. Numerical investigation of transitional andweak turbulent flowpastasphere. Journal of Fluid Mechanics,416:45–73, 2000. [60] J.S. Turner. Buoyancy Effects in Fluids. Cambridge University Press, New York, NY, 1973. [61] B. Voisin. Lee waves from a sphere in a stratified flow. Journal of Fluid Mechanics, 574:273–315, 2007. [62] D. Wilcox. Turbulence Modeling for CFD. DCW Industries, Inc., 1998. [63] J.-S. Wu and G.M. Faeth. Sphere wakes in still surroundings at intermediate Reynolds numbers. AIAA Journal, 31(8):1448–1455, 1993. 99 Appendix A Procedure to Initialize Sphere-less Simulations The procedure as well as relevant equations are reproduced here for convenience. At each node: Step 1. Calculate L σv and L σh by: L σvFr = A(Fr)x 0.55 0.94+3.2/(Fr−4) 2.25. x. 3 2 Fr (A.1) L σhFr = B(Fr)x 0.55 2.25. x. 3 2 Fr (A.2) A(Fr) = 0.11sech[0.24(Fr−7.42)]+0.37 B(Fr) = 0.46exp(−0.31Fr)+0.37 100 Step 2. Use values of L σv and L σh from Step 1 to calculate U o , u ′ o , v ′ o , w ′ o , and ρ ′ o by: U oFr L σvFr =9.68x −2.19 2.25. x.Fr (A.3) u ′ oFr L σvFr =0.87x −1.49 2.25. x. 3 2 Fr (A.4) v ′ oFr L σvFr =3.95f v (Fr)x −1.43 2.25. x. 3 2 Fr (A.5) w ′ oFr L σvFr =2.9f w (Fr,x)x −1.49 2.25≤ x≤ 3 2 Fr (A.6) ρ ′ oFr =0.056f ρ ′(Fr)x 17.5Fr −1 2.25≤ x. 3 2 Fr (A.7) f v (Fr) = 0.26tanh 0.39 Fr 2 f w (Fr,x) = 0.15 " 2+sin 1.16 π 2 2x Fr 0.94 !# f ρ ′(Fr) = tanh Fr 6 exp −15Fr −2 Step 3. Use centerline values from Step 2 to calculate U(x,y,z), u ′ (x,y,z), v ′ (x,y,z),w ′ (x,y,z),andρ ′ (x,y,z)byusing(orinterpolating)thecoefficients from Tables 3.1 and 3.2 along with the form of Eqs. (3.36) and (3.37) to complete the spatially varying flow field: q(x,z) q o (x) Fr = exp " − (z/L σv −B z ) 2 C 2 z + B 2 z C 2 z # +exp " − (z/L σv +B z ) 2 C 2 z + B 2 z C 2 z # q(x,y) q o (x) Fr = exp " − (y/L σh −B y ) 2 C 2 y + B 2 y C 2 y # +exp " − (y/L σh +B y ) 2 C 2 y + B 2 y C 2 y # Step 4. Complete any remaining standard simulation initialization procedures. 101
Abstract (if available)
Abstract
A numerical investigation of the near wake of a sphere moving horizontally through a linearly stratified fluid. Simulations are performed at a low Reynolds number of 200 for a range of internal Froude number. Simulations at Reynolds number of 1000 across a range of Froude number greater than 1 provide a description and parameterization of the near wake, including the density field. The effects of utilizing two different averaging techniques in the unsteady near wake region are discussed. Perturbation quantities in the stratified near wake are anisotropic, and the Froude number at which the stratified flow can be treated as a uniform density fluid is suggested to be on the order of 100. Parameterization of the near wake is accomplished using the parameterized wake height, downstream distance from sphere, and Froude number as parameters.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Wake modes of rotationally oscillating circular cylinder in cross-flow and its relationship with heat transfer
PDF
Pattern generation in stratified wakes
PDF
Large eddy simulations of laminar separation bubble flows
PDF
On the simulation of stratified turbulent flows
PDF
Near wake characteristics of towed bodies in a stably stratified fluid
PDF
Passive flight in density-stratified fluid environments
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
A Boltzmann model for tracking aerosols in turbulent flows
PDF
Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction
PDF
Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders
PDF
Reconstruction and estimation of shear flows using physics-based and data-driven models
PDF
Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles
PDF
Tracking and evolution of compressible turbulent flow structures
PDF
Boundary layer and separation control on wings at low Reynolds numbers
PDF
Model based design of porous and patterned surfaces for passive turbulence control
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Flow and thermal transport at porous interfaces
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
Asset Metadata
Creator
Orr, Trevor Stuart
(author)
Core Title
Numerical simulations of linearly stratified flow around a sphere
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Aerospace and Mechanical Engineering (Computational Fluid and Solid Mechanics)
Publication Date
09/24/2014
Defense Date
09/02/2014
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
anisotropic wake,density field,density perturbation,detached eddy simulation,direct numerical simulation,direct numerical simulation of stratified flow,DPIV,Froude number,internal wave,linearly stratified,linearly stratified flow,near wake,NEQ regime,NEQ region,non-equilibrium regime,non-equilibrium region,OAI-PMH Harvest,pancake vortices,spatial averaging,Sphere,sphere wake,sphere wakes,sphere-less initialization,stratified flow,stratified near wake,stratified wakes,thermally stratified,vortex shedding,Wake,wake parameterization,wake parametrization
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Domaradzki, Julian A. (
committee chair
), Eliasson, Veronica (
committee member
), Nakano, Aiichiro (
committee member
), Redekopp, Larry G. (
committee member
), Spedding, Geoffrey R. (
committee member
)
Creator Email
trevoror@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-484651
Unique identifier
UC11287214
Identifier
etd-OrrTrevorS-2988.pdf (filename),usctheses-c3-484651 (legacy record id)
Legacy Identifier
etd-OrrTrevorS-2988.pdf
Dmrecord
484651
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Orr, Trevor Stuart
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
anisotropic wake
density field
density perturbation
detached eddy simulation
direct numerical simulation
direct numerical simulation of stratified flow
DPIV
Froude number
internal wave
linearly stratified
linearly stratified flow
near wake
NEQ regime
NEQ region
non-equilibrium regime
non-equilibrium region
pancake vortices
spatial averaging
sphere wake
sphere wakes
sphere-less initialization
stratified flow
stratified near wake
stratified wakes
thermally stratified
vortex shedding
wake parameterization
wake parametrization