University of Southern California Dissertations and Theses

On the simulation of stratified turbulent flows

(USC Thesis Other)

On the simulation of stratified turbulent flows

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content ON THE SIMULATION OF STRATIFIED TURBULENT FLOWS by Yuncheng Lin 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 2010 Copyright 2010 Yuncheng Lin Dedication This work is dedicated to my wife Cindy Yuanchun Huang. ii Acknowledgements I want to thank my wife for her endless support and love during the time I was weak, poor and grumpy. It would take infinite amount of money, although not enough, to pay back the debt I owed. I want to thank Professor Peter Diamessis forhisencouragementandguidanceasIwasconfusingandstrugglingwiththenew environmentsatUSC.IalsowanttothankProfessorJulianDomaradzkiforputting up with me who naively and crazily wanted to learn everything without knowing anything. I also want to thank for the door he opened for me to explore and enjoy the beauty and the mistery of turbulence, a subject I still eager to know more. Finally, I want to thank Professor Redekopp for his endless patience and under- standing in directing a student who is constantly struggling with reckless mistakes and finishing jobs on time. I also want to thank all my friends in the AME depart- ment for the friendship and support. I want to also thank my families and friends on the other side of the earth for the support and patience. Finally, I want to thank God for making it possible to love again. iii Table of Contents Dedication ii Acknowledgements iii List Of Tables vi List Of Figures vii Abstract xi Part I − Effective Numerical Viscosity on SMPM 1 Part I, Chapter 1: Literature Review and Introduction 1 Part I, Chapter 2: Problem Statement 7 2.1 DNS Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Numerical Stabilization Techniques . . . . . . . . . . . . . . . . . . 12 2.2.1 Penalty Method . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Spectral Filtering . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Effective Numerical Viscosity . . . . . . . . . . . . . . . . . . . . . 17 Part I, Chapter 3: Results 23 3.1 Horizontal Center Plane . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Model Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Interfacial Subdomains . . . . . . . . . . . . . . . . . . . . . . . . . 27 Part I, Chapter 4: Conclusions 35 Part II − The Wave-Induced Boundary Layer Under In- ternal Solitary Waves 39 Part II, Chapter 1: Introduction 39 1.1 Internal Solitary Waves . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.2 The Wave-Induced Boundary Layer . . . . . . . . . . . . . . . . . . 43 iv Part II, Chapter 2: Problem Definition 49 2.1 Base Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.1.1 Solitary Wave Profile . . . . . . . . . . . . . . . . . . . . . . 50 2.1.2 Solitary Wave Characteristics . . . . . . . . . . . . . . . . . 55 2.1.3 Inviscid Flow Field . . . . . . . . . . . . . . . . . . . . . . . 59 2.1.4 Boundary layer Basics . . . . . . . . . . . . . . . . . . . . . 63 2.2 Turbulent Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 2.2.1 Reynolds-averaged equation . . . . . . . . . . . . . . . . . . 71 2.2.2 k−ω model . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Part II, Chapter 3: Numerical Issues 78 3.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.2 Grid and Accuracy Tests . . . . . . . . . . . . . . . . . . . . . . . . 81 3.3 Velocity Profiles and Eddy Viscosity . . . . . . . . . . . . . . . . . 83 3.3.1 Sample Velocity Profiles . . . . . . . . . . . . . . . . . . . . 84 3.3.2 Eddy Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . 85 Part II, Chapter 4: Boundary layer Properties 91 4.1 Skin friction coefficient . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.2 Boundary Layer Thickness . . . . . . . . . . . . . . . . . . . . . . . 97 4.3 Separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Part II, Chapter 5: Concluding Perspectives 110 References 112 Appendix A Spectral Energy Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 116 Appendix B B.1 Wave Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 B.2 Dimensionless Pressure Gradient . . . . . . . . . . . . . . . . . . . 120 B.3 Compressible To Incompressible Corrections . . . . . . . . . . . . . 120 v List Of Tables 3.1 Filter Configurations for Re=5000 . . . . . . . . . . . . . . . . . . 23 2.1 ζ o : dimensionless extreme wave amplitude . . . . . . . . . . . . . . 57 3.1 GCI Report for C f . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.2 GCI Report for δ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.3 GCI Report for δ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.1 Shape factor H =δ 1 /δ 2 for typical boundary layer. . . . . . . . . . 99 vi List Of Figures 2.1 A sphere of diameter D is towed with a speed U in the x direction. 8 2.2 Spectral filters used in Fourier space. See equation (2.17).. . . . . . 17 2.3 The decay of total kinetic energy over the course of 8 timesteps with minmum filtering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1 Effective numerical viscosity normalized by the molecular viscosity, i.e., ν/ν o , atthehorizontalcenterplane. k c representsthemaximum wavenumber, namely, k c = p k 2 max,x +k 2 max,y . The energy spectrum E(k,z) used in equation (2.29) and (2.30) is corresponding to the baseline case, i.e., minimum dissipation. . . . . . . . . . . . . . . . 30 3.2 Refer to the caption in figure (3.1). The energy spectrum E(k,z) usedinequation(2.29)and(2.30)iscorrespondingtothedissipative case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Thecomparisonbetweentheresultscalculatedfromequations(3.3a), (i.e., the enhanced lines), and (3.3b). . . . . . . . . . . . . . . . . . 31 3.4 The comparison of Effective numerical viscosity between the model simulation(withoutNS)andtheactualsimulationwithonlyFourier filter turned on. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.5 Effective numerical viscosity, ν/ν o , calculated at bottom and top subdomain interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.6 The effective numerical viscosities at the neighboring points of the bottom subdomain interface. . . . . . . . . . . . . . . . . . . . . . . 34 vii 2.1 A definition sketch of a stationary solitary wave and its boundary conditions. The depth ratio, h 2 /h 1 = h, is greater than 1 for a de- pression wave. For stable stratification, ρ 2 > ρ 1 . The wave-induced velocity fields in the streamwise and transverse direction are u(x ∗ ) and v(x ∗ ,y ∗ ) respectively. u, v, x ∗ , y ∗ , η, h 1 , h 2 , ρ 1 and ρ 2 are all dimensional variables.. . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.2 Solitary wave profiles for depth ratio h = 8. x c is the center of the wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.3 The variation of the extreme wave amplitude and the corresponding wave speed with respect to the depth ratio. . . . . . . . . . . . . . . 58 2.4 The variation of characteristic wave length with respect to the wave amplitude ζ o (left panel) as well as the wave amplitude ratio ζ o /ζ oe (right panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.5 (a).ζ o /ζ oe ∼ 0.5. The legend shows the depth ratio h = h 2 /h 1 . (b). h=8. ; The legend shows the amplitude ratio ζ o /ζ oe . . . . . . . . . 62 2.6 Refer to equation (2.20). h=8, ζ o =0.5,2.0,3.3. . . . . . . . . . . 64 2.7 The velocity difference under the peak of the wave. See equation (2.21). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.8 Re w,c :the wave-center Reynolds number, defined in equation (2.22). 67 2.9 The variation of F w andF wc with respect to the amplitude ratio. F w and F wc are defined in equation (2.22) and (2.24). . . . . . . . . . . 67 2.10 Theareaofcomputationisconfinedbythedashlinesandthebottom wall, i.e., y ∗ c ·L ∗ c k 1 =k 2 =6,y ∗ c >>δ(x ∗ ). The characteristic wave length L ∗ w is computed from equation (2.13), and is shown here as merely an visual estimate. . . . . . . . . . . . . . . . . . . . . . . . 77 3.1 Iterative convergence of inlet flow (see equation (3.1)). h = 8, ζ o = 3.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.2 (a):◦ ; (b):Δ ; (c):¦ ; (d):¤ (See figure (2.2)). . . . . . . . . . . 85 3.3 (a):◦ ; (b):Δ ; (c):¦ ; (d):¤ (See figure (2.2)). . . . . . . . . . . 86 viii 3.4 The variation of maximum viscosity ratio with respect to ζ o . Re w = 10 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.5 transversevariationofviscosityratioatwavecenter. h=8;ζ o /ζ oe = 2.0/3.42,2.5/3.42,3.0/3.42 and 3.4/3.42. Re w =10 8 . . . . . . . . . . 88 3.6 Thestreamwisevariationofviscosityratiocorrespondingtothewaves in figure (2.2). h=8 ; Re w =10 8 . . . . . . . . . . . . . . . . . . . . 89 4.1 top: wave profiles ; bottom: the corresponding C f profiles. Re w = 10 8 ; numbers in the legends are the amplitude ratio, ζ o /ζ oe . . . . . . 103 4.2 The variation of the maximum and total skin friction coefficients with respect to the wave Reynolds number. . . . . . . . . . . . . . . 104 4.3 Replotofskinfrictionvariationinlog-logscale: −:ζ o /ζ oe =0.5/3.42; −.−:ζ o /ζ oe =2.0/3.42.. . . . . . . . . . . . . . . . . . . . . . . . . 104 4.4 The variation of skin friction coefficient against local Reynolds num- ber, Re δ 2 ,U c . h=8. ζ oe =3.42. . . . . . . . . . . . . . . . . . . . . . 105 4.5 The variation of maximum and total skin friction coefficients with respect to amplitude ratio. Re w =10 8 . . . . . . . . . . . . . . . . . 105 4.6 ThevariationofskinfrictioncoefficientswithrespecttolocalReynolds number for three different wall roughnesses. δ ∗ 2 . : momentum thick- ness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.7 The variation of displacement and momentum thicknesses at the center of the wave with respect to the amplitude ratio. δ 1 and δ 2 are scaled by h 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.8 Thevariationofboundary-layerthicknesseswithrespecttothewave Reynolds number Re w . h = 8. δ 1 and δ 2 are scaled by the upper layer depth h 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.9 The variation of shape factor at the wave center with respect to the wave amplitude ζ o and the amplitude ratio ζ o /ζ oe . Re w =10 8 . . . . 107 4.10 The variation of shape factor as a function of amplitude ratio at the leading inflection point (L.I.), wave center (C.) and rear inflection point (R.I.). Re w =10 8 . h=8. . . . . . . . . . . . . . . . . . . . . 108 ix 4.11 The variation of D SR (defined in equation (4.7)) with respect to the amplitude ratio. Re w =10 8 . . . . . . . . . . . . . . . . . . . . . . . 108 4.12 D SR :(x sep −x R.I. )/L w . Re w =10 8 . . . . . . . . . . . . . . . . . . . 109 x Abstract In the first part of this report, the effects of numerical dissipation presented in a turbulence Direct Numerical Simulation (DNS) code called SMPM is investigated through the use of the effective numerical eddy viscosity, a concept pioneered by Domaradzki and Xiao in 2003. This SMPM algorithm, which is inherently free of any intrinsic truncation error, employs three numerical stabilizers, i.e., Fourier filtering,LegendrefilteringandPenaltymethodtoensurethenumericalstabilityin the calculations of turbulent flows with very high Reynolds number, for which the simulationisoftenunder-resolved. Numericaldataextractedfromthesimulationof aturbulentwakeflowisemployedtoquantifytheeffectivenumericaleddyviscosity for all three explicit numerical stabilizers. The results shown that the effects of the stabilizers can be significant, if not dominated, compared to the physical molecular viscosity. Away from the vertical subdomain interfaces, all three stabilizers behave expectedly to remove the kinetic energy from the system, i.e., dissipative. At the subdomain interfaces, the Legendre filtering and the Penalty scheme exhibit the unexpected anti-dissipative character, i.e., to increase the energy of the system. Such counter-intuitive behavior is attributed to the discontinous formulism of the xi Penaltyscheme, whichisabsolutelynecessaryifmulti-domainmethodisemployed. This anti-dissipative behavior diminishes quickly only a few grid points away from the subdomain interfaces. Theoriginalformulismofeffectivenumericalviscositywasdedicatedtoaturbu- lence DNS code with triply-periodic boundary conditions. It has been successfully extended to our SMPM code whose boundary condition is not periodic in the ver- tical direction, the direction aligns with the only non-zero external force field, i.e., gravity. Thesuccessofthisextensiongivesustheconfidencethatthesameconcept can be applied to any numerical algorithm for under-resolved simulation, so long as the baseline case can be defined. xii The second part of this report is dedicated to discovering the nature of bound- ary layer flows induced by the propagation of the internal solitary waves. Long internal waves are common features on the continental shelf and in lakes, but their dissipation via benthic boundary layer drag is largely unknown, particularly when the wave amplitudes are large and the boundary layer corrections based on the linear theory are clearly invalid. In general, the boundary layer induced by a soli- tary wave of depression experiences a continuous favorable-to-adverse variation of the pressure gradient, undergoes transition, may reach a strongly turbulent state, and frequently separates near the point of maximum adverse pressure gradient in the lee of the wave. In this study a model for fully-nonlinear solitary waves of depression in a two-layer stratification is employed as the inviscid base state, and a RANS solver with k− ω turbulence model is used to compute the stationary boundary layer under the wave. Local friction coefficients and eddy viscosities are computed in the footprint of the wave. Locations of boundary layer separation are computed as well as the integrated frictional drag over the region of attached boundary layer flow. Boundary layer characteristics are presented for a range of environmental conditions, Reynolds numbers, and surface roughness in an attempt to provide a quantitative measure of the frictional drag of long internal waves in realistic, shallow environments. xiii Part I - Effective Numerical Viscosity on SMPM Part I, Chapter 1 Literature Review and Introduction Incompressible fluid flow is described by Navier-Stokes equations: ∂u i ∂x i =0, (1.1a) ∂u i ∂t +u j ∂u i ∂x j =− 1 ρ ∂p ∂x i +ν ∂ 2 u i ∂x j ∂x j +f i , (1.1b) whereu i isvelocitycomponentintheidirection, pisthepressure, ρisthedensity, ν is the kinematic viscosity and f i is the force per unit mass in i direction, (i.e. i = 1,2,3). The only external force being considered in this study is the uniform gravitational force, i.e., f i =(0,0,g), (1.2) 1 whereg istheconstantgravitationalacceleration. Theaboveequations(1.1)canbe re-castedintoanon-dimensionalformwiththeproperchoiceofcharacteristiclength and time scales. The resulting equations are mainly dependent on two parameters, namely, the Reynolds number Re and the Froude number Fr, defined as Re= UL ν , Fr = U √ gL , (1.3) where U and L are the characteristic velocity and length scales of the flow system. The dynamics of the fluid flow varies dramatically with the variation of these pa- rameters. Most notably, turbulence can arise as a consequence of the structure of Navier-Stokes equations when the Re is sufficiently high. It is then overwhelmingly expensive to resolve all the active modes computationally for the large-Re flow and as a result, a high-Re flow is often under-resolved. In other words, the dynamics of the molecular viscosity, which corresponds to the smallest physical scale, is not fully captured numericallly. A numerically under-resolved simulation is not only physically inaccurate, but also it might render itself to be numerically unstable to theextentthatnousefulresultcanbeobtainedatall. Inordertoavoidthisdespair situation, two mechanisms are invented to compensate the unresolved dynamics. Oneistouseamodeltorepresenttheunresolveddynamicsbasedontheknowledge of the physics of turbulence. The analysis of averaged-Navier-Stokes equation, i.e., RANS,withtheconceptofeddyviscosityisoneofthemodelscommonlyusedwith high-Re turbulent flows. Over the past two decades, the Large Eddy Simulation 2 (LES) and Detached-Eddy Simulation (DES) have been the focus of turbulence modeling community and are the promising candidates to capture the accurate physics for a large number of practical applications of high-Re flows. Some reviews are in [Pop00], [Les96], [Pio99] and [Spa09]. The other method to deal with the issue of under-resolution is to introduce certain kind of numerical dissipation to account for the unresolved molecular dissipation. The truncation error obtained by directly discretizing Navier-Stokes equations serves as an example of implicit numerical dissipation. This implicit nature poses a challenge for the researchers to answer the question whether or not the unresolved physics is properly repre- sented by the implicit numerical dissipation. To resolve this issue, Domaradzki and Xiao [DX03] pioneered a concept called ”effective numerical eddy viscosity” in an attempt to quantify the implicit numerical dissipation inevitably presented in a Direct-Numerical-Simulation (DNS) solver called ”monotonically integrated large eddysimulations”orMILES.Thisconceptisformulatedbyusingaspectralkinetic energy equation for isotropic turbulence, ∂ ∂t E(k,t)=−2νk 2 E(k,t)+T(k,t), (1.4) where E(k,t)Δk is the total spectral kinetic energy corresponding to a band of wavenumbers defined by k−Δk/2 ≤ k ≤ k +Δk/2, and k = (k 2 x +k 2 y +k 2 z ) 1/2 , the magnitude of the wavenumber. E(k,t) then is regarded as the energy density function. T(k,t) describes the nonlinear energy transfer through the interactions 3 between modes of different wavenumbers (eddies of different sizes). Notice that the integration of T(k,t) over all possible k should be zero for all timesteps, that is, Z ∞ 0 T(k,t)dk =0. (1.5) The mathematical derivation of equation (1.4) utilizes the Fourier transform of the Navier-Stokes equations and is outlined in Appendix (A). With the evolution equation for the spectral energy kinetic (equation (1.4)) at hand, the effective viscosity is numerically computed as follows. First, the velocity components in the spectral space ˜ u(k) are obtained by the application of discrete Fouriertransformuponthevelocityfieldscalculatedfromthesolverinspatialspace, i.e., u(x). Secondly, the spectral kinetic energy E(k) is computed by E(k)= 1 2 ˜ u i (k)˜ u ∗ i (k), (1.6) where ˜ u ∗ i (k) is the complex conjugate of ˜ u i (k). Finally, three timesteps of E(k) are saved to calculate the time rate of change of E(k), i.e., dE(k)/dt, using the third-order backward differential scheme, that is, µ dE dt ¶ N = E(k,t−2Δt)−4E(k,t−Δt)+3E(k,t) 2Δt . (1.7) 4 The subscript N stands for the numerical calculation. Owing to the existence of the implicit numerical dissipation originated from the truncation error, terms in equation (1.4) will not balance one another numerically and a residual is defined as the effective numerical dissipation, ² n (k,t)= µ dE dt ¶ N +2νk 2 E(k,t)−T(k,t). (1.8) Notice that ² n (k,t) ≈ 0 if the truncation error is negligible, i.e., fully-resolved simulation. As a result, a k-dependent numerical viscosity can be defined through the analogy with the dissipation due to molecular viscosity, the first term on the right hand side of equation (1.4), that is, ν n (k)= ² n (k) 2k 2 E(k) . (1.9) This quantity is then compared with the molecular viscosity used in actual simula- tions, ν, in part to estimate the performance and accuracy of the numerical solver. The equation (1.9) is valid with the use of Fourier analysis, i.e., triply-periodic boundary condition. It might have a different expression if different numerical schemewithdifferentboundaryconfigurationistobeanalyzed. Forinstance, when ahigh-orderdiscretizationschemeisusedwithorthogonalpolynomials,i.e.,Cheby- shev or Legendre, the artificial dissipation originated from the truncation error is 5 negligiblewithuseofadequatepolynomialorder. Tobemorespecific, aDirectNu- mericalSimulation(DNS)withuse ofthe abovehigh-order discretizationschemeis well-resolved without any implicit numerical dissipation, i.e., ² n (k,t)≈0, provided that the Reynolds number is not too high. However, when the simulation is under- resolved, i.e., high-Re flow, there is no intrinsic or implicit numerical mechanism to make up the uncaptured dynamics associated with the molecular viscosity. There- fore, one has to introduce either an explicit numerical dissipative mechanism, i.e., spectral filtering, or penalty multidomain collocation scheme, if an element-based discretizationisemployed,toobtainthecorrectphysicsofthesmall-scaledynamics. The aim of this study is to reveal the role played by these explicit dissipative mechanisms employed in a DNS solver which uses a spectral multidomain penalty method model to simulate high-Re turbulent flows [DDH05]. The remaining chap- ters are arranged as follows. The associated DNS solver, the stabilization tech- niques as well as the formulation of the associated effective eddy viscosity will be discussed in Chapter 2, followed by the results and the conclusions in Chapter 3 and 4 respectively. 6 Part I, Chapter 2 Problem Statement Thenumericalsimulationunderinvestigationisthemodelingoftheturbulentwake generated by towing a sphere of diameter D at a characteristic velocity U. The physical geometry is described in Figure (2.1). The numerical solver used to simu- late the wake is introduced first, followed by the detailed description of the explicit numerical stabilizers whose performance is to be analyzed. The last section of this chapter is dedicated to the formulation of the spectral energy equation based on which the effecitve numerical viscosities are defined. 2.1 DNS Solver TheDNSsolverthatisbeinganalyzedwasdevelopedbyDiamessisetal. [DDH05]. A general review is given in this section. This solver uses a spectral multidomain penalty method (SMPM) to model incompressible stratified turbulent flows with 7 Figure 2.1: A sphere of diameter D is towed with a speed U in the x direction. highReynoldsnumber. TheincompressiblestratifiedflowisgovernedbytheNavier- Stokes equations with Boussinesq approximation [KC08], ∂u ∂t =− 1 2 {u·∇u+∇(u·u)}+F g − 1 ρ o ∇p 0 +ν∇ 2 u. (2.1a) ∂ρ 0 ∂t =−∇·[u(ρ 0 +ρ(z))]+κ∇ 2 ρ 0 . (2.1b) ∇·u=0. (2.1c) F g =−g ρ 0 ρ o ˆ e k , (2.1d) where g is the gravitational constant, ν is the kinematic viscosity and κ is the thermal diffusivity. The nonlinear term, i.e., the second term on the left hand side of the equation (1.1b), is rearranged as a skew-symmetric form, namely, the first two terms on the right hand side of equation (2.1a), to improve numerical stability 8 [DFM02]. The Boussinesq Approximation governs the relation between the mean quantities and perturbations of density and pressure, p=p(x,y,z)+p 0 (x,y,z,t), (2.2a) ρ=ρ o +ρ(z)+ρ 0 (x,y,z,t), (2.2b) ∂p ∂z =−(ρ o +ρ)g. (2.2c) Equation (2.2c) represents the hydrostatic balance between the mean pressure, p, and the reference density, ρ o +ρ. The proposed Navier-Stokes equations (2.1) are solved with periodic boundary conditions in both x and y directions. In the vertical direction, the top boundary is a shear-stress free surface and a rigid wall at the bottom. That is, (u;v;w;p 0 ;ρ 0 )(x,y,z,t)=(u;v;w;p 0 ;ρ 0 )(x+L x ,y,z,t). (2.3a) (u;v;w;p 0 ;ρ 0 )(x,y,z,t)=(u;v;w;p 0 ;ρ 0 )(x,y+L y ,z,t). (2.3b) u(x,y,0,t)=0, v(x,y,0,t)=0, w(x,y,0,t)=0, ρ 0 (x,y,0,t)=0. (2.3c) ∂u ∂z | (x,y,Lz,t) =0, ∂v ∂z | (x,y,Lz,t) =0, w(x,y,L z ,t)=0, ρ 0 (x,y,L z ,t)=0. (2.3d) The numerical solution is proceeded in time using pressure projection scheme developed by Karniadakis et al. [KIO91]. This particular scheme splits the time 9 interval, Δt = t n+1 −t n , into three fractional timesteps, each of which deals with nonlinear interaction, pressure correction and viscous dissipation respectively. ˆ u− P J i −1 q=0 α q u n−q Δt = J e −1 X q=0 β q N(u n−q ). (2.4a) ˆ ˆ u− ˆ u Δt =−∇φ n+1 . (2.4b) γ o u n+1 − ˆ ˆ u Δt =νL(u n+1 )=∇ 2 u n+1 , (2.4c) where ˆ u and ˆ ˆ u are the intermediate velocities between fractional timesteps. The equation (2.4a) is the first fractional timestep and is solved using explicitly numer- ical scheme. N(u) is the combination of the skew-symmetric nonlinear term and the buoyance force, i.e., N(u)=− 1 2 {u·∇u+∇(u·u)}+F g . (2.5) The coefficients, J i ,J e ,α q ,β q are corresponding to a third-order backward differen- tialschemeandcanbefoundin[DFM02],or[KIO91]. Theintermediatescalarfunc- tion φ n+1 is there to ensure the incompressiblity of the solutions of full timesteps, that is,∇·u n+1 =0. The scalar field is related to the pressure through Z t n+1 t n ∇p 0 dt=Δt∇φ n+1 . (2.6) 10 By taking the divergence of the equation (2.4b), one obtains, ∇· ˆ ˆ u−∇·ˆ u Δt =−∇·∇φ n+1 . (2.7) With the assumption that ∇· ˆ ˆ u = 0, the above Poisson equation is solved with boundary condition proposed by Karniadakis et al. [KIO91], ∂φ n+1 ∂z ¯ ¯ ¯ ¯ z=0,Lz = J e −1 X q=0 β q N(w n−q ) ¯ ¯ ¯ ¯ ¯ z=0,L z − J e −1 X q=0 β q {ν∇×(∇×w)} n−q ¯ ¯ ¯ ¯ ¯ z=0,L z . (2.8) to ensure high numerical accuracy for the second fractional timestep φ n+1 , which is then used in the final viscous step (2.4c) to obtainu n+1 with fully implicit scheme. The detailed derivation of equation (2.8) is available in [KIO91]. The initial condi- tion of the simulation is obtained by using a noise generator generating velocities corresponging to a turbulent wake described in Figure (2.1). Velocity components of 8 consecutive timesteps are saved to calculate the energy decay rate as well as the effective numerical viscosity. FourierspectraldiscretizationisusedwithN x andN y Fouriermodesinperiodic horizontaldirections. Intheverticaldirection,thecomputationaldomainisdecom- posed into M number of subdomains and the Legendre spectral discretization is 11 employed within each subdomain. This solver, SMPM, is designed as a collocation method in which any function g(z) is approximated in nodal form, i.e., g(z)= N k −1 X j=0 g(z j )l j (z), (2.9) where l j (z) is the j th -order Cardinal function of Legendre kind with the use of N k Legendre-Gauss-Lobatto collocation points [Boy01]. See [DDH05] for more infor- mation about this DNS solver. 2.2 Numerical Stabilization Techniques Owing to the numerical nature of under-resolution of high-Re flow, explicit in- troduction of numerical dissipation is necessary to stabilize the simulation. Two different numerical dissipative mechanisms, i.e., Penalty scheme and spectral filter- ing, are employed at different stages of the simulation. 2.2.1 Penalty Method ThemainconceptofPenaltymethodorPenaltyschemecanbefoundin[Hes97]. It is applied at the interfaces of subdomains at both nonlinear and viscous fractional 12 steps. Theequations(2.4a)and(2.4c)fornonlinearandviscousstepscanbewritten symbolically with only one component, u, that is, ∂u ∂t =N(u). (2.10a) ∂u ∂t =νL(u). (2.10b) Other components of the velocity field, (v,w), are treated in a similar way. Equa- tion (2.10a) is regarded as a hyperbolic equation and the corresponding penalty formulation for the k th subdomain is ∂u k ∂t =N(u k )−τ k 1 δ i0 [αu k 0 −αu k−1 N ]−τ k 2 δ iN [γu k N −γu k+1 0 ], (2.11) where δ ij is the Kronecker delta tensor with the subscript i corresponding to the collocation point z k i selected for each vertical subdomain. τ k 1 is the penalty coef- ficient at the bottom subdomains and τ k 2 is associated with the top subdomain of 13 k th element. τ k 1 , τ k 2 , α, γ are determined by the vertical interfacial velocities, W k 0 , W k N computed from the previous timesteps, 8 > > > > > > > > > > > > < > > > > > > > > > > > > : W k 0 ·ˆ n o ≥0: α =0, τ k 1 =0, W k 0 ·ˆ n o <0: α =W 0 , τ k 1 = 1 2ω 2 H k , W k N ·ˆ n o ≥0: γ =0, τ k 2 =0, W k N ·ˆ n o <0: γ =|W N |, τ k 2 = 1 2ω 2 H k , (2.12) where ˆ n o is the outward normal unit vector of the interfaces, H k is the height of the k th subdomain, and ω = 2 N(N +1) , (2.13) where N is the number of collocation points in a subdomain. g 1 (t) and g 2 (t) are associated with the velocities of the neighboring points. The penalty method used for viscous step of parabolic nature is intrinsically different than what is described above in a way that it is applied to the Fourier modes in the spectral space, i.e., ∂ 2 ˜ u ∂z 2 − ³ k 2 x +k 2 y − γ o νΔt ´ ˜ u=− ˜ ˆ ˆ u νΔt , (2.14) where ˜ u=F © u n+1 ª , ˜ ˆ ˆ u=F n ˆ ˆ u o . (2.15) 14 The second-order term, ∂ 2 ˜ u/∂z 2 , in the above Helmholtz equation (2.14) is negligi- ble if νΔt << Δz 2 . Without this second-order term, this equation can not match both boundary conditions in the interfaces of the vertical domains and it then be- comes a singular boundary-layer problem with unmatching boundary conditions. There exists solution to this singular perturbation equation only if νΔt≈ Δz 2 . In other words, the equation has to be rescaled to resolve the dynamics within the numerical boundary layer of thickness O( √ νΔt) for high-Re flows. As a result, the corresponding penalty method is formulated with ², ²≡ ³ k 2 x +k 2 y − γ o νΔt ´ −1 . (2.16a) g k 0 (t)=αu k 0 −β² ∂u k 0 ∂z k , g k−1 N (t)=αu k−1 N −β² ∂u k−1 N ∂z k−1 . (2.16b) g k N (t)=γu k N −δ² ∂u k N ∂z k , g k+1 0 (t)=γu k+1 0 −δ² ∂u k+1 0 ∂z k+1 . (2.16c) ² ∂ 2 ˜ u ∂z 2 − ˜ u+² ˜ ˆ ˆ u νΔt −τ k 1 δ i0 [g k 0 −g k−1 N ]−τ k 2 δ iN [g k N −g k+1 0 ]=0. (2.16d) The coefficients α,β,γ,δ are set to 1 at all subdomain interfaces and the physical boundary conditions. The definitions of τ k 1 and τ k 2 can be found in [DDH05]. In the actual simulations, the penalty scheme is said to be turned ”on” or ”off” by multiplying the above two parameters, τ k 1 and τ k 2 , with the Re-dependent ”small” or ”large” numerical numbers. 15 2.2.2 Spectral Filtering Unlike the Penalty method being applied on both spatial and spectral spaces, a low-pass spectral exponential filter is used in both Fourier and Legendre spaces. In Fourier space, the filter function is defined as σ f (k)= 8 > > > < > > > : 1, 0≤k≤k c , exp · −α µ k−k c k max −k c ¶ p f ¸ , k >k c , (2.17) where p f is the order of the filter, k c is the cut-off wavenumber and is set to be 0 throughout our study, k max = p k 2 x,max +k 2 y,max . α =−ln² M denotes the machine presicion. The magnitude of each complex Fourier mode is then ”filtered” by being multiplied with a k-dependent factor σ f (k). σ f (k) of different orders are shown in Figure (2.2). An analogous expression for the filtering operation in Legendre space is expressed as σ l (n)= 8 > > > < > > > : 1, 0≤n≤N c , exp · −α µ n−N c N−N c ¶ p l ¸ , n>N c , (2.18) where n indicates the n th mode in Legendre space, N c is the cut-off Legendre mode and is also set to be 0. N is the total number of Legendre modes. The filtered quantity becomes g F (z i )= N k −1 X j=0 g(z j )σ l (j)l j (z i ). (2.19) 16 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 k/k max σ f (k) p f =6 p f =8 p f =10 p f =20 Figure 2.2: Spectral filters used in Fourier space. See equation (2.17). 2.3 Effective Numerical Viscosity The numerical dissipation is computed as the difference between the energy de- cay rate with minimum numerical dissipation and the energy decay rate with the numerical dissipative mechanism being analyzed. That is, ² n = ½ ∂E ∂t ¾ min−diss − ½ ∂E ∂t ¾ diss . (2.20) The two terms on the right hand side of equation (2.20) are referred to as the baseline and the dissipative cases respectively. Owing to the under-resolution, the baseline case can not be completely free from explicit numerical dissipation for a numericalschemewithnegligibletruncationerror. AminmumamountofLegendre spectral filtering, (P L =20 and P L =10 for Re=5000 and Re=10 5 respectively), 17 is required to stabilize the 8-timestep simulation and to conserve the total kinetic energy for the pseudo-inviscid simulations, i.e., ν/ν o = 0.1 in figure (2.3). Notice that if all the active modes were fully captured by use of adequate grid resolution, theuseofthebaselinecasecouldbesaved. Itisvitaltoensuretheuseofminimum filteringshasnegligibleeffectsontheevolutionofenergyofthesystem. Figure(2.3) shows the decay of total kinetic energy subject to minmum numerical dissipations, i.e., the baseline case. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 t E(t)/E(0) ν/ν 0 =0.1 ν/ν 0 =1. ν/ν 0 =10 (a) Re=5000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.9998 0.99985 0.99990 0.99995 1.00000 1.00005 t E(t)/E(0) ν/ν 0 =0.1 ν/ν 0 =1. ν/ν 0 =10. (b) Re=10 5 Figure 2.3: The decay of total kinetic energy over the course of 8 timesteps with minmum filtering. The Reynolds number Re used in the long-time simulation is defined based on the value of ν o . As is indicated in both figures, the total kinetic energy is almost conserved for an pseudo-inviscid flow, (ν/ν o = 0.1), which justifies the validity of the choice of minimum filtering. Once the numerical value of ² n is obtained, it is normalized with proper factor to obtain the expression of the effective numerical viscosity as described in (1.9). We will now describe how to obtain an equation for 18 the decay rate of an anisotropic spectral energy from which the proper factor can be determined. The discrete Fourier transform of u is defined as F[u]= 1 N x N y N x −1 X m=0 Ny−1 X n=0 u(x m ,y n ,z,t)e −ikxxm e −ikyyn , (2.21) where k x =2πj/L x , k y =2πl/L y , −N x /2≤j≤N x /2−1, −N y /2≤l≤N y /2−1. (2.22) N x and N y are the number of grid points in x and y direction respectively. Under this transformation, the Navier-Stokes equations in the spectral space are written as ∂ ∂t ˜ u n (k,z,t)=N n (k,z,t)+ν µ ∂ 2 ∂z 2 −|k| 2 ¶ ˜ u n (k,z,t)+F n (k,z,t), (2.23) where k≡(k x ,k y ),|k| 2 =k 2 x +k 2 y , (2.24) 19 and N n represents the contributions of nonlinear term and the pressure term. The spectral energy function is defined as E(k,t)= 1 2 ˜ u n (k,z)˜ u ∗ n (k,z). (2.25) Theevolutionequationofspectralenergyisobtainedbyaddingtheequation(2.23) being multiplied by the complex conjugate of ˜ u n , i.e., ˜ u ∗ n , with the complex conju- gateofequation(2.23)beingmultipliedbythevelocityamplitude, ˜ u n . Theresulting evolution equation of spectral energy is ∂ ∂t E(k,z)=R[˜ u ∗ n (k,z)F n (k,z)]+R[˜ u ∗ n (k,z)N n (k,z)]+ν ∂ 2 ∂z 2 E(k,z) −2ν|k| 2 E(k,z)−ν µ ∂ ∂z ˜ u n (k,z) ¶µ ∂ ∂z ˜ u ∗ n (k,z) ¶ . (2.26) R[] stands for taking the real part of the enclosed quantity. Terms on the right hand side of equation (2.26) describe, in the order of occurence, energy production, energy redistribution by nonlinear interactions, energy redistribution by viscous stresses, dissipation caused by the horizontal gradients of the velocity field, and dissipation due to the vertical gradients. The numerical viscosity ² n , calculated from the velocity amplitude ˜ u n is expressed as a function of the magnitude of the wavenumber in spectral space k = p k 2 x +k 2 y through the formula given as ˜ u(k,z)= X |k|− 1 2 Δk≤|k|<|k|+ 1 2 Δk ˜ u(k,z). (2.27) 20 This operation will remove the dependence on k x and k y . As a result, the im- portant physical quantities, ² n , ˜ u n , etc., can then be expressed as functions of the wavenumber magnitude k, that is, ∂ ∂t E(k,z)=R[˜ u ∗ n (k,z)F n (k,z)]+R[˜ u ∗ n (k,z)N n (k,z)]+ν ∂ 2 ∂z 2 E(k,z) −2ν|k| 2 E(k,z)−ν µ ∂ ∂z ˜ u n (k,z) ¶µ ∂ ∂z ˜ u ∗ n (k,z) ¶ . (2.28) Recall that the analogous equation in the case of triply-periodic turbulence, i.e., equation (1.4), is a lot more simpler than equation (2.28), which is different by having only one non-periodic boundary condition. The significance of the equation (2.28) is to provide different formulations of effective numerical viscosity when the boundary conditions are not triply-periodic. In the same spirit as equation (1.9), two expressions of effective numerical viscosity can be defined: ν n,1 (k,z)= ² n (k,z) 2|k| 2 E(k,z) , (2.29) and ν n,2 (k,z)= ² n (k,z) 2|k| 2 E(k,z)+| ∂ ∂z ˜ u(k,z)| 2 . (2.30) Onethingthathastobeclearlystatedisthateventhoughtherearetwoformulaeto quantify numerical viscosity, it doesn’t mean that the total numerical dissipation is different from one formula to the other. Also, the energy spectrum used for 21 normalization in both equation (2.29) and (2.30) are taken to be the baseline case utilizingonlyminimumstabilizerswithoutfurtherindication. ThechoiceofE(k,z) has no particular reason. We also observe that there is no significant change to the energy spectrum E(k,z) over the 8 timesteps of simulation between the baseline runs and the dissipative runs. It is the time rate of change of energy dE/dt that varies significantly from one case to the other. 22 Part I, Chapter 3 Results Theratioofeffectivenumericalviscositytomolecularviscosityisshownatthehor- izontal center plane as well as the adjacent top and bottom subdomain interfaces. Results of two Reynolds numbers, Re = 5000 and Re = 10 5 , are selected for com- parison. Four different stabilizer conditions are listed in table (3.1) for Re=5000. The filter parameters for Re = 10 5 can be found in [DLD08]. In the calculation of numerical viscosity, equation (2.29) is employed unless mentioned otherwise. All the results presented here were calculated at the 5 th of a total 8-timestep run. No significant difference was observed at prior or subsequent timesteps. Table 3.1: Filter Configurations for Re=5000 Cases Fourier Legendre Penalty Symbol Fourier P f =20 P l =20 Minimum o Legendre Off P l =8 Minimum ∗ Penalty Off P l =20 On 4 Total P f =20 P l =8 On − 23 3.1 Horizontal Center Plane Figures (3.1) are the ratio of the effective numerical viscosity to the molecular viscosity ν/ν o at the horizontal center plane, i.e., z = L z /2. As mentioned earlier in chapter (2), penalty scheme is designed to be only effective at the subdomain interfaces, and therefore should have absolutely no effect on the simulated fluid flows in the central region (far away from the subdomain interfaces). This feature has been shown clearly in all the panels in Figure (3.1). It is also very clear that the Fourier filter is the major stabilizing mechanism at the high-wavenumber region whereas the Legendre filter is effective throughout the entire wavenumber space. Curves in the panels of (3.1(a)) and (3.1(b)) were computed using equation (2.29)whereasthecurvesinpanels(3.1(c))and(3.1(d))wereobtainedwithformula (2.30). No qualitative difference has been observed between numerical viscosity calculated using equation (2.29) and (2.30). There are two fundamental differences between the results of Re = 5000 and Re = 10 5 . First, the magnitude of ν/ν o for Re = 10 5 is more than two orders in magnitude larger than that of Re = 5000 in the high-wavenumber region. Secondly, the maximum dissipation occurs at the highest wavenumber in Re=5000 while it occurs at k/k c =0.8 for Re=10 5 . The later is mainly due to the choice of the energy spectrum E(k,z) that is used to calculate the numerical viscosity ν n in equation (2.29) and (2.30). As mentioned earlier in section (2.3), the energy spectrum E(k,z) was chosen to be the one with minmumnumericaldissipation,i.e.,baselinecase. IftheenergyspectrumE(k,z)is 24 replaced by the one corresponding to the dissipative case, this discrepancy between Reynoldsnumberdisappears,seefigure(3.2). Thedistributionshowninfigure(3.2) is what one would expect as having numerical dissipation being only effective at small-scale dynamics. The shift of peak dissipation also has to do with the order of Fourier filter and will be explained in the following subsection. The large variation in magnitude has to be explained by the nature of the Fourier filter whose order is 10 for Re = 10 5 and 20 for Re = 5000. We have designed a ”model simulation” to demonstrate the possible reason for the large variation in magnitude. 3.2 Model Simulation A ”model simulation” is a simulation in which only the Fourier filter is applied. No Navier-Stokes solver is involved. It is just the successive applications of Fourier filteronavelocityfieldwithinitialenergyspectrumas E o (k,t). Thecorresponding non-dissipative configuration is the constant energy. Clearly, after n timesteps, the energy spectrum becomes E(t)=E(nΔt)=E o e −nα( k km ) P f , (3.1) 25 where k m is the maximum wavenumber and P f is the order of the filter. The time derivative of equation (3.1) is · ∂E ∂t ¸ diss = · ∂E ∂n ∂n ∂t ¸ diss = E o Δt · −α( k k m ) P f ¸ e −nα( k km ) P f . (3.2) The corresponding effective numerical viscosities are ν n,1 (k)= 1 Δt · α( k k m ) P f ¸ e −nα( k k m ) P f . (3.3a) ν n,2 (k)= 1 Δt · α( k k m ) P f ¸ . (3.3b) Noticethat(dE/dt) min−diss =0sincetheenergyisconstantforthenon-dissipative case. The comparison between the equation (3.3a) and (3.3b) are shown in Figure (3.3). It clearly demonstrates, at least partially, the reason why the peak dissipa- tion is not located at the highest wavenumber, (equation (3.3a)), when the energy spectrum corresponding to the simulation with minimum stabilizers, i.e., baseline case, isusedtocalculatetheeffectivenumericalviscosity. Furthermore, thesmaller the order of the filter (or the stronger the filter is), the larger the scale on which the filterinfluences (seeFigure(2.2)). Wealso makeacomparison of the numerical viscositybetweenamodelsimulationandaNSsolverwithonlyFourierfilterturned on in Figure (3.4). The effective numerical viscosity is calculated using the energy spectrum corresponding to the baseline case, i.e., minmum dissipation. The differ- 26 ence between two curves on both Figures can be regarded as the effect of nonlinear Navier-Stokes equation solver. The left panel of Figure (3.4) has shown that the effect of the Fourier filter has been enhanced (nearly double) by the Navier-Stokes operator as k/k c ≥ 0.9. The right panel, however, shows that the peak value of the ratio ν n /ν is about 400 with the absence of Navier-Stokes equations solver. It impliesthatν/ν o =400doesn’tliteratelymeanthatthenumericalviscosityisactu- ally 400 times ”stronger” than the molecular viscosity, since there is no ”molecular viscosity” presented in the model simulation at all. We conclude that the large value of ν n /ν is mostly attributed to the nature of the high-order Fourier filter. One last comment about the large value of ν n /ν is that even though the numerical viscosity is extremely strong in the range of high-wavenumber spectrum, most dis- sipated energy in this range is corresponding to the noises of very high frequency and is not likely having much influence on the dynamics of large eddies. 3.3 Interfacial Subdomains The scale-dependent numerical viscosity ν n (k,z) at both bottom and top cen- tral subdomain interfaces is shown in Figure (3.5). It is clear that the effect of Fourier filter is the same as that in the interior region (see Figure (3.1)). It is also shown that when all the stabilizers are applied, the curves behave very differently across Reynolds numbers. The Fourier filter still dominates over other stabilizers in high-Re case as it does in the subdomain interior. However, in low-Re flow (i.e., 27 Re = 5000) all the stabilizers contribute roughly equally and the overall result is the nonlinear interaction of the three. We here explain what appears to be very counter-intuitive about the observation that the numerical viscosity becomes neg- ative in high-wavenumber region when only Legendre filter or penalty scheme is applied. Assuggestedbyitsname”dissipative”,onemightexpectthatthestabilizersare taking energy out of the system but not the other way around. However, the very purpose of introducing of penalty scheme at the subdomain interfaces is to satisfy both boundary and patching conditions simultaneously. This partial enforcement of continuity will increase or decrease the energy of the system depending dynami- cally on what has to be done in order to stabilize the simulation. The enforcement of this weak discontinuity at the subdomain interfaces generates an error which is compatible to the scheme itself [Hes97]. This weak interfacial discontinuity is also suspective to lead to an increase of energy at the highest Legendre modes accord- ing to [Boy01]. This increase of energy is then transformed nonlinearly to other Fourier modes and manifests the overall numerical viscosity to be mostly negative in high-wavenumber region. Finally, for the purpose of accessing the spatial effect of this anti-dissipative character, theeffectivenumericalviscositycorrespondingtothecasewhenonlythe 28 penaltymethodisactiveareshownfromthesubdomaininterfacetofivegridpoints away for both high- and low-Re simulations, see Figures (3.6). The result shows that only two or three grid points away from the interface does the panelty scheme become negligible again. 29 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5 0 5 10 k / k c ν n / ν Fourier Legendre Penalty Total (a) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −200 0 200 400 600 800 1000 1200 1400 1600 k / k c ν n / ν Fourier Legendre Penalty Total 0 0.1 0.2 0.3 0.4 −2 0 2 4 6 8 10 12 14 16 0 < k/kc < 0.4 (b) Re=10 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5 0 5 10 k / k c ν n / ν Fourier Legendre Penalty Total (c) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −200 0 200 400 600 800 1000 1200 1400 1600 k / k c ν n / ν Fourier Legendre Penalty Total 0 0.1 0.2 0.3 0.4 0 5 10 15 0 < k/kc < 0.4 (d) Re=10 5 Figure 3.1: Effective numerical viscosity normalized by the molecular viscosity, i.e., ν/ν o , at the horizontal center plane. k c represents the maximum wavenumber, namely, k c = p k 2 max,x +k 2 max,y . The energy spectrum E(k,z) used in equation (2.29) and (2.30) is corresponding to the baseline case, i.e., minimum dissipation. 30 (a) (b) Figure 3.2: Refer to the caption in figure (3.1). The energy spectrum E(k,z) used in equation (2.29) and (2.30) is corresponding to the dissipative case. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 k/k m ν ν n,1 ν n,2 (a) P f =20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 k/k m ν ν n,1 ν n,2 (b) P f =8 Figure 3.3: The comparison between the results calculated from equations (3.3a), (i.e., the enhanced lines), and (3.3b). 31 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 0 1 2 3 4 5 6 k / k c ν n / ν without NS with NS (a) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −200 0 200 400 600 800 1000 1200 1400 1600 k / k c ν n / ν without NS with NS (b) Re=10 5 Figure 3.4: The comparison of Effective numerical viscosity between the model simulation (without NS) and the actual simulation with only Fourier filter turned on. 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −10 −5 0 5 k / k c ν n / ν Fourier Legendre Penalty Total (a) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −200 0 200 400 600 800 1000 1200 1400 1600 k / k c ν n / ν Fourier Legendre Penalty Total (b) Re=10 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −10 −5 0 5 k / k c ν n / ν Fourier Legendre Penalty Total (c) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −200 0 200 400 600 800 1000 1200 1400 1600 k / k c ν n / ν Fourier Legendre Penalty Total (d) Re=10 5 Figure 3.5: Effective numerical viscosity, ν/ν o , calculated at bottom and top sub- domain interfaces. 33 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5 −4 −3 −2 −1 0 1 2 k / k c ν n / ν Z=0.4583333 dz=5.0875e−04 dz=1.69727e−03 dz=3.54233e−03 dz=6.0133e−03 (a) Re=5000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −6 −5 −4 −3 −2 −1 0 1 k / k c ν n / ν Z = 0.4888888 dz = 4.9698555e−05 dz = 1.663122e−04 dz = 3.487677e−04 dz = 7.070611e−04 (b) Re=10 5 Figure 3.6: The effective numerical viscosities at the neighboring points of the bottom subdomain interface. 34 Part I, Chapter 4 Conclusions When a numerical simulation is under-resolved, implicit or explicit artificial dissi- pation is required to model the uncaptured small-scale physics in order to stabilize the simulation. A Direct Numerical Simulation (DNS) with the introduction of the artificial numerical dissipation, makes it possible to realize a stable simulation whose Re number is in general at least 1 to 2 orders higher than the ones without any stabilizing mechanism. Much more detailed structure of turbulence as well as morerealistic(Reformostrealisticgeophysicalturbulentflowsareextremelyhigh) simulations are possible with this increase of attainable Re by applying artificial dissipation. Therefore,itiscrucialtodeterminetheeffectsoftheemployedartificial numerical dissipation, implicit or explicit, on the actual dynamics of the simulated flow. The concept of effective numerical viscosity proposed by Domaradzki et al. was used to quantify the implicit artificial dissipation embedded in MILES, a DNS solver for homogeneous, isotropic turbulence (triply-periodic boundary conditions) [DX03]. In our study, their primary formulation has been successfully extended 35 to the SMPM code, (Spectral Multi-domain Penalty Method), which has one non- periodic boundary condition in z direction. With this success, we believe that the same procedure can be introduced in principle to quantify any numerical code as long as a stable and ”minimum” dissipative simulation can be defined. The ”minimum” dissipative condition in SMPM is corresponding to the case when only minimum Legendre and Penalty scheme are applied while the Fourier filter is com- pletely absent. The numerical viscosities are calculated using the evolution equation of kinetic energyinspectralspace,i.e.,(2.29)or(2.30). Thisquantityisdefinedinthespectral space and it is normalized by the molecular viscosity used in the actual simulations inallthepresentedfigures. Noqualitativedifferenceisobservedbetweentheresults calculated by equation (2.29) and (2.30). Three explicit stabilizers, i.e., Fourier fil- ter, Legendre filter and penalty scheme have been analyzed with four stabilizer conditions (see table (3.1)). The results show that these three stabilizers are only effective either at the spatial subdomain interfaces or spectral high-wavenumber modes. In the horizontal center plane, the Fourier filter dominates the high-wavenumer regionwhereasthelow-wavenumberdynamicsismostlyinfluencedbytheLegendre filter. Also, for both low- and high-Re simulations, it is proved that the Penalty 36 scheme has negligible effect on the energy content in this region since it is designed to be effective only at the subdomain interfaces. A model simulation is employed to explain the unexpected large ratio of ν n /ν o and the shift of peak dissipation for the high-Re simulation. This ”model simulation” simply takes a initial data and applies the Fourier filter five times successively. In this sense, the calculated numerical viscosity is free from the influence of Navier-Stokes equations solver. The results successfully demonstrate that the order of the filter is the main cause of the large ratio viscosities (see figures (3.3) and (3.4)). In other words, it is not appropriate to regard ν n /ν o =1000 as having the numerical dissipation to be 1000 times larger than the ”physical” molecular dissipation. The unexpected large ratio is mainly originated from the order of the Fourier filter as well as the lack of a scale-dependent molecular viscosity which is only active at small-scale dynamics. The shift of the peak dissipation has to do with the order of the Fourier filter as well as the energy spectrum, i.e., E(k,z), that is used to calculate the effective numerical viscosity (see figures (3.3) and (3.4)). In the subdomain interfaces, the Fourier filter still dominates the region of high-wavenumber as it does in the interior region. The effect of Legendre filter is not negligible compared to the Fourier filter at high-wavenumber region. Both the Penalty scheme and the Legendre filter show the unexpected anti-dissipative character, i.e., the stabilizers are not removing the energy out of the system but adding the energy into it. The anti-dissipative character can be attributed to the 37 discontinuous nature of Penalty scheme at the subdomain interfaces. It appears that the partial enforcement of the continuity at the subdomain interface provides the Penalty scheme the freedom to either increase or decrease the energy content in order to avoid the Gibbs’ oscillations and stabilize the simulations. We suspect that this anti-dissipative effect in the ”Legendre” case (see table (3.1)) is somehow magnifiedbyLegendrefilterthroughthenonlinearinteraction. ThePenaltyscheme quickly becomes negligible only a few points away from the subdomain interfaces. 38 Part II - The Wave-Induced Boundary Layer Under Internal Solitary Waves Part II, Chapter 1 Introduction 1.1 Internal Solitary Waves When two opposite effects, namely, nonlinearity and dispersion balance each other, apropagatingwavewithpermanentwaveformispossible. ThefamousKorteweg-de Vries equation (KdV), derived by Korteweg and De Vries in 1895 in their study of long surface gravity waves in a rectangular canal [KdV95], governs the evolution of waves of weakly-nonlinear (small but finite amplitude), weakly dispersive nature. The KdV equation and the corresponding solitary wave solution are, 39 η t +(c+α 1 η)η x +βη xxx =0, η(x)=η o sech 2 ( x−Vt β 2 ), (1.1) where the coefficients α 1 ,α 2 ,β,η o , c, β 2 and V corresponding to a plane wave in a two-fluid, internal wave system can be found in [OS05]. Characteristic of a nonlin- ear wave, the wave amplitude η o , enters the definition of both the wave length β 2 and the wave speed V. One of the features of the solitary wave solution admitted by equation (1.1) is that the properties (e.g., waveform, speed, etc.) are preserved upon interaction with one another. Such dynamics is granted the name ”soliton” for this almost particle-like behavior [JD65]. KdV solitons have been observed in many nonlinear physical systems: optical-fiber solitons, lattice dynamics with non- linear atomic interactions, blood flow in arteries, internal solitary waves, just to nameafew. Thecurrentstudyfocusesonnonlinearinternalwavesobservedonthe continental shelf or in stratified lakes, and particularly on the boundary layer that must form along the bottom surface under the footprint of such waves. Inenvironmentswherethe stratification issuchthat the pycnoclineis relatively thin compared to the depth of the fluid column, long nonlinear waves appear either as waves of elevation or waves of depression. In a wave of elevation the pycnocline is disturbed upward into the upper mixed layer as it propagates. In a wave of depression the pycnocline is distorted downward into the lower region of the water 40 column. For long waves in a two-layer stratification, the profile shape always takes the form of a wave of elevation with respect to the most shallow layer. Thus, if the upper layer is shallower than the lower layer, the profile shape will be a wave of depression. Alternatively, if the lower layer is the controlling depth (most shallow layer),thewavewillbeoneofelevation. Thedistinctionbetweenwavesofelevation and depression is very important when it comes to the development of the wave- induced boundary layer. The distinction between wave types is important because it determines the streamwise pressure field impressed on the benthic region by a propagating wave. Extremely energetic packets of internal solitary waves of depression were ob- served over a three-week period on a shallow and strongly stratified pycnocline off Northern Oregon by Stanton and Ostrovsky [SO98]. Even though many studies, e.g. [Lam94], suggest that the steep-edged internal tidal bore is prone to disinte- grate into packets of weakly or moderately nonlinear solitons by the presence of weak dispersion, Stanton and Ostrovsky measured a highly nonlinear soliton train. AstheKdVequationrepresentsalow-orderasymptoticapproximationofnonlinear propagation of long waves in a wave guide, larger amplitude dynamics can be mod- eledbyanextendedKorteweg-deVries,(eKdV)equation. Includingasecond-order approximation for the nonlinear advection, the eKdV equation has the form η t +(c+αη 1 +α 2 η 2 )η x +βη xxx =0. (1.2) 41 The additional second-order nonlinear term accounts for the strongly-nonlinear ef- fects that the regular KdV equation (1.1) does not have. Stanton and Ostrovsky’s observations show the necessity to research solitons of more than weakly-nonlinear form, and motivates use of even a fully-nonlinear wave profile. An investigation of fully-nonlinear, long internal waves in a two-fluid system was undertaken by Choi and Camassa [CC99]. Asymptotic approximations were used to formulate evolution models for weakly-nonlinear and fully-nonlinear waves. These authors compared their asymptotic results with solutions of the full Euler equations. The nonlinear wave profiles used in the present study are based on an alternate approach introduced by Sakai and Redekopp [SR07]. They developed a KdV model that employs an exact (i.e., fully nonlinear) expression for the speed of long waves moving along a single characteristic coordinate. Details of their work is deferred to section (2.1.1). The particular objective here is to compute the fully- nonlinear dynamics of the viscous boundary layer that must necessarily develop in the benthic layer under strongly-nonlinear, internal solitary waves in shallow seas. The propagation of a soliton near the bottom of the continental shelf is likely to induce large bottom shear stresses. Furthermore, as a wave form necessarily in- volves a region of adverse pressure gradient, there is potential for the wave-induced boundary-layer to separate as it moves over the seabed. An energetic soliton (i.e., 42 large amplitude wave form) is expected to trigger transition to turbulence in the wave-inducedboundarylayer,aswellasotherpronounceddynamicsassociatedwith instabilities in the separating boundary layer. Such dynamics may even contribute to episodes of marked sediment resuspension. Further aspects of the instability of separated boundary layers and the potential for resuspension and mixing are now reviewed briefly. 1.2 The Wave-Induced Boundary Layer The observation of internal solitary waves propagating upstream along a strongly stratified bottom layer below a well-mixed upper layer on the California shelf was reported by Bogucki et al. [BDR97]. The observation site was the Palos Verdes shelf where it is estimated that about 200 tons of DDT remains mixed in sediment deposits[NY91]. Thedetailedphysicalprocessformeasuredsedimentresuspension from this site is quite complicated and might result from several different mecha- nisms. Both the near-bottom steady currents and the long internal waves moving over the sedimentary layer are likely to impose shear stresses on the seabed that are larger than the critical value for the dislodging of particulates and, as a conse- quence, resuspension might occur [NW94]. Opposite to the traditional belief that sediment resuspension occurs mostly when the immediate steady current outside 43 theboundarylayerexceedsathresholdvalue(e.g., [DM86]),thesoliton-inducedre- suspensionaccountsforamarkedincreaseofparticulatematterwhenthelow-speed bottom current is still very quiet. Bogucki et al. obtain a good correlation between theincreased waterturbidityowingto the increase of resuspension and the passage of a soliton train with three solitary waves. Their results support the assumption that the resonantly-generated internal solitary wave is indeed an effective agent for theobservedsedimentresuspension and mixing of the bottom boundary-layer flow. Their simulations show that up to 73% of the energy of the internal wave field is carried upstream by the internal solitary waves. The unsteady dynamics of boundary-layer flow under the footprint of a long in- ternalwaveofelevationgeneratedbytopographicresonance(i.e.,resonancebetween a mean current and weak topographic forcing) was studied by Wang and Redekopp [jWR00]. A two-layer stratification was used to model the flow over a localized to- pographicforcingtoexploreamechanismbywhichthebottomboundary-layerflow ispronetoseparateandtotriggertheonsetofaspontaneous(global)instability. In theirstudythecreatedlongwavepropagatedupstreamagainstashearedcurrent,a flowbasedonobservationsreportedbyBoguckietal. [BR99]. Whentheamplitude of the wave exceeded a critical value, the induced reversed flow in the separation bubble under the wave became sufficiently strong that a global instability ensued, causingeruptionsoftheseparatedregionandafiner-scalestructureofconcentrated 44 vortices. Duetotheslowly-movingandspatially-extensivenatureofthelongwave, the induced synchronous dynamics of the global instability creates an enhanced spatio-temporal coherence of the flow field near the bottom layer over a significant spatial range. According to Bogucki and Redekopp [BR99], the gradients of both normal and tangential stresses obtained on the bottom surface are increased ten- fold over the region associated with the global instability. This dramatic increase in surface stress is more likely to be related to some intrinsic mechanism triggered by the onset of global instability rather than any bursts and sweeps generated by thefluctuatingstressfieldobservedinatypical,attachedturbulentboundarylayer. Bogucki and Redekopp [BR99] employed a 2-D direct numerical simulation en- compassing a domain between the lower portion of the waveguide and the bottom wall in order to probe the stress field within the boundary layer induced by the propagation of a soliton of elevation. This configuration closely mimics the flow state observed in [BDR97]. The solution of the weakly-nonlinear KdV equation governing a soliton propagating along the density interface in a two-layer stratifi- cation is used to determine the boundary conditions on the top of computational domain. When the amplitude of the soliton reached a critical value, intrinsic and spontaneousdynamics(i.e.,globalinstability)abruptlyappeared. Thesynchronous nature of the global instability was manifested through the coherence of the power 45 spectra of the horizontal velocity taken at two different locations within the sepa- ration bubble. Bogucki and Redekopp proposed [BR99] that the prime mechanism for the resuspension observed in [BDR97] involve a discrete frequency pulsation plus the accompanying reduction of the horizontal scale of the coherent dynamics in the separated region under the wave. Theboundary-layerflowinducedbyaweakly-nonlinearsolitarywaveofdepres- sion was examined by Diamessis and Redekopp [DR06] using a DNS code with high spectral accuracy and fine-scale resolution near the benthic boundary. Their results showed that, as the wave amplitude exceeds a critical value, a coherent and spontaneous oscillation extends considerably beyond the height of the separation bubble. This synchronous dynamics, resulting from global instability, causes the separation bubble to break up as it grows. The result is that vortices are shed upward into the water column. Their discovery for waves of depression, as well as that of Bogucki and Redekopp [BR99] for waves of elevation in a sheared current, indicate that it is this synchronous mechanism relating to global instability that forms the flow structure that is most likely to contribute to sediment resuspension. Carr and Davies [CD06] investigated the velocity field induced by the propaga- tion of a large-amplitude solitary wave of depression in a laboratory tank following the ideas put forward by Diamessis and Redekopp. Their measurements showed some semblance of agreement with Diamessis and Redekopp’s simulation, but the 46 method used to form the large amplitude wave initiated strong disturbances near the bottom boundary layer in the tank, and this contaminated the flow and made correspondence with the simulation problematic. With the above reviews as motivation, this report is aimed toward understand- ing the turbulent dynamics of the boundary-layer flow under the footprint of a solitary wave of depression. Properties such as boundary-layer thicknesses, wall shear stress, the separation point, etc., are studied with the solitary wave ampli- tude and boundary layer Reynolds number as parameters. A linearized theory for the wave-induced boundary layer formed under weakly-nonlinear solitary waves described by the first-order KdV equation has been developed by Ott and Sudan [OS70] and Miles [Mil76]. This theory requires that the wave-induced flow (par- ticle speed) differs only marginally relative to the wall. These analyses have been extended by a weakly-nonlinear boundary layer analyzed by Liu [LPC07]. Neither of the approaches can address the issue of transition to turbulence in the wave- induced boundary layer nor the issue of boundary layer separation. The present work is aimed at finding the boundary layer drag exerted on propagating internal solitary waves with reasonable accounting for the onset of transition to turbulent flow. This objective necessarily involves a direct numerical simulation with turbu- lence modeling. 47 In the next chapter the properties of a family of fully-nonlinear solitary wave solutions of the Euler equation are discussed. This is followed by a description of the numerical model and a summary of computed results. 48 Part II, Chapter 2 Problem Definition 2.1 Base Flow Theinviscidflowfieldinducedbyaleft-runninginternalsolitarywaveofdepression in a two-layer system is shown in Figure (2.1). The components of the dimensional velocity field, u(x ∗ ) and v(x ∗ ,y ∗ ), are referred to a steady reference frame that is fixed to the moving wave. Relative to this moving frame, the solitary wave is stationary,andthetopandbottomboundariesofthewaveguide,plustheupstream water column, are moving uniformly to the right at the wave speed C ∗ w . That is, u(x ∗ = ±∞) = C ∗ w . The wave profile itself serves as a dividing line between the two immiscible layers. As a result, the entire lower-layer appears as a convergent- divergent channel in which the flow undergoes accelerating-decelerating phases to conserve the mass flow rate. The solitary wave profile is derived in section (2.1.1), followedbytheexactformulaeforthevelocityfield, u(x ∗ )andv(x ∗ ,y ∗ ). Thereader 49 ρ 1 ρ 2 h 1 h 2 u(x * ,0) = C * w v(x * ,y * ) u(x * ) η o η(x * ) u(x * =−∞)=C * w Figure 2.1: A definition sketch of a stationary solitary wave and its boundary conditions. The depth ratio, h 2 /h 1 = h, is greater than 1 for a depression wave. Forstablestratification,ρ 2 >ρ 1 . Thewave-inducedvelocityfieldsinthestreamwise and transverse direction are u(x ∗ ) and v(x ∗ ,y ∗ ) respectively. u, v, x ∗ , y ∗ , η, h 1 , h 2 , ρ 1 and ρ 2 are all dimensional variables. willnotethat,consistentwithshallowwatertheory(alt.,longwaves),thehorizontal velocity component is a function of the streamwise coordinate alone. 2.1.1 Solitary Wave Profile In long-wave theory, the wave length of any motion of interest must be long com- pared to the controlling length scale. For internal wave motions the controlling length scale for waves of depression is the upper mixed layer depth. In the case of the two-layer model for the density stratification used here, that length scale is h 1 . It is for this reason that we adopt h 1 as the length scale in forming a non- dimensional representation of the motion of the interface depicted in Figure (2.1). Also, we use the linear phase speed for infinitely long waves, c o , as a velocity scale. For the flow model studied here, the value of this reference speed is given by 50 c 2 o = g(ρ 2 −ρ 1 )h 1 h 2 ρ 2 h 1 +ρ 1 h 2 ∼ = ˜ g h 1 h 2 h 1 +h 2 , ˜ g =g( ρ 2 −ρ 1 ρ 1 ). (2.1) The later approximated expression gives the value of c o in the limit of the Boussi- nessq approximation. Using these reference scales, a dimensionless amplitude ζ is introduced, along with a dimensionless streamwise coordinate x and time t, by the following definitions: ζ = η h 1 , x= x ∗ h 1 , t= c o t ∗ h 1 . (2.2) In addition, we define the dimensionless depth ratio h= h 2 h 1 , (2.3) which becomes the primary dimensionless parameter in the inviscid solution which follows. In these expressions x ∗ and t ∗ represent dimensional independent variables for distance and time, respectively. The governing equation for long-wave motions on the interface is based on the analysis by Sakai and Redekopp, [SR07]. Essential elements of that analysis are reviewedinthissection. ThestartingpointistheextendedKdVequationbasedon the exact nonlinear phase speed c E (ζ) for infinitely long, plane waves moving along 51 a single characteristic in a two-layer fluid model. That is, the interface evolution is described by the equation, ζ t +c E (ζ)ζ x +β o {β(ζ)c E (ζ)ζ xx } x =0. (2.4) The nonlinear phase speed c E was derived by Ostrovsky and Stepanyants [OS05] and, when scaled by the reference speed c o , is given as c E =1+3 h 2 −h 1 (h 2 +h 1 ) 2 (h 2 −h 1 −2η) 8 < : s (h 1 +η)(h 2 −η) h 1 h 2 − h 2 −h 1 −2η h 2 −h 1 9 = ; (2.5a) =1+3 h−1 (h+1) 2 (h−1−2ζ) r 1 h (1+ζ)(h−ζ)−3( h−1−2ζ h+1 ) 2 . (2.5b) Sakai and Redekopp discuss several strongly-nonlinear models based on different polynomial approximations to the fully-nonlinear wave speed c E (ζ). Models with polynomial approximations are useful as they can be employed to derive conser- vation laws of the evolution system. Such conserved densities are often useful in analysis, and are important to assessing the accuracy of numerical approximations. In this study, however, we choose to use equation above without approximation. 52 The dispersive parameter β o and the nonlinear dispersive function β(ζ) appear- ing in the extended KdV equation are given by the following expressions: β o = 1 6 h 1 h 2 h 2 1 = 1 6 h, β(ζ)= (h 1 +η)(h 2 −η) h 1 h 2 = 1 h (1+ζ)(h−ζ). (2.6) The reader should note that the interface displacement function η (alt., ζ) is taken to be positive in the downward direction since we consider exclusively waves of de- pression. Hence,positivevaluesofthewaveamplitudecorrespondtotheamplitude of waves with a depressed thermocline. The mathematical model presented above admits a family of solitary wave so- lutions. Since these solutions will be used later to define the inviscid flow field adjacent to the bottom (benthic) surface, the profile characteristics of the family of solitary wave solutions are developed here. Such solutions of (2.4) have the functional form ζ(x,t)=Z(X)=Z( x−Vt √ β o ). (2.7) Substituting this form into the evolution equation (2.4), one obtains {β(Z)c E (Z)Z XX } X ={(V −1)−(c E (Z)−1)}Z X . (2.8) 53 Considering the case of a solitary waveform such that Z → 0 as X → −∞, and integrating once, a nonlinear eigenvalue problem emerges for the waveform Z(X): Z XX = λZ−UZ β(Z)c E (Z) =R(Z), (2.9) where the function U(Z) is U(Z)= Z Z 0 {c E (Z)−1}dZ. (2.10) In the foregoing equation, the eigenvalue λ is related to the solitary wave speed V through the expression λ=V −1. Before discussing solutions of the eigenvalue problem we note that equation (2.9) can be multiplied by dZ/dX and written in the form d dX ½ 1 2 ( dZ dX ) 2 − Z Z R(θ)dθ ¾ =0. (2.11) Noting that dZ/dX vanishes at the center of a solitary wave, the equation can be integrated from the location at X = −∞ where the wave amplitude is Z = 0 to the center of the wave X = 0 where the wave amplitude has its maximum value 54 Z = ζ o . The result is the following integral relation for the eigenvalue λ in terms of the wave amplitude ζ o (i.e., the peak displacement of the interface): λ= q 1 (ζ o ) q 2 (ζ o ) , q 1 (ζ o )= Z ζo 0 U(Z) β(Z)c E (Z) dZ, q 2 (ζ o )= Z ζo 0 Z β(Z)c E (Z) dZ. (2.12) It is evident from this discussion that the wave amplitude ζ o depends on the sim- ple environmental parameter h = h 2 /h 1 , when the Boussinessq approximation is invoked. Herein, Z and ζ are interchangable without any confusion. 2.1.2 Solitary Wave Characteristics Before considering basic aspects of the boundary layer flow under a solitary wave, andthenumericalsimulationofthatwall-boundedflow, someessentialcharacteris- ticsofthewavefielddrivingtheboundarylayerarereviewed. Thewaveprofilesfor two typical solitary waves of depression at a fixed depth ratio are shown in Figure (2.2). These profiles are obtained as solutions of the eigenvalue problem described in section (2.1.1). The reader will note that the width of a solitary wave in this family varies considerably with amplitude. The properties of a host of solutions of this type are summarized in order to define the general nature of the ambient flow near the bottom surface. The properties of the flow at this level define the character of the boundary layer that develops along the bottom surface. We adopt the perspective of singular perturbation theory where the outer flow, given by the 55 inviscid solitary wave solution in equation (2.16,2.17), is considered first, and then the inner boundary layer flow driven by this outer flow is examined through nu- merical simulation. −100 −50 0 50 100 0 0.1 0.2 0.3 0.4 0.5 (x−x c ) ζ ζ o /ζ oe = 0.5/3.42 −100 −50 0 50 100 0 0.5 1 1.5 2 2.5 3 3.5 ζ o /ζ oe = 3.4/3.42 (x−x c ) ζ Figure 2.2: Solitary wave profiles for depth ratio h = 8. x c is the center of the wave. Thefirstwavepropertyweconsideristhevariationofthepeakwaveamplitude, denotedbyζ o ,asafunctionofthedepthratioh. Foraspecificdepthratio,thepeak amplitude of a solitary increases until a critical (alt., extreme) value ζ oe is reached. As this critical amplitude is approached, the width of the solitary wave begins to increaserapidly,andincreasesindefinitelypreciselyatthecriticalvalue. Thiseffect isevidentintheprofilesexhibitedinFigure(2.2). Atthisextremewaveamplitude, the fully-nonlinear, inviscid equations of motion yield a solution resembling a finite jump of the interface. This solution connects two distinct parallel flows, two flows referred to as conjugate states. This critical, or extreme, amplitude ζ oe is dictated 56 by equation (2.9), and its value as a function of the depth ratio is shown in the left panel of Figure (2.3), as well as in the Table (2.1). It is evident from this figure that the extreme wave amplitude varies almost linearly with the depth ratio. The solitary wave speed is C w = 1+λ and reaches an extreme value when the amplitude ζ o reaches its extreme value ζ oe . Denoting this extreme speed as C oe , we exhibit its variation as a function of the depth ratio in the right panel of Figure (2.3). The wave speed C w for any depth ratio increses monotonically from unity to the extreme value C oe as the amplitude ζ o increases from zero to its extreme value ζ oe . The variation of C oe is near linear, increasing from a value of unity at the depth ratio of h=1 to values that are substantially greater than the reference value of linear, long waves at the same depth ratio. When the depth ratio is unity, and the Boussinesq approximation is invoked, the linear solution becomes exact and the wave speed corresponds precisely to the linear long-wave speed. h 2 3 4 5 6 7 8 9 10 ζ o 0.5 0.98 1.47 1.96 2.45 2.93 3.42 3.91 4.4 Table 2.1: ζ o : dimensionless extreme wave amplitude An important solitary wave property relevant to the development of the wave- induced boundary layer is a characteristic measure of the wave length (alt., the 57 0 2 4 6 8 10 0 1 2 3 4 5 h ζ oe 2 4 6 8 10 1 1.2 1.4 1.6 1.8 h C oe Figure 2.3: The variation of the extreme wave amplitude and the corresponding wave speed with respect to the depth ratio. solitary wave width). In this work we introduce the characteristic wave length L w defined as L w = Z ∞ −∞ ζ(x) ζ o (x) dx. (2.13) NotethatthewavelengthL w isscaledwiththeupperlayerdepthh 1 . Thevariation of L w with respect to the amplitude ζ o and the amplitude ratio ζ o /ζ oe is depicted in Figure (2.4). The characteristic wave length first decreases as the wave ampli- tude increases, consistent with first-order KdV theory, but then increases again as the critical amplitude is approached. Also, for a fixed amplitude ratio, waves in environments with a larger depth ratio are, in general, broader. These proper- ties are evident in the left panel of Figure (2.4). However, when the amplitude is scaled by the critical value for that depth ratio, a considerably more universal vari- ationofthewavelengthisrealized. ThisisshownintherightpanelofFigure(2.4). 58 0 1 2 3 4 10 12 14 16 18 20 22 ζ o L w 0 0.2 0.4 0.6 0.8 1 10 12 14 16 18 20 22 ζ o /ζ oe L w h=3 h=6 h=8 Figure 2.4: The variation of characteristic wave length with respect to the wave amplitude ζ o (left panel) as well as the wave amplitude ratio ζ o /ζ oe (right panel). 2.1.3 Inviscid Flow Field We are now in a position to calculate the inviscid flow field induced by a solitary wavewhoseprofileisdefinedby(2.9). Asnotedearlier, theflowintheentirelower- layer,whenviewedinawave-fixedframe,issimilartothatinaconvergent-divergent channel. Conservation of mass in this channel flow requires that the streamwise velocity satisfies the relation u(x ∗ )(h 2 −η(x ∗ ))=C ∗ w h 2 . (2.14) Thevelocityfieldinthetransversedirectioncanthenbeobtaineddirectlyfromthe conservation of mass constraint through integration. Anticipating the later boundary layer calculation, we re-cast this equation in an appropriate dimensionless form. To this end we note that the linear long-wave 59 phasespeedc o dependsonthe depthratio h. Inorder to revealthefull dependence ofthevelocityfieldonthedepthratiohexplicitly,wehenceforthdefineanalternate reference speed c oo where c 2 oo =g( ρ 2 −ρ 1 ρ 1 )h 1 = ˜ gh 1 , (2.15) where ˜ g isthereducedgravityfortheinterfacialconfiguration. Thisreferencespeed dependssolelyonthereferencelengthh 1 , andthedimensionlessexpressionsforthe velocity field nowreveal explicitly their full dependence on the depth ratio h. With thischangein velocityscaling, theambient, inviscid velocityfield in the lower layer is given by the component relations U ∞ (x)= u(x) c oo =C w r h 1+h µ h h−Z(X) ¶ , (2.16) V ∞ (x,y c )= v(x,y c ) c oo =−C w r h 1+h µ h (h−Z(X)) 2 ¶ dZ dX y c , (2.17) For purposes of defining the ambient flow in the calculation of the boundary layer adjacent to the wall, an upper limit of the distance from the wall, say y = y c , is chosen with y c much larger than the boundary-layer thickness δ(x) for all x. Refer to Appendix (B.1) for detailed derivations for equation (2.16) and (2.17). 60 Another quantity of primary relevance to the development of a boundary layer is the streamwise pressure gradient. In Prandtl’s original conception of the bound- ary layer, and in the asymptotic sense of inner and outer expansions valid in the limit of large Reynolds numbers, the pressure gradient is prescribed by the inviscid outer flow. In this case, the solitary wave profile determines the inviscid flow field near the bottom wall, and thus it sets the pressure gradient. It is emphasized that the slope of the wave profile, dζ/dx (alt., dZ/dX), relates directly to the pressure gradientoftheinviscidflowfieldand,therefore, playsadecisiveroleindetermining the physical characteristics of the induced boundary-layer flow. The relation between the pressure gradient and the solitary wave profile is derived inAppendix(B.2). Thedimensionlessexpressionfortheexternalpressuregradient is ∂P ∂x =− (C w h) 2 (h−ζ) 3 dζ dx . (2.18) This expression can be rewritten in the alternate form ∂P ∂x =− ζ oe C 2 w h 2 [h−ζ oe (ζ/ζ oe )] 3 d(ζ/ζ oe ) dx . (2.19) Since ζ oe is a unique function of h (i.e., an essentially linear function), and since the wave speed C w depends solely on h and ζ o /ζ oe , the maximum pressure gradient necessarily scales directly in terms of these two variables: the depth ratio h and 61 the amplitude ratio ζ o /ζ oe . −8 −7 −6 −5 −4 −0.1 −0.05 0 0.05 0.1 (x * −x * c )/L * w ∂ P/∂ x (a) 3 6 8 −8 −7 −6 −5 −4 −0.1 −0.05 0 0.05 0.1 (x * −x * c )/L * w ∂ P/∂ x (b) 0.292 0.438 0.584 Figure 2.5: (a).ζ o /ζ oe ∼ 0.5. The legend shows the depth ratio h = h 2 /h 1 . (b). h=8. ; The legend shows the amplitude ratio ζ o /ζ oe . The spatial variation of the pressure gradient for different wave amplitudes are compared in Figure (2.5). The pressure gradient is negative (favorable) in the front half of the wave, passes through zero in the center of the wave, and becomes positive (adverse) in the back half of the wave. The boundary layer begins with zero velocity difference and zero pressure gradient, moves to a maximum velocity difference at the center where the pressure gradient vanishes, and then moves to a region of decreasing velocity difference and opposing pressure gradient aft of the wavecenter. Theregionoffavorablepressuregradientgenerallyinhibits,oratleast delays, transition to turbulence. Hence, it is expected that the boundary layer will tend to be laminar in the front portion of the wave. As the Reynolds number increases along the flow direction, the boundary layer under the front of the wave 62 may reach a state where the boundary layer undergoes transition, even though the pressure gradient is favorable. However, as transition is generally encouraged in regions of adverse pressure gradient, onset of turbulent motion is expected to almost always occur somewhere along the aft portion of the wave, even when the boundary layer is laminar at the wave center. Considerations of this sort dictated theselectionofaturbulencemodelingcodeasthetoolforstudyingthedevelopment of this wave-induced boundary layer. 2.1.4 Boundary layer Basics As to the variation of the speed of the solitary waves, another velocity measure in addition to the wave speed C w defined previously is of prime importance. In particular, we refer to the magnitude of the velocity difference across the boundary layer (alt., the strength of the driving stress for the boundary layer). This measure is given by the relation ΔU(x)= ΔU ∗ (x) c oo = u ∗ −C ∗ w c oo =C w r h 1+h ζ(x) (h−ζ(x)) . (2.20) 63 −5 0 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 (x * −x * c )/L * w Δ U (x) 0.5 2.0 3.3 Figure 2.6: Refer to equation (2.20). h=8, ζ o =0.5,2.0,3.3. A typical variation of ΔU(x) under the wave is shown in Figure (2.6). The figure shows, forinstance, thattheinviscidvelocitynearthecenterofthewaveisasignif- icant fraction greater than the wave speed, approaching even 100% for the highest wave amplitude case considered in Figure (2.6). This demonstrates that waves of large amplitude induce ambient flow states with velocities quite large relative to the wall. As the velocity difference increases from the far left to the center of the wave, the boundary-layer flow develops accordingly, starting from a state of zero difference at the front of the wave, reaching a maximum at the wave center, and then diminishing again to zero difference downstream of the wave. After the flow passes the center of the wave, it slows down owing to the presence of an adverse pressure gradient. At any cross-section normal to the bottom boundary, ΔU(x) serves as the true effective velocity scale for the boundary layer at any station. 64 Asimpliedbythedefinitionoftheeffectivevelocity,equation(2.20), theboundary- layer flow is driven by the velocity difference between the wall and the inviscid flow fieldinducedbythesolitarywave. Usingequation(2.20),wedefinethepeakdriving stress (velocity difference) across the boundary layer as ΔU c =ΔU(x c ), (2.21) where x c denotes the position of the wave center. The magnitude of this driving velocity difference, scaled relative to the wave speed C w , is shown as a function of the velocity ratio for several values of the depth ratio in Figure (2.7). The wavecenterisparticularlyinterestingasitdivides theupstream region offavorable pressure gradient from the downstream region of adverse pressure gradient. At the wave center, the pressure gradient is identically zero. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ζ o /ζ oe Δ U(x c ) 2 3 6 8 Figure2.7: Thevelocitydifferenceunderthepeakofthewave. Seeequation(2.21). It isconvenient todefineaReynolds number by using thepeak drivingvelocitydif- ference found at the centerof the wave; that is, we define the wave-center Reynolds 65 number Re w,c = ΔU ∗ (x c )L ∗ w ν =ΔU(x c )L w Re o =F wc Re o . (2.22) The Reynolds number Re o is based on the environmental parameters c oo and h 1 , which are fixed for any simulation (See Appendix (B.1)). It is defined as Re o = c oo h 1 ν . (2.23) Incontrasttothewave-centerReynoldsnumber, awaveReynoldsnumber Re w can bedefinedbasedonthewavespeed(i.e.,theoncomingstreamspeedinawave-fixed frame) by the relation Re w = C ∗ w L ∗ w ν =C w L w r h 1+h Re o =F w Re o . (2.24) ThefactorsF wc andF w relatingtherespectiveReynoldsnumbersto Re o areshown inFigure(2.9). Thewave-centerandwaveReynoldsnumberaresignificantlylarger thantheReynoldsnumberbasedonenvironmentalconditions. Ingeneral,thewave Reynoldsnumberisconsiderablylargerthanthewave-centerReynoldsnumber,ex- cept when the wave amplitude becomes large. This is shown in Figure (2.8) where the wave-center Reynolds number is plotted as a function of the wave amplitude 66 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 x 10 7 ζ o /ζ oe Re w,c h=2 h=3 h=6 h=8 Figure 2.8: Re w,c :the wave-center Reynolds number, defined in equation (2.22). 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 ζ o /ζ oe F wc 0 0.2 0.4 0.6 0.8 1 2 3 4 5 6 7 ζ o /ζ oe F w h=2 h=3 h=6 h=8 Figure 2.9: The variation of F w and F wc with respect to the amplitude ratio. F w and F wc are defined in equation (2.22) and (2.24). ratio for a wave Reynolds number Re w fixed at the value of 10 8 . The wave-center Reynolds number Re w,c is the preferred indicator of the de- velopment of the boundary-layer; for example, such properties as boundary-layer thickness, skin friction coefficient, transition, etc. In addition to the two Reynolds numbersdefinedabove, wedefineanalternatewave-centerReynoldsnumberwhere 67 the length scale is taken as the boundary layer momentum thickness δ ∗ 2 . This is the most relevant Reynolds number for characterizing the onset of transition, and other turbulence quantities, and is defined as: Re δ 2 ,ΔU c = ΔU ∗ c δ ∗ 2 ν =ΔU c δ 2 Re o . (2.25) In this relation δ 2 is the dimensionless momentum thickness of the boundary layer. The displacement thickness δ 1 and the momentum thickness δ 2 are defined by the dimensionless relations δ 1 (x)= Z y s 0 Ã 1− u(x,y)− c C w u(x,y s )− c C w ! dy, (2.26) δ 2 (x)= Z ys 0 Ã u(x,y)− c C w u(x,y s )− c C w !Ã 1− u(x,y)− c C w u(x,y s )− c C w ! dy, (2.27) wherelengthsarescaledbyh 1 andvelocitiesbyc oo . Thewavespeed c C w is,therefore, related to C w by c C w = r h 1+h C w . (2.28) Recall that u(x,y) is measured from a reference frame fixed to a left-running wave. The level y s (x) is a location sufficiently above the boundary layer so that the hor- izontal speed at this level u(x,y s (x)) matches the ambient speed. Of course, since the flow field (i.e., the (inner) boundary layer plus the (outer) external flow) in 68 0 < y < y s is computed via the Navier Stokes equations (with turbulence mod- eling), the computed outer flow does not necessarily correspond precisely to the inviscid flow prescribed by the solitary wave solution. That is, u(x,y s ) does not necessarily equal the value given by U ∞ (x). It is for this reason that we use the computed speed u(x,y s ) as the ambient reference speed in computations of the dis- placement and momentum thickness measures of the boundary layer profile. Atwo-dimensionalNavier-Stokesequationsolverwithturbulencemodel,mainly k−ω model, is then utilized to find a statistically steady-state solution of the flow field within the bottom boundary layer. The computational domain is setup to be a rectangular box of size L ∗ c ×y ∗ c measured from the bottom wall as outlined by dashline in figure (2.10). For all computational results presented in this study, the length of the computational domain was fixed by using k 1 = k 2 = 6. These values were chosen both for optimal resolution and computational efficiency. 69 The boundary conditions applied for the numerical simulation are, u=C w , v =0, @y =0. (2.29a) u=U ∞ (x), v =V ∞ (x,y c ), @y =y c . (2.29b) u=C w , v =0, @x=0. (2.29c) ∂u ∂x =0, ∂v ∂x =0, @x=L c . (2.29d) 2.2 Turbulent Model Owing to the fact that many long, internal waves appearing in shallow waters, as on the continental shelf and in lakes, are strongly nonlinear, the wave-induced boundary layer flow is likely to be turbulent over portions of its development. As a result, aNavier-StokessolverwithaturbulencemodelforahighReynoldsnumber, wall-bounded flow is necessary. In the following sections the fundamentals of the turbulencemodelingemployedinthisworkarepresented. Inparticular, wepresent specifics of the k−ω closure model. 70 2.2.1 Reynolds-averaged equation Weconsidertheplanemotionofaviscousflowwhichisgovernedbythedimensional Navier-Stokes equations for an incompressible fluid: ∂u i ∂x i =0, (2.30a) ∂u i ∂t +u j ∂u i ∂x j =− 1 ρ ∂p ∂x i +ν ∂ 2 u i ∂x j ∂x j , i = 1,2. (2.30b) In these equations u i represents the velocity component in the i direction, p is the thermodynamic pressure, ρ is the (constant) density of the fluid, and ν is the kinematic viscosity of the fluid. It is commonly believed that the Navier- Stokesequations(2.30)containalloffullyturbulentdynamicsoffluidmotion[Fri95, p.1]. Since turbulence consists of statistical fluctuations, the velocity and pressure fields in equation (2.30) are nominally decomposed into a mean component plus a fluctuating component. That is, the velocity and pressure fields are represented as u i =U i +u 0 i , p=P +p 0 , where U i =u i , P =p. (2.31) The overline indicates some sort of averaging operation. Theoretically, ensemble averages (i.e., averaging over a large number of identical experiments) are most 71 often implied. However, in typical engineering applications, a time averaging is frequently used; namely, U i ≡u i = 1 τ Z to+τ t o u i (t)dt, u 0 i =0. (2.32) The averaging time τ must be chosen to be much larger than any time scale for the fluctuating portion of the dynamics. The decompositions (2.31) are then sub- stituted into equations (2.30), and the equations are subjected to the time average operation, leading to equations for the mean-flow. The equations for the mean flow quantities, i.e., U i , P, etc., take the form ∂U i ∂x i =0, (2.33a) ∂U i ∂t +U j ∂U i ∂x j =− 1 ρ ∂P ∂x i + 1 ρ ∂ ∂x j ³ 2μS ji −ρu 0 j u 0 i ´ . (2.33b) The set of equations (2.33) are called the Reynolds-averaged Navier-Stokes equa- tions, or the RANs model. Detailed derivations can be found in [TL72, Sec. 2.1] or [Wil06, Sec. 2.3]. The quantity represented by S ij is mean stress tensor, defined as S ij = 1 2 µ ∂U i ∂x j + ∂U j ∂x i ¶ . (2.34) The rightmost term in equation (2.33b), ρu 0 j u 0 i , is known as Reynolds stress tensor and is denote herein by ρτ ij . It is quite obvious that, unless the equations for the 72 Reynolds stress terms can be obtained, or their correlations with the mean-flow quantities can be defined, the number of unknowns is greater than the number of equations. This is the classical closure problem of turbulence, and is the focus of different models of turbulence. A specific model for the Reynolds stress terms, the one employed in this study, is discussed in the next section. 2.2.2 k−ω model Oneobservesfirstthatthetworightmosttermsinequation(2.33b)canbecombined if the entire stress tensor can be expressed as 2μS ji −ρu 0 j u 0 i =2S ji ρν t , (2.35) where ν t =μ t /ρ is called eddy viscosity. The concept of eddy viscosity is based on the assumption that the mechanism of turbulent diffusion is analogous to that of moleculardiffusion. Thisassumptionitselfhasmoretodowithdimensionalanalysis than fundamental physics. Since the Reynolds stress tensor varies in both space and time, it is clear that the eddy viscosity is also, in general, a function of both time and space. In order that the RANs set of equations form a closed system for the mean flow state, it is necessary that ν t be expressed in terms of other dynamic variables. Wechoosehereintoemploythek−ω modeldevelopedbyWilcox[Wil06] to resolve the closure problem. In the k−ω model, the eddy viscosity is calculated as the ratio of the turbulent kinetic energy k to the specific energy dissipation rate 73 ω. Twoadditionalequations,(2.36b)and(2.36c),areintroducedtocalculate k and ω in terms of the mean flow state. The entire k−ω model reprinted from [Wil06, Sec. 4.3.1] is as follows: ν t = k ˜ ω , ˜ ω =max ( ω, C lim s 2S ij S ij β ∗ ) , C lim = 7 8 . (2.36a) ∂k ∂t +U j ∂k ∂x j =τ ij ∂U i ∂x j −β ∗ kω+ ∂ ∂x j ·µ ν +σ ∗ k ω ¶ ∂k ∂x j ¸ . (2.36b) ∂ω ∂t +U j ∂ω ∂x j =α ω k τ ij ∂U i ∂x j −βω 2 + σ d ω ∂k ∂x j ∂ω ∂x j + ∂ ∂x j ·µ ν +σ k ω ¶ ∂ω ∂x j ¸ . (2.36c) α = 13 25 , β =β oo f β , β ∗ = 9 100 , σ = 1 2 , σ ∗ = 3 5 , σ do = 1 8 . (2.36d) σ d = 8 > > > < > > > : 0, ∂k ∂x j ∂ω ∂x j ≤0 σ do , ∂k ∂x j ∂ω ∂x j >0. (2.36e) β oo =0.0708, f β = 1+85χ ω 1+100χ ω , χ ω ≡ ¯ ¯ ¯ ¯ Ω ij Ω jk S ki (β ∗ k) 3 ¯ ¯ ¯ ¯ . (2.36f) Ω ij = 1 2 µ ∂U i ∂x j − ∂U j ∂x i ¶ , S ij = 1 2 µ ∂U i ∂x j + ∂U j ∂x i ¶ . (2.36g) ThismodelwasfirstproposedbyKolmogorov[Kol42]andlatermodifiedbyothers (e.g., Saffman [Saf70], Launder and Spalding [EB72], and Wilcox [Wil06]). Devel- opers of the model differ somewhat in the addition of selected terms and different calculations of the closure coefficients (e.g., α, β ∗ , etc.). Equation (2.36c) for ω 74 has a form similar to that for k. The expression for ν t and the equation for ω are formulatedbasedondimensionalanalysis,physicalreasoning,andmeasuredresults [Wil06, Sec. 4.3]. It is worth noting that when Kolmogorov first postulated the equation for the dissipation rate, (2.36c), he did not include any production term; that is, the first term on the right hand side of (2.36c). This term is vital as it pro- vides the means through which the turbulence extracts energy from the mean flow. He may have reasoned that, since dissipation occurs mainly at the smallest scales, the dissipation rate ω shouldn’t be related to the large-scale mean flow. However, it is commonly accepted today that even though the dissipative process is mostly active at small scales, the rate of dissipation is determined by the interaction be- tween the large eddies and mean flow (see [Sch00, Sec. 16.5.2]). Therefore, the equation for the dissipation rate ω must include a production term. The closure coefficients are introduced to replace the various high-order moments (e.g., p 0 u 0 j , u 0 i u 0 j u 0 j , etc.) that arise in the process of deriving equation (2.36). Their values are decided through comparison of simulations with experimental observations, as well as other theoretical calculations and constraints. The particular values must nec- essarily be optimized according to different generic flow applications; for example, boundary-layer flows, jets, wakes, etc. Equation (2.36) is referred to as the Wilcox 2006 k− ω model. Equations (2.33) and (2.36) are coupled and must be solved simultaneously with a proper numerical scheme. A particular numerical scheme for 75 solution of these equations is discussed in the next chapter. For a general review of different turbulent models, see [Wil06] or [Sag01]. 76 y * c k 1 L * w k 2 L * w L * c L * w Figure 2.10: The area of computation is confined by the dash lines and the bottom wall, i.e., y ∗ c ·L ∗ c k 1 = k 2 = 6,y ∗ c >> δ(x ∗ ). The characteristic wave length L ∗ w is computed from equation (2.13), and is shown here as merely an visual estimate. 77 Part II, Chapter 3 Numerical Issues 3.1 Initialization The numerical program, called EDDY2C, used to solve the compressible versions of equations (2.33) and (2.36), was developed originally by Wilcox [Wil06]. The program employs the MacCormack (1985) fully-implicit flux-splitting method with Gauss-Seidellinerelaxation; see[Wil06,Chap. 3]. Thealgorithmpossessessecond- orderaccuracyinbothspaceandtime. ThefundamentalconceptsofMacCormack’s Predictor-Correctorschemecanbefoundin[And94],withfurtherdetailin[Mac82]. What is emphasized here is that, in order to simulate the incompressible flow with a compressible numerical code, some scaling issues must be addressed to simulate an incompressible state relevant to ocean applications. 78 From both physical and theoretical points of view, the flow field ahead of the propagating wave, when viewed in the wave-fixed frame, should consist of a uni- form, horizontal stream moving to the right with the speed C w of the wave. In other words, when the numerically-generated steady-state solution is obtained, the boundary-layer flow velocity ahead of the wave (left) should be zero relative to the wall. However, owing to the compressible-to-incompressible correction made in the current study (see Appendic (B.3)), and with the expressions for the velocity field derived for an incompressible flow, this theoretical condition is not fully satisfied withoutsomecorrectionfor(small)compressibilityeffectsimplicittothenumerical solver. As noted in Appendix (B.3), we fix the Mach number of the ambient flow to be M a = 0.05. Then, as weak compressibility admits upstream propagation of pressure disturbances, a modification of the upstream state is required to realize the desired velocity at the inlet to the computational domain. A simple iterative scheme is implemented to adjust the inlet flow to match the theoretically-derived state corresponding to any particular solitary wave solution forincompressibleflow. Thisisaccomplishedbyuseofiterativeadjustmentscheme specified by u i+1 (−∞,y ∗ )=u i (−∞,y ∗ )−(u i (x ∗ o ,y ∗ )−C ∗ w ), (3.1a) u 0 (−∞,y ∗ )=C ∗ w , (3.1b) 79 where u i (x o ,y) is the computed, stationary velocity field at a streamwise position x o ahead of the wave where the local wave amplitude is less than one percent of the maximum wave amplitude at the wave center. U(−∞) is the speed of the upstream inlet flow needed to yield the correct ambient state ahead of the wave. The iterative scheme (3.1) is applied consecutively until the velocity difference u(x o ,y)−C w caused by the compressible-to-incompressible correction is reduced below a prescribed level. Figure (3.1) shows the ratio of velocity difference to the wave speed C w for the first three iterations: i.e., i=0,1,2. −0.15 −0.1 −0.05 0 0.05 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 (u i (x * o ,y * )−C * w )/C * w y * i=0 i=1 i=2 Figure 3.1: Iterative convergence of inlet flow (see equation (3.1)). h=8, ζ o =3.4. The figure clearly demonstrates that there is a relatively rapid convergence to the required inlet flow state. The overall ratio reduces from 12% to about 2% with merely three iterations. We choose to discontinue the iterative process when the following condition is met: ½ u i (x ∗ o ,y ∗ )−C ∗ w C ∗ w ¾ max ≤2×10 −2 . (3.2) 80 Forallthenumericalcomputationsdescribedhere,threeiterationsprovedsufficient to meet the criterion. 3.2 Grid and Accuracy Tests The numerical accuracy of the simulations is demonstrated with the help of a tech- nique called Richardson Extrapolation (see [Roa98]). Based on this technique, a Grid Convergence Index (GCI) is introduced by Roache [Roa98] to estimate possi- ble computational errors, where GCI is defined as GCI≡F s ² h r p −1 , ² h = φ h −φ rh φ h , (3.3) and F s is recommended to be 1.25 if three different grids are used. The variable φ represents any computed quantity whose error is to be estimated; e.g., δ 1 , C f , etc. If we are using space-doubling grids r =2, and the order of accuracy p can be estimated based on φ and r: p= ln|(φ rh −φ r 2 h )/(φ rh −φ h )| lnr . (3.4) Since EDDY2C is second-order accurate, p can be simply taken as 2. Tables (3.1), (3.2)and(3.3)reportvaluesofGCIcorrespondingtotheskinfrictioncoefficient,the displacement thickness, and the momentum thickness, respectively, for a flow case where the wave-induced boundary layer transitions to turbulent flow considerably 81 before the peak velocity difference across the boundary layer is reached (i.e., h=3, ζ o /ζ oe =0.93). N b showninthetableisthenumberofgridpointsusedtoresolvethe thinattachedboundarylayerinwhichtheviscosityisimportant. Itsvalueisalways more than half of the total grid points used in transverse direction, i.e., N y . The numerical resolution in the performed calculations is 400 points in the streamwise direction and 260 points in the transverse direction, and with 180 points dedicated to resolving the boundary layer. This means that the value of GCI is no more than 5, even for the flow possessing a reasonably high level turbulence (i.e., high eddy viscosity). 82 Table 3.1: GCI Report for C f N x N y N b C f p ² h GCI 128 64 40 2.092588e-4 None None None 256 128 80 2.1350812e-4 None 0.0199021 -0.0274925 512 256 160 2.1310396e-4 3.39421 -0.0018965 -0.0026198 Table 3.2: GCI Report for δ 1 N x N y N b δ 1 p ² h GCI 128 64 40 3.1755038e-3 None None None 256 128 80 3.3779012e-3 None 0.0599180 -0.08723208 512 256 160 3.4065199e-3 2.82216014 0.0084011 -0.01223087 Table 3.3: GCI Report for δ 2 N x N y N b δ 2 p ² h GCI 128 64 40 1.9161831e-3 None None None 256 128 80 2.0154557e-3 None 0.0492556 -0.07621129 512 256 160 2.0345280e-3 2.3799187 -0.0018965 -0.01450447 3.3 Velocity Profiles and Eddy Viscosity In this section some boundary-layer velocity profiles and the variation of eddy viscosity will be presented to reveal the fundamental flow features obtained by use of the k−ω model and EDDY2C. The kinematic viscosity ν is chosen to be 10 −6 (m 2 /sec). In other words, we define a perfect gas with properties that match the primary state and constitutive properties of liquid water. The values of other required parameters can be found in Appendix (B.1). 83 3.3.1 Sample Velocity Profiles The downstream evolution of the mean velocity profiles corresponding to the wave profile shown in the left panel in Figure (2.2) is shown in Figure (3.2). Panels (a) and (b) in Figure (3.2) show profiles corresponding to the upstream portion of the solitary wave when the boundary layer flow develops under a favorable pressure gradient. The wave-induced boundary layer experiences maximum driving stress (i.e., velocity difference) near the center of the wave, at which location the pressure gradient vanishes. The velocity profile at the wave center is shown in panel (c) in Figure (3.2). Aft of the wave center the driving stress (velocity difference) dimin- ishes and the pressure gradient is adverse. As the boundary layer develops further in the presence of the adverse pressure gradient, it eventually detaches from the wall surface and a separation bubble is likely to form. A velocity profile following boundary layer separation is shown in panel (d) of Figure (3.2). Asimilarsetofmeanvelocityprofilesformedunderanothermuchmoreenergetic solitary wave forcing is shown in Figure (3.3). Owing to the larger wave amplitude (i.e., ζ o /ζ oe =0.99), the favorable and adverse pressure gradients impressed by the externalinviscidflowaresomuchlarger, asistheimpresseddrivingstress(velocity difference), that the boundary layer flow undergoes transition. This is evident by the nearly linear profile variation of the mean velocity outside the thin sub-layer regionadjacenttothewall,acharacteristicoffullyturbulentboundarylayers. Such 84 a characteristic turbulent profile is seen in panel (b) of Figure (3.3) corresponding to a streamwise location very near the wave center. 0 0.005 0.01 0.015 −5 0 5 10 x 10 −3 u/u ∞ y * (a) 0 0.01 0.02 0.03 0.04 −5 0 5 10 x 10 −3 u/u ∞ y * (b) 0 0.02 0.04 0.06 −5 0 5 10 x 10 −3 u/u ∞ y * (c) −0.01 0 0.01 0.02 0.03 −5 0 5 10 x 10 −3 u/u ∞ y * (d) Figure 3.2: (a):◦ ; (b):Δ ; (c):¦ ; (d):¤ (See figure (2.2)). 3.3.2 Eddy Viscosity In turbulence modeling the strength of the turbulence is often defined in terms of an eddy viscosity ν t which measures the strength of the diffusion arising from the turbulent motion. Different RANs models employ different assumptions as to the 85 0 0.05 0.1 0.15 0.2 −0.01 0 0.01 0.02 0.03 u/u ∞ y * (a) 0 0.2 0.4 0.6 0.8 −0.01 0 0.01 0.02 0.03 u/u ∞ y * (b) 0 0.2 0.4 0.6 0.8 −0.01 0 0.01 0.02 0.03 u/u ∞ y * (c) −0.2 0 0.2 0.4 0.6 −0.01 0 0.01 0.02 0.03 u/u ∞ y * (d) Figure 3.3: (a):◦ ; (b):Δ ; (c):¦ ; (d):¤ (See figure (2.2)). definition and calculation of ν t . The definition of ν t in the k−ω model is given in (2.2.2). When the flow is purely laminar, the ratio of the eddy viscosity to the molecular kinematic viscosity is very small (i.e., ν t /ν << 1). In a fully-developed turbulent flow the diffusion of momentum occurs predominately by the turbulent eddy motions. As a consequence, ν t is usually several orders of magnitude larger than the molecular viscosity ν. In our study the highest ratio of ν t /ν, which occurs for wave amplitudes near the critical value (i.e., ζ o /ζ oe ∼1), is about 20. 86 The numerical value of the eddy viscosity is determined purely by the turbulence modeling. Figure(3.4)showsthevariationofthemaximumviscosityratio(ν t /ν) max with respect to ζ o , where the maximum corresponds to the peak value over the en- tire streamwise and transverse development of the boundary layer. It shows that transition occurs when the amplitude ratio is between 0.3 to 0.4. As seen in Figure (3.4), the maximum eddy viscosity becomes (essentially) independent of the depth ratio h for h& 5. The highest viscosity ratio for h = 6 is about twice as large as that for h=3 (not visible from Figure (3.4)). A comparison of a typical transverse variation across the boundary layer of the viscosity ratio ν t /ν at a streamwise lo- cation at the wave center for different amplitude ratios at h=8 is shown in Figure (3.5). The viscosity ratio rises rapidly outside the sublayer, and remains at an elevated level across the major extent of the boundary layer. 0 1 2 3 4 10 −1 10 0 10 1 ζ o (ν t /ν) max 0 0.2 0.4 0.6 0.8 1 10 −1 10 0 10 1 ζ o /ζ oe h=2 h=3 h=6 h=8 Figure3.4: Thevariationofmaximumviscosityratiowithrespecttoζ o . Re w =10 8 . 87 10 −8 10 −6 10 −4 10 −2 10 0 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 y ν t /ν 2.0/3.42 2.5/3.42 3.0/3.42 3.4/3.42 Figure 3.5: transverse variation of viscosity ratio at wave center. h = 8 ; ζ o /ζ oe = 2.0/3.42,2.5/3.42,3.0/3.42 and 3.4/3.42. Re w =10 8 . We choose to define an average across the thickness of the boundary layer of the viscosity ratio (ν t /ν) a (x) as follows: ³ ν t ν ´ ave = 1 y ∗ c Z y ∗ c y ∗ =0 ν t ν dy ∗ (3.5) The position y c is the transverse location where the value of (ν t /ν) is negligible compared to its maximum value at the same streamwise location. This function, which varies in the streamwise direction, indicates in general how the turbulence develops in the streamwise direction. Figure (3.6) shows this variation for the solitary wave profiles given in Figure (2.2). In a purely laminar flow (left panel in Figure (3.6)), ν t is negligible compared to ν. In a turbulent boundary-layer 88 (right panel in Figure (3.6)) the computed average viscosity ratio clearly shows that transition occurs a short distance ahead of the peak velocity difference for the wave. After the flow passes the peak of the wave, the strength of the Reynolds stress begins to diminish as the boundary layer develops further into the region of adverse pressure gradient. As the peak wave amplitude increases, and so also the velocity difference driving the boundary layer, transition moves upstream, even into the region of strongly favorable pressure gradient. An example is shown in the second panel of Figure (3.6). −100 −80 −60 −40 −20 0 0 2 4 6 8 10 12 14 x * −x * c (ν t /ν) a *10 7 ζ o /ζ oe =0.5/3.42 −100 −50 0 50 10 −2 10 −1 10 0 10 1 x * −x * c (ν t /ν) a ζ o /ζ oe =3.4/3.42 Figure 3.6: The streamwise variation of viscosity ratio corresponding to the waves in figure (2.2). h=8 ; Re w =10 8 . With values for the viscosity ratio as computed using the present numerical model, we note the following characteristic features of the development of turbulence in the wave-induced boundary layer. The boundary layer flow is always turbulent over some streamwise extent of its development for all four depth ratios exam- ined, h = 2,3,6,8, so long as the amplitude ratio ζ o /ζ oe exceeds a value of about 89 0.3∼0.4. This is seen, for example, in Figure (3.4). In general, as the depth ratio h is held fixed, the extentof the turbulent region, along with the level of the turbu- lence intensity, increases with increasing wave amplitude. When we compare flows atdifferentdepthratios,wefindthatitisadvantageoustocorrelateboundarylayer properties with the amplitude ratio ζ o /ζ oe , and not with the amplitude ζ o alone. Scaling with the amplitude ratio accounts for the depth ratio effect since the value oftheextremewaveamplitudeζ oe variesalmostlinearlywiththedepthratio h(see the first panel of Figure (2.3)). 90 Part II, Chapter 4 Boundary layer Properties Usingtheturbulencek−ωmodel,boundarylayervelocityprofiles,variousboundary- layerthicknesses,localskinfrictiondistributions,totalfrictionaldrag,andthepoint ofboundarylayerseparationcanbecomputed. Resultsarepresentedforfourdiffer- ent depth ratios, h 2 /h 1 ≡h=2,3,6,8, for comparison, and to obtain a reasonably comprehensive estimate of the main characteristics of the wave-induced boundary layer under long internal waves in shallow seas. As noted earlier, a given solitary wave is labeled by its amplitude ζ o , and there exists an unique value of the extreme amplitudeζ oe for each depth ratio. The correlation between the extreme amplitude ζ oe and the depth ratio h is (nearly) linear, as shown previously in Figure (2.3). Hence, the boundary layer characteristics depend on the following set of param- eters: the wave amplitude ζ o , the depth ratio h, and the Reynolds number. As noted earlier in section (2.1.4), several different Reynolds numbers can be defined depending on the choice of velocity scale and length scale used in its definition. Correlations will be developed using two different Reynolds numbers, one based on 91 ambient wave properties and another based on boundary layer properties. These distinctions will be emphasized later in this section. This chapter is arranged as follows. In Section (4.1), the variation of the skin frictioncoefficientwithrespecttothewaveamplitudeandthewaveReynoldsnum- ber are presented. This is followed in Section (4.2) by an examination of the different measures of the boundary layer thickness. Lastly, in Section (4.3), aspects of boundary-layer separation and the effects of bottom roughness are presented. 4.1 Skin friction coefficient The skin friction, or wall shear stress, is defined as τ o (x)=μ ∂u ∂y | y=0 , (4.1) where μ is the dynamic viscosity and we have positioned the origin of the vertical coordinate at the bottom boundary located at y = 0. A dimensionless measure of the wall shear stress is obtained by defining the local skin friction coefficient, C f , and the total skin friction coefficient, C f,t . These are defined as C f (x)= τ o (x) q ∞ , C f,t = Z x s xo C f (x)dx, q ∞ = 1 2 ρ ∞ C 2 w . (4.2) 92 The local stress is scaled by the dynamic pressure q ∞ of the freestream state. The total frictional drag coefficient is integrated from an upstream position x o where C f (x o ) < 10 −5 to the downstream location of the separation point x s where the local skin friction coefficient vanishes. In the interval x o < x < x s , the local skin friction coefficient increases monotonically to a maximum value C f,m , and then de- creasesmonotonicallyuntilitbecomesnegativedownstreamoftheseparationpoint denoted by x s . Typical profiles of C f (x) for several solitary waves at two depth ratios are shown in Figure (4.1). As expected intuitively, waves with larger amplitudes experience greater local frictional drag induced by the no-slip boundary condition on the bottom wall. The insert in the upper panels of Figure (4.1) give the wave amplitude ratio ζ o /ζ oe for each profile. It is clear that the local stress varies dramatically as the wave am- plitude increases toward its extreme value. For example, C f,m increases by nearly 50% when the amplitude ratio increases by only about 14%. This implies that the correlation between the energy of the wave (∼ζ 2 ) and the induced shear stress on the wall surface is decidedly nonlinear. Severalcompoundinginfluencescontributetothisstrongincreaseinmaximumskin friction. First, as thewaveamplitudeincreasesthe driving stress acrossthe bound- ary layer increases, as seen Figure (2.6). Second, as the wave amplitude increases 93 the characteristic width of the wave L w first decreases, but then increases as the extreme amplitude is approached. This means that the development length of the boundarylayerinitiallydecreasesasthedrivingstressincreases, andthenincreases as the velocity difference continues to increase. These coupled effects contribute to a quite nonlinear dependence of C f,m on wave amplitude. Next, the variation of both C f,t and C f,m with respect to the Reynolds number is presented in Figure (4.2) for a fixed depth ratio h=8. The boundary-layer flow correspondingtolowamplitudecasesshownintheuppertwopanelsinFigure(4.2) is laminar. In the bottom two panels of Figure (4.2) the boundary layer clearly un- dergoestransitionwhenthewaveReynoldsnumberisabout7×10 6 . BothC f,m and C f,t decay continuously with increasing Reynolds number, as expected, when the boundary layer is laminar over its entire extent (upper two panels). The variation with Reynolds number for the two amplitude cases shown in Figure (4.2) has been replotted on log-log scales to expose the functional dependence. As evident in Figure (4.3), the variation for the laminar (low amplitude) case followsastraightlinewithnegativeslope. Analysisofthedatashowsthattheslope is−0.497 ∼ =−0.5,consistentwiththeknownvaraitionfromlaminarboundarylayer theory. Thus,thefunctionaldependenceinthelaminarregionsatisfiestherelations 94 C f,m = K m (h,ζ o /ζ oe ) Re 1/2 w , C f,t = K t (h,ζ o /ζ oe ) Re 1/2 w . (4.3) What remains is to discover the functions K m (h,ζ o /ζ oe ) and K t (h,ζ o /ζ oe ) over the rangeofparameterswhentheboundarylayerremainslaminar. Inasimilarmanner, the slope of decay when the boundary layer is turbulent is -0.65. The implies the functional dependence when the boundary layer is turbulent sastifies the relations C f,m = d K m (h,ζ o /ζ oe ) Re 0.65 w , C f,t = c K t (h,ζ o /ζ oe ) Re 0.65 w . (4.4) Again, a set of simulations could be performed to reveal the functional variation of d K m and c K t . Possiblereasonswhytheboundary-layerflowremainslaminarevenforhighReynolds numbersarethefollowing. First,thevalueofthewaveReynoldsnumberRe w =10 8 may exceed the parameter range for which the numerical code is designed. Second, the Reynolds number required for transition to occur for a low amplitude ratio (e.g., ζ o /ζ oe ∼ 0.14) may be as high as Re w = 10 9 , or larger. What this implies is that a more relevant Reynolds number needs to be defined for the wave-induced boundary layer, and that any comparison with other experimental or numerical data of nominal boundary layer flows is difficult because of the strongly varying pressure field impressed by the wave, first favorable and then adverse. To this end, 95 we plot in Figure (4.4) the skin friction data as a function of the Reynolds number basedonthepeak velocitydifference acrossthe boundary layerusing the boundary layer momentum thickness as the length scale. The variation of C f,m and C f,t with respect to this Reynolds number is shown for three different wave amplitudes. One notes from Figure (4.4) that the local Reynolds number based on the momentum thickness and local velocity difference is a good indicator of when transition to a turbulent boundary layer flow occurs. In this wave-induced boundary layer, transition occurs at approximately the same value of the local Reynolds number (Re δ 2 ,ΔU c ≈10 3 ) for all amplitude ratios. The variation of the skin friction with respect to ζ o /ζ oe is presented in Figure (4.5). Overall, the skin friction is proportional to the the amplitude ratio. When the wave amplitude is small, the correlation between the wave amplitude ratio, ζ o /ζ oe , and the variation of total or maximum skin friction coefficients, C f,t or C f,m are almost linear. This linearity does not exist when the flow becomes turbulent at higher amplitude ratios. The turbulence k−ω model [Wil06] utilizes the concept of roughness height, h r in the prescription of the wall boundary condition on the specific dissipation rate ω. The roughness height affects both the Reynolds number for transition to 96 turbulence and the turbulence intensity. Figure (4.6) depicts the variation of the skinfrictioncoefficientwithrespecttothelocalReynoldsnumberforthreedifferent scaled roughness heights (i.e., h r =h ∗ r /δ ∗ 2 ); h r =0.308,1.206,2.914. The roughness heighth r =0.308isusedthroughoutthisstudyasthebase-line(smoothwall)value. The k−ω model predicts earlier transition and enhanced skin friction coefficient as the roughness height increases. 4.2 Boundary Layer Thickness The displacement thicknes δ 1 (x) and the momentum thickness δ 2 (x) were defined earlier in equations (2.26) and (2.27). The variations of δ 1 and δ 2 with respect to the wave amplitude at the wave center are shown in Figure (4.7). Basedontheresultsshowninthisfigure,wecannotethefollowinggeneralchar- acteristics. When the boundary-layer is laminar (i.e., the amplitude ratio is small), the displacement and momentum thicknesses decrease with increasing amplitude ratio because the impressed favorable pressure gradient is increasing. The decrease in these thicknesses with increasing favorable pressure gradient on the front of the wave for low to moderate amplitude ratio is consistent with theory and experiment of nominal boundary layer flows. The increase in these thicknesses with ampli- tuderatioathighervalues, eventhoughthepressuregradientbecomesincreasingly 97 favorable, comes about because the boundary layer development length is also in- creasing (see Figure (2.4)). The increase in these thicknesses for amplitude ratios with ζ o /ζ oe >0.5 is the consequence of a more extensive upstream region of a more fully developed turbulent boundary layer. Variations of δ 1 (x) and δ 2 (x) with respect to the wave Reynolds number are shown in Figure (4.8). For laminar flow (i.e., the top two panels for a wave with low amplitude ratio), the thicknesses vary (essentially) as the inverse square-root of the Reynolds number. For a larger wave amplitude case (i.e., the bottom two panelsinFigure(4.8)),transitionoccursataboutRe w ∼3×10 6 ,andtheboundary layer is essentially turbulent when Re w ∼7×10 6 . The variation with respect to wave amplitude ζ o of the shape factor H =δ 1 /δ 2 is presented in Figure (4.9), where results are plotted for the shape factor at the wave-center location. Values of the shape factor for some typical boundary layer flows (i.e., boundary layers with fixed, constant pressure gradient) are useful as a context for interpreting computed results and are summarized in Table (4.1). At fixed pressure gradient, and fixed flow state (i.e., laminar or turbulent), the shape factor is a constant independent of Reynolds number. That is, the boundary layer develops in a similarity form. In the present case, as already noted, the pressure 98 gradient is variable and the boundary layer development length also varies. Fur- thermore, there is no distinct origin for the boundary layer, starting from a free origin where the driving velocity difference increases from zero in the downstream direction. Hence, there is no typical boundary layer for direct comparison. Flow state Laminar Turbulent flat plate H =2.6 H =1.3 ∂P/∂x< 0 H <2.6 H <1.3 ∂P/∂x> 0 H >2.6 H >1.3 separation point H ∼ = 3.5 H ∼ = 2.4 Table 4.1: Shape factor H =δ 1 /δ 2 for typical boundary layer. The shape factor for low amplitude ratios increases with amplitude, but is always below the flat plate (Blasius; zero pressure gradient) value of H = 2.6. Since the computedresultssatisfyH <2.25forsmallamplitudes, theyareatleastconsistent with known data for typical boundary layers developing under a favorable pressure gradient. In addition, the shape factor is seen to decrease significantly for higher amplitude ratios. This rapid decrease in H is also consistent with the fact that the boundary layer has transitioned to a turbulent state. In Figure (4.10) the shape factor at three streamwise positions in the boundary layer are tracked as a function of wave amplitude for a fixed depth ratio h = 8 and a fixed Reynolds number Re w = 10 8 . The shape factor was evaluated at the inflection point on the leading face (L.I.) and the rear face (R.I.) of the wave profile. It was also computed at the wave center (C.). The shape factor at the 99 rear inflection point is consistently above the flat plate (zero pressure gradient) value H =2.6 for a laminar boundary layer, and considerably above the separation value H ∼ = 2.4 for a turbulent boundary layer. Since the rear inflection point is the location of maximum adverse pressure gradient, these values conform to expected trends. Although the boundary layer at the left (front) inflection point is (almost always) laminar, the strong favorable pressure gradient gives rise to a somewhat surprisingly low shape factor. Of course, this inflection point is near the origin of the boundary layer where the driving velocity difference builds slowly from a zero value. Hence, inferences drawn from nominal boundary layer behavior may not be relevant. That the shape factor at the wave center decreases slightly with wave amplitude is likely a manifestation of the upstream movement of transition and the lengthening of the region of turbulent flow. Even so, the value of 1.8<H c <1.9 is considerably above the nominal value of H ∼ = 1.3 for turbulent boundary layers on a flat plate (zero pressure gradient). 4.3 Separation In all the cases that with amplitude ratios in the range 0.03≤ ζ o /ζ oe ≤ 0.95, and wave Reynolds numbers in the range 10 6 ≤ Re w ≤ 10 9 , boundary layer separation 100 alwaysoccurs. Thepointofseparation x sep isdefinedasthepointofvanishingskin friction under the aft portion of the wave. That is, µ ∂u ∂y | y=0 ¶ x sep =0. (4.5) In some cases, the flow separates upstream of the inflection point along the rear half of the wave profile, and in some cases the separation point occurs behind the rear inflection point. The variation of the point of separation with respect to the wave amplitude for various depth ratios is shown in Figure (4.11). The normalized separationpointx s employedinFigure(4.11)isdefinedrelativetothewavecenter. That is, x s =x sep −x c (4.6) where x sep is the point of separation defined by (4.5). It is interesting to note that the different curves for various depth ratios coalesce into somewhat of a universal variation when plotted as a function of the amplitude ratio. The normalized distance between the separation point and the rear inflection point is defined as D SR = x sep −x R.I. L w . (4.7) ThevariationofD SR withrespecttotheamplituderatioisplottedinFigure(4.12) for h = 3,8 at fixed Reynolds number Re w = 10 8 . It is clear that the point of 101 separation moves upstream of the rear inflection point as the amplitude increases. Thisimpliesthatthelengthofanyseparationbubbleaftofthewavecenterincreases as the wave amplitude increases. 102 −1 0 1 0 0.5 1 1.5 2 2.5 3 3.5 h=8 0.99 0.87 0.61 0.234 −1 0 1 0 0.2 0.4 0.6 0.8 ζ h=3 0.81 0.71 0.51 0.30 −1 0 1 −5 0 5 10 15 20 x 10 −5 (x * −x * c )/L * w C f −1 0 1 −6 −4 −2 0 2 4 6 8 x 10 −4 (x * −x * c )/L * w Figure 4.1: top: wave profiles ; bottom: the corresponding C f profiles. Re w =10 8 ; numbers in the legends are the amplitude ratio, ζ o /ζ oe . 103 10 6 10 7 10 8 10 9 0 0.2 0.4 0.6 0.8 1 1.2 x 10 −4 C f,m ζ o /ζ oe = 0.5/3.42 10 6 10 7 10 8 10 9 0 2 4 6 8 x 10 −3 C f,t ζ o /ζ oe = 0.5/3.42 10 6 10 7 10 8 10 9 0 0.5 1 1.5 2 x 10 −3 C f,m ζ o /ζ oe = 3.0/3.42 10 6 10 7 10 8 10 9 0 0.02 0.04 0.06 0.08 0.1 C f,t ζ o /ζ oe = 3.0/3.42 Figure 4.2: The variation of the maximum and total skin friction coefficients with respect to the wave Reynolds number. 10 6 10 7 10 8 10 9 10 −6 10 −5 10 −4 10 −3 C f,m 10 6 10 7 10 8 10 9 10 −4 10 −3 10 −2 10 −1 C f,t Figure 4.3: Replot of skin friction variation in log-log scale: − : ζ o /ζ oe = 0.5/3.42; −.−:ζ o /ζ oe =2.0/3.42.. 104 10 3 10 4 0 0.5 1 1.5 2 x 10 −3 Re δ 2 ,Δ U c ζ o = 2.0, 2.5, 3.0 C f,m 10 2 10 3 10 4 10 5 0 0.02 0.04 0.06 0.08 0.1 Re δ 2 ,Δ U c ζ o = 2.0, 2.5, 3.0 C f,t 2.0 2.5 3.0 Figure4.4: ThevariationofskinfrictioncoefficientagainstlocalReynoldsnumber, Re δ 2 ,Uc . h=8. ζ oe =3.42. 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 x 10 −4 ζ o /ζ oe C f,m 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 ζ o /ζ oe C f,t h=2 h=3 h=6 h=8 Figure 4.5: The variation of maximum and total skin friction coefficients with respect to amplitude ratio. Re w =10 8 . 105 10 3 10 4 0 2 4 6 8 x 10 −4 Re δ 2 ,Δ U c C f,m h r = h * r /δ * 2 10 3 10 4 0 0.01 0.02 0.03 0.04 Re δ 2 ,Δ U c C f,t h r = h * r /δ * 2 2.914 1.206 0.308 Figure 4.6: The variation of skin friction coefficients with respect to local Reynolds number for three different wall roughnesses. δ ∗ 2 . : momentum thickness. 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 x 10 −3 ζ o /ζ oe δ 1 0 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 x 10 −4 ζ o /ζ oe δ 2 h=2 h=3 h=6 h=8 Figure 4.7: The variation of displacement and momentum thicknesses at the center of the wave with respect to the amplitude ratio. δ 1 and δ 2 are scaled by h 1 . 106 10 6 10 7 10 8 10 9 0 0.002 0.004 0.006 0.008 0.01 δ 1 ζ o = 0.5 10 6 10 7 10 8 10 9 0 1 2 3 4 5 x 10 −3 δ 2 ζ o = 0.5 10 6 10 7 10 8 10 9 0 2 4 6 8 x 10 −3 δ 1 ζ o = 3.0 10 6 10 7 10 8 10 9 0 1 2 3 4 x 10 −3 δ 2 ζ o = 3.0 Figure 4.8: The variation of boundary-layer thicknesses with respect to the wave Reynolds number Re w . h=8. δ 1 and δ 2 are scaled by the upper layer depth h 1 . 0 1 2 3 4 1.8 1.9 2 2.1 2.2 2.3 ζ o δ 1 /δ 2 0 0.2 0.4 0.6 0.8 1 1.8 1.9 2 2.1 2.2 2.3 ζ o /ζ oe δ 1 /δ 2 h=3 h=6 h=8 Figure 4.9: The variation of shape factor at the wave center with respect to the wave amplitude ζ o and the amplitude ratio ζ o /ζ oe . Re w =10 8 . 107 0 0.2 0.4 0.6 0.8 1 1.6 1.8 2 2.2 2.4 2.6 2.8 3 ζ o /ζ oe δ 1 /δ 2 L.I. C. R.I. Figure 4.10: The variation of shape factor as a function of amplitude ratio at the leading inflection point (L.I.), wave center (C.) and rear inflection point (R.I.). Re w =10 8 . h=8. 0 1 2 3 4 0.25 0.3 0.35 0.4 0.45 0.5 ζ o x s 0 0.2 0.4 0.6 0.8 1 0.25 0.3 0.35 0.4 0.45 0.5 ζ o /ζ oe x s h=3 h=6 h=8 Figure 4.11: The variation of D SR (defined in equation (4.7)) with respect to the amplitude ratio. Re w =10 8 . 108 0 0.2 0.4 0.6 0.8 1 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 ζ 0 /ζ oe D SR h=3 h=8 Figure 4.12: D SR :(x sep −x R.I. )/L w . Re w =10 8 . 109 Part II, Chapter 5 Concluding Perspectives Thewave-inducedboundarylayerunderafamilyofinternalsolitarywaveshasbeen studied through use of a two-dimensional, RANs code employing the k−ω turbu- lence model. The principle motivation is to compute the frictional drag, both in its spatial distribution and its integrated value, for a range of relevant environmental conditions. Models for the evolution and propagation of long, nonlinear internal waves in shallow seas has matured considerably in the last decade, but charac- teristics of the frictionally-induced drag acting on these waves as they propagate shoreward on the continental shelf is lacking. Long internal waves with quite large amplitude are ubiquitous features on continental shelves, and they are believed to have important implications for mixing and transport in coastal regions. Hence, understanding the decay of these waves, as well as their generation and propaga- tion, is an essential element in characterizing the strength of these waves in the nearshore region. 110 The data presented here can be used, together with evolution models such as (2.4) for the wave dynamics, to derive estimates of the frictional decay of the wave energy as it evolves while propagating in shallow waters. The derivation of such estimatesisclearlyanimportantnextstepbeyondthepresentstudy. Realisticesti- matesofthedecayrateincorporatingeffectsoftransitionalandturbulentboundary layers, and the effect of benthic roughness, will provide a significant advance in the application of long wave models to realistic wave conditions. With quantitative drag data as a function of key environmental parameters, the temporal decay of internal solitary waves of depression is now essentially in hand for the first time. Afurtherextensionofthisworkwouldbetoemploythepresentwavemodeland the RANs code to develop an equivalent of the well-known Moody diagram (e.g., [Wil10], [Whi03]) for the present family of solitary waves. Such a result would provide a valuable tool for both modelers and field investigators. The steps toward producing such a diagram are clearly set forth in this study. What remains is an extensivesetofsimulationscoveringthefullrangeoftherelevantparameters: depth ratio, amplitude ratio, and Reynolds number. The generation of the needed data is seen as much more of a production study in contrast to a research investigation. The theoretical and numerical development discussed here clearly points the way to cover the parameter space for the necessary production runs. 111 References [And94] John David Anderson. Computational fluid dynamics : the basics with applications. McGraw-Hill, first edition, 1994. [Bat67] G. K. Batchelor. An Introduction To Fluid Dynamics. Cambridge Uni- versity Press, first edition, 1967. [BDR97] D. Bogucki, T. Dickey, and Larry G. Redekopp. Sediment resuspension and mixing by resonantly generated internal solitary waves. Journal of Physical Oceanography, 27:1181–1196, 1997. [Boy01] J.P.Boyd. Chebyshev and Fourier Spectral Methods. Dover,firstedition, 2001. [BR99] Darek J. Bogucki and Larry G. Redekopp. A mechanism for sediment resuspension by internal solitary waves. Geophysical Research Letters, 26:1317–1320, 1999. [CC99] Wooyoung Choi and Roberto Camassa. Fully nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics, 396:1–36, 1999. [CD06] M. Carr and P. A. Davies. The motion of an internal solitary wave of depression over a fixed bottom boundary in a shallow, two-layer fluid. Physics of Fluids, 18(01), 2006. [DDH05] P. J. Diamessis, J. A. Domaradzki, and J. S. Hesthaven. A spectral mul- tidomainpenaltymethodmodelforthesimulationofhighreynoldsnum- ber localized stratified turbulence. Journal of Computational Physics, 202:298–322, 2005. [DFM02] M.O. Deville, P.F. Fischer, and E.H. Mund. High Order Methods for Incompressible Fluid Flow. Cambridge University Press, first edition, 2002. [DLD08] P.J. Diamessis, Y.C. Lin, and J.A. Domaradzki. Effective numerical vis- cosity in spectral multidomain penalty method-based simulations of lo- calized turbulence. Journal of Computational Physics, 227:8145–8164, 2008. 112 [DM86] Grant W. D. and O. S. Madsen. The continental-shelf bottom boundary layer. Annu. Rev. Fluid Mech., 18:265–305, 1986. [DR06] P. J. Diamessis and L. G. Redekopp. Numerical investigation of solitary internal wave-induced global instability in shallow water benthic bound- ary layers. Journal of Phys. Oceanogr., 36:784–812, 2006. [DX03] J.A.DomaradzkiandZ.Xiao. Effectiveeddyviscositiesinimplicitlarge eddy simulations of turbulent flows. Physics of Fluids, 15:3890–3893, 2003. [EB72] Launder B. E. and Spalding D. B. Mathematical Models of Turbulence. Academic Press, London, first edition, 1972. [Fri95] Uriel Frisch. Turbulence, the legacy of A. N. Kolmogorov. Cambridge University Press, 40 West 20th Street, New York, NY 10011-4211, first edition, 1995. [Hes97] J. S. Hesthaven. A stable penalty method for the compressible navier- stokes equations: Ii. one-dimensional domain decomposition scheme. SIAM J. Sci. Comput., 18:658–685, 1997. [JD65] ZabuskyN.J.andKruskalM.D. Interactionsofsolitonsinacollisionless plasma and the recurrence of initial states. Phys. Rev. Lett., 15:240–243, 1965. [jWR00] BangjiWangandLarryG.Redekopp. Longinternalwavesinshearflows: topographic resonance and wave-induced global instability. Dynamics of Atmospheres and Oceans, 33:263–302, 2000. [KC08] Pijush K. Kundu and Ira M. Cohen. Fluid Mechanics. Elsevier Inc., 525 B. Street, Suite 1900. San Diego, California 92101-4495,USA, fourth edition, 2008. [KdV95] D. J. Korteweg and G. de Vries. On the change of form of long waves advancing in a rectangular canal, and on a new type of long standing waves. Philosophical Magazine, 39:422–443, 1895. [KIO91] G. E. Karniadakis, M. Israeli, and S. A. Orszag. High-order splitting methodsfortheincompressiblenavier-stokesequations. Journal of Com- putational Physics, 97:414–443, 1991. [Kol42] A. N. Kolmogorov. Equations of turbulent motion of an incompressible fluid. Izvestia Academy of Sciences, USSR, 1942. 113 [Lam94] Kevin G. Lamb. Numerical experiments of internal wave generation by strong tidal flow across a finite amplitude bank edge. Journal of Geo- physical Research, 99:843–864, 1994. [Les90] M. Lesieur. Turbulence in Fluids. Kluwer Academic, second edition, 1990. [Les96] M Lesieur. New trends in large-eddy simulations of turbulence. Annu. Review Fluid Mech., 28:45–82, 1996. [LPC07] Philip L.-F. Liu, Yong Sung Park, and Edwin A. Cowen. Boundary layer flow and bed shear stress under a solitary wave. Journal of Fluid Mechanics, 574:449–463, 2007. [Mac82] R. W. MacCormack. A numerical method for solving the equations of compressible visous flows. AIAA Journal, 1982. [Mil76] John W. Miles. Korteweg-devries equation modified by viscosity. The Physics of Fluids, 19:1063–1063, 1976. [NW94] Charles A. Nittrouer and L. Donelson Wright. Transport of particles across continental shelves. Reviews of Geophysics, 32:85–113, 1994. [NY91] M.N.NicholsandS.Young. The Amazing LA Environment:A Handbook for change. Living Planet Press, first edition, 1991. [OS70] E. Ott and R. N. Sudan. Damping of solitary waves. The Physics of Fluids, 13:1432–1434, 1970. [OS05] L.OstrovskyandY.A.Stepanyants. Internalsolitonsexistinlaboratory experiments : Comparison with theoretical models. Chaos, 15, 2005. [Pio99] U. Piomelli. Large-eddy simulations:achievements and challenges. Prog. Aero. Sci., 35:335–412, 1999. [Pop00] S. B. Pope. Trubulent Flows. Cambridge University Press, first edition, 2000. [Roa98] P. J. Roache. Verification and Validation in Computational Science and Engineering. Hermosa Publishers, first edition, 1998. [Saf70] P. G. Saffman. A model for inhomogeneous turbulent flow. Proc. R. Soc., Lond., A317, 1970. [Sag01] P. Sagaut. Large Eddy Simulation for Incompressible Flows: An Intro- duction. Springer Verlag, first edition, 2001. 114 [Sch00] Hermann Schlichting. Boundary-Layer Theory. McGraw-Hill, Inc., 1221 Ave.of theAmericas, NewYork, NY10020-1095 USA,8th edition, 2000. [SO98] T. P. Stanton and L. A. Ostrovsky. Observations of highly nonlinear internal solitons over the continental shelf. Geophysical Research Letters, 25:2695–2698, 1998. [Spa09] Philippe R. Spalart. Detached-eddy simulation. Annu. Review Fluid Mech., 41:181–202, 2009. [SR07] T. Sakai and L. G. Redekopp. Models for strongly-nonlinear evolution of long internal waves in a two-layer stratification. Nonlinear Processes in Geophysics, pages 31–47, 2007. [TL72] H. Tennekes and J. L. Lumley. A First Course in Turbulence. The MIT Press, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, first edition, 1972. [Whi03] Frank M. White. Fluid mechanics. McGraw-Hill, 1221 Avenue of the Americas, New York, NY 10020, 5th edition, 2003. [Wil06] David C Wilcox. Turbulence Modeling for CFD. DCW Industries, Inc., 5354 Palm Drive, La Canada, California 91011, third edition, 2006. [Wil10] David C Wilcox. Basic Fluid Mechanics. DCW Industries, Inc., 5354 Palm Drive, La Canada, California 91011, forth edition, 2010. 115 Appendix A Spectral Energy Equation In order to avoid confusion about indices and imaginary unit, i, the Navier-Stokes equation (1.1) are written using different indices, ∂u n ∂x n =0 (A.1a) ∂u n ∂t +u m ∂u n ∂x m =− 1 ρ ∂p ∂x n +ν ∂ 2 u n ∂x m ∂x m (A.1b) The Fourier and Inverse Fourier Tranforms are defined as, ˜ u n (k)≡F{u n }= 1 (2π) 3 Z u n (x)e −ik·x d 3 x (A.2a) u n (x)≡F −1 {u n }= Z ˜ u n (k)e ik·x d 3 k (A.2b) Z d 3 x= ZZZ dx 1 dx 2 dx 3 (A.2c) Under the condition that |u n (x)|→ 0 as x→±∞. The Fourier transform of the derivatives with respect to time and space are as follows. F ½ ∂u n ∂t ¾ = ∂ ∂t F{u n }= ∂ ∂t ˜ u n (k) (A.3) 116 F ½ ∂u n ∂x 2 ¾ = 1 (2π) 3 ZZZ ∂u n ∂x 2 e −ikx d 3 x = 1 (2π) 3 ZZ dx 1 dx 3 e −i(k 1 x 1 +k 3 x 3 ) Z ∂u n ∂x 2 e −ik 2 x 2 dx 2 = 1 (2π) 3 ZZ dx 1 dx 3 e −i(k 1 x 1 +k 3 x 3 ) ½ u n e −ik 2 x 2 ª ∞ −∞ −(−ik 2 ) Z u n e −ik 2 x 2 dx 2 ¾ = ik 2 (2π) 3 ZZZ u n e −i(kx) d 3 x (A.3.1) Therefore, F ½ ∂u n ∂x j ¾ =ik j F{u n }=ik j ˜ u n (k) (A.3.2) F ½ ∂ 2 u n ∂x 2 j ¾ =(ik j ) 2 F{u n }=−k 2 j ˜ u n (k) (A.3.3) F ½ ∂ 2 u n (x) ∂x j ∂x j ¾ =F ½ ∂ 2 u n ∂x 2 1 + ∂ 2 u n ∂x 2 2 + ∂ 2 u n ∂x 2 3 ¾ =−(k 2 1 +k 2 2 +k 2 3 )F{u n }=−k 2 ˜ u n (k) (A.3.4) With these formulae derived, the nonlinear term of the Navier-Stokes equations, (1.1), can be Fourier-transformed as follow, F ½ u m ∂u n ∂x m ¾ =F ½ ∂u n u m ∂x m ¾ =ik m F{u n u m }=ik m Z ˜ u n (p)˜ u m (k−p)d 3 p (A.4) Notice that the convolution theorem is used in the last expression. The pressure term is taken care of by first taking the divergence of equation (A.1b), then apply incompressibility condition to remove the time-derivative term and the dissipation term, the remaining equation is, ∂ 2 ∂x m ∂x n (u m u n )=− 1 ρ ∂ 2 ∂x 2 n p(x) (A.5) Next, apply Fourier transform and use equation (A.3.2) and (A.3.4), one obtains −k n k m Z ˜ u n (p)˜ u m (k−p)d 3 p= k 2 ρ ˜ p(k) (A.6) Finally, the Fourier transform of equation (A.1b) is ∂ ∂t ˜ u n =iP lnm Z ˜ u l (p)˜ u m (k-p)d 3 p−νk 2 ˜ u n (k) (A.7) 117 where P lnm =−k m ½ δ nl − k n k l k 2 ¾ (A.8) The evolution equation for the spectral energy, E(k)= 1 2 u n (k)u ∗ n (k) (A.9) can be obtained by using equation (A.7) as ∂E(k) ∂t =−2νk 2 E(k)+T(k) (A.10) where T(k)= 1 2 Im ½ ˜ u ∗ n (k)P lnm Z ˜ u l (p)˜ u m (k-p)d 3 p ¾ (A.11) Refer to Pope [Pop00] or Lesieur [Les90] for more details. 118 Appendix B B.1 Wave Parameters The derivation of equations (2.16) and (2.17), dimensionless inviscid velocities, are given in this appendix. The starting point is the expression of streamwise velocity from equation (2.14), i.e., u(x ∗ )= C ∗ w h 2 h 2 −η(x ∗ ) , (B.1) which becomes C w c o µ h h−ζ ¶ , (B.2) with the use of equations (2.1), (2.2), (2.3), (2.15) as well as C w = C ∗ w c o . (B.3) Equation (2.16) is readily obtained by dividing equation (B.2) by c oo . The velocity in the tranverse direction V ∞ can be derived in a similar way with the help of continuity equation. We also list here the scaling laws used in this study: ˜ g = ρ 2 −ρ 1 ρ 1 g, c oo = p ˜ gh 1 (B.4a) c 2 o = ˜ gh 1 h 2 h 1 +h 2 , C ∗ w =c oo c o c oo C ∗ w c o C w = C ∗ w c o , (B.4b) ζ = η h 1 , x= x ∗ h 1 , y = y ∗ h 1 . (B.4c) 119 The variables with superscript∗ are all dimensional. The realistic values for ˜ g and c oo in the environmental flows are ˜ g =0.3 (m/sec 2 ) ρ 2 −ρ 1 ρ 1 ∼ = 0.03 g ∼ = 10 (m/sec 2 ) (B.5a) h 1 =5 (m), c oo = p ˜ gh 1 = √ 1.5 (B.5b) B.2 Dimensionless Pressure Gradient In this section, we derive the dimensionless pressure gradient for the external (in- viscid) flow with the assumption that the dependence of the streamwise velocity u(x ∗ ,y ∗ ) on y ∗ is only weakly and can be neglected. With this assumption, the momentum equation for streamwise direction can be expressed as −ρu(x ∗ ) ∂u(x ∗ ) ∂x ∗ = ∂P ∗ ∂x ∗ . (B.6) Substitute equation (2.14) into (B.6), we obtain ∂P ∗ ∂x ∗ =−ρ (C ∗ w h 2 ) 2 (h 2 −η(x ∗ )) 3 dη(x ∗ ) dx ∗ . (B.7) The above equation can be re-casted into its dimensionless form with use of equa- tions (2.2) and (2.3), as well as P = P ∗ ρc 2 o . (B.8) We then arrive at the formula for dimensionless pressure gradient of the external inviscid flow, ∂P ∂x =− (C w h) 2 (h−ζ) 3 dζ dx . (B.9) B.3 CompressibleToIncompressibleCorrections In order to simulate the incompressible flow with a compressible numerical code, there are some engineering tricks or corrections that we have to apply. The nu- merical errors caused by this corrections is not significant. The reduction of the 120 errors by some iterative mechanism will be discussed in the next chapter. The cor- rection starts with a choice of a very low Mach number, e.g., M a =0.05, to realize incompressibility. The Mach number is defined as M a = C ∗ w a ∞ = C ∗ w √ γRT ∞ , γ = C p C v (B.10) Where a ∞ is the speed of sound and γ is the ratio of specific heats. Refer to any classical thermodynamics for the isentropic propagation of sound wave in an ideal gas, e.g., [Bat67] or [KC08]. Then, with the help of equation (B.4b), it is straightforward to show that M a =C w c oo √ γRT ∞ r h 1+h or R = C w 2 c 2 oo M 2 a γT ∞ µ h 1+h ¶ (B.11) The gas constant, R, is then calculated with γ = 1.4, T ∞ = 290 o K. The density of the fluid, ρ 2 , and kinematic viscosity, ν are chosen to be 10 −6 (m 2 /sec) and 1000(kg/m 3 ). In other words, we are somehow working with a kind of gas which possesses properties of liquid water. The values of other used parameters can also be found in Appendix (B.1). 121

Abstract (if available)

Abstract In the first part of this report, the effects of numerical dissipation presented in a turbulence Direct Numerical Simulation (DNS) code called SMPM is investigated through the use of the effective numerical eddy viscosity, a concept pioneered by Domaradzki and Xiao in 2003. This SMPM algorithm, which is inherently free of any intrinsic truncation error, employs three numerical stabilizers, i.e., Fourier filtering, Legendre filtering and Penalty method to ensure the numerical stability in the calculations of turbulent flows with very high Reynolds number, for which the simulation is often under-resolved. Numerical data extracted from the simulation of a turbulent wake flow is employed to quantify the effective numerical eddy viscosity for all three explicit numerical stabilizers. The results shown that the effects of the stabilizers can be significant, if not dominated, compared to the physical molecular viscosity. Away from the vertical subdomain interfaces, all three stabilizers behave expectedly to remove the kinetic energy from the system, i.e., dissipative. At the subdomain interfaces, the Legendre filtering and the Penalty scheme exhibit the nexpected anti-dissipative character, i.e., to increase the energy of the system. Such counter-intuitive behavior is attributed to the discontinous formulism of the Penalty scheme, which is absolutely necessary if multi-domain method is employed. This anti-dissipative behavior diminishes quickly only a few grid points away from the subdomain interfaces.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

Large eddy simulations of laminar separation bubble flows

PDF

Numerical simulations of linearly stratified flow around a sphere

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

A subgrid-scale model for large-eddy simulation based on the physics of interscale energy transfer in turbulence

PDF

Wave induced hydrodynamic complexity and transport in the nearshore

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

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

Open channel flow instabilities: modeling the spatial evolution of roll waves

PDF

RANS simulations for flow-control study synthetic jet cavity on NACA0012 and NACA65(1)412 airfoils.

PDF

Grid-based Vlasov method for kinetic plasma simulations

Asset Metadata

Creator Lin, Yuncheng (author)

Core Title On the simulation of stratified turbulent flows

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace

Publication Date 12/08/2010

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag aliasing effect,Boussinesq approximation,DNS,Fourier filtering,internal solitary waves,k-\omega model,KdV equation,Legendre polynomial,numerical dissipation,OAI-PMH Harvest,penalty method,Poisson equation,SMPM,turbulent boundary-layer flows,under-resolved simulations

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Domaradzki, Julian Andrzej (committee chair), Redekopp, Larry G. (committee chair), Kukavica, Igor (committee member), Wilcox, David C. (committee member)

Creator Email yunchenglin@hotmail.com,yunchenglin@ymail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m3589

Unique identifier UC173842

Identifier etd-Lin-4196 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-429818 (legacy record id),usctheses-m3589 (legacy record id)

Legacy Identifier etd-Lin-4196.pdf

Dmrecord 429818

Document Type Dissertation

Rights Lin, Yuncheng

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu