University of Southern California Dissertations and Theses

Large eddy simulations of laminar separation bubble flows

(USC Thesis Other)

Large eddy simulations of laminar separation bubble flows

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content LARGE EDDY SIMULATIONS OF LAMINAR SEPARATION BUBBLE FLOWS by Francois Cadieux 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) May 2015 Copyright 2015 Francois Cadieux Acknowledgments AheartfeltthankstomysupervisorAndrzejforhisguidance,supportandsense of humour. To Giacomo, thank you for countless in-depth and useful discussions. I am also grateful to those who shared their code with me without which this endeavour might have taken yet another year: Tawan, Brian, Tak, and Peter. I am deeply indebted to Vina, who was always by my side when I needed her most. The support of my parents Johanne and Yves and my sister Genevieve and their eagerness for me to join them in their travels kept me sane and focused. Finally, a special thanks to Dr Philippe Spalart for sharing his DNS data and answering my many questions. This research was supported by the National Science Foundation through grant CBET-1233160. ii Contents Acknowledgments ii List of Tables v List of Figures vi Abstract viii 1 Introduction 1 1.1 Laminar Separation Bubble Flows . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Background 7 2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Reynolds-averaged Navier-Stokes . . . . . . . . . . . . . . . . . . . 9 2.3.1 RANS for Laminar Separation Bubble Flows . . . . . . . . . 11 2.4 Large Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.1 Filtered Navier-Stokes Equations . . . . . . . . . . . . . . . 14 2.4.2 Subgrid-scale Models . . . . . . . . . . . . . . . . . . . . . . 15 2.4.3 LES Results for Laminar Separation Bubble Flows . . . . . 21 2.5 Wall-modeled LES and Hybrid RANS-LES . . . . . . . . . . . . . . 23 3 Methodology 25 3.1 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Flow Specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.1 Center for Turbulence Research Code . . . . . . . . . . . . . 30 3.3.2 Spectral code in vorticity form . . . . . . . . . . . . . . . . . 32 3.3.3 Spectral code in skew-symmetric form . . . . . . . . . . . . 36 3.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 iii 4 Center for Turbulence Research Results 46 4.1 Numerical Dissipation . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.1 Estimating Numerical Dissipation due to Filtering . . . . . . 52 4.1.2 Quantifying Numerical Dissipation . . . . . . . . . . . . . . 54 5 Spectral Results I 59 5.1 LES at 3% of DNS Resolution . . . . . . . . . . . . . . . . . . . . . 61 5.2 LES at 1% of DNS Resolution . . . . . . . . . . . . . . . . . . . . . 63 5.3 Spanwise Structures . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6 Spectral Results II 83 6.1 LES-V at 3% of DNS Resolution . . . . . . . . . . . . . . . . . . . 83 6.2 LES-V at 1% of DNS Resolution . . . . . . . . . . . . . . . . . . . 85 7 Conclusions 96 7.1 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 96 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Reference List 101 iv List of Tables 3.1 Validation using linear stability theory: vorticity spectral code . . . 40 3.2 Validation using linear stability theory: skew-symmetric spectral code 45 4.1 CTR simulation parameters . . . . . . . . . . . . . . . . . . . . . . 46 4.2 CTR numerical dissipation with filtering . . . . . . . . . . . . . . . 57 4.3 CTR numerical dissipation without filtering . . . . . . . . . . . . . 58 5.1 Spectral simulation parameters . . . . . . . . . . . . . . . . . . . . 60 5.2 Spectral 1% LES performance . . . . . . . . . . . . . . . . . . . . . 69 6.1 Vorticity spectral simulation parameters . . . . . . . . . . . . . . . 84 6.2 Spectral 1% LES-V performance . . . . . . . . . . . . . . . . . . . . 90 v List of Figures 1.1 Laminar separation bubble flow sketch . . . . . . . . . . . . . . . . 3 1.2 Laminar separation bubble on an airfoil . . . . . . . . . . . . . . . . 4 3.1 C p of laminar separation bubble on airfoil and flat plate . . . . . . . 41 3.2 Anatomy of a flat plate laminar separation bubble flow . . . . . . . 42 3.3 Visualization of a flat plate laminar separation bubble flow . . . . . 43 3.4 Flat plate laminar separation bubble flow computational setup . . . 43 3.5 Numerical dissipation example . . . . . . . . . . . . . . . . . . . . . 44 3.6 Spectral LES boundary conditions . . . . . . . . . . . . . . . . . . . 44 4.1 CTR boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Visualization of CTR DNS . . . . . . . . . . . . . . . . . . . . . . . 49 4.3 CTR mean velocity contour plot . . . . . . . . . . . . . . . . . . . . 50 4.4 CTR time-averaged C p and C f . . . . . . . . . . . . . . . . . . . . . 51 4.5 CTR energy rate-of-change in turbulent region . . . . . . . . . . . . 53 5.1 Spectral 3% LES flow visualization . . . . . . . . . . . . . . . . . . 61 5.2 Spectral 3% LES C p and C f . . . . . . . . . . . . . . . . . . . . . . 64 5.3 Spectral 3% LES boundary layer thicknesses . . . . . . . . . . . . . 65 5.4 Spectral 3% LES mean and RMS velocity profiles . . . . . . . . . . 66 vi 5.5 Spectral 1% LES C p and C f . . . . . . . . . . . . . . . . . . . . . . 70 5.6 Spectral 1% LES subgrid-scale dissipation . . . . . . . . . . . . . . 75 5.7 Spectral 1% LES boundary layer thicknesses . . . . . . . . . . . . . 76 5.8 Spectral 1% LES mean velocity profiles in wall units . . . . . . . . 77 5.9 Spectral 1% LES mean and RMS velocity profiles . . . . . . . . . . 78 5.10 Spectral 1% LES maximum RMS velocity . . . . . . . . . . . . . . 79 5.11 Spectral 1% LES velocity auto-correlation functions . . . . . . . . . 80 5.12 Spectral 1% LES velocity auto-correlations for wider domain . . . . 81 5.13 Spectral 3% LES velocity auto-correlation functions . . . . . . . . . 82 6.1 Spectral 3% LES-V C p and C f . . . . . . . . . . . . . . . . . . . . . 86 6.2 Spectral 3% LES-V boundary layer thicknesses . . . . . . . . . . . . 87 6.3 Spectral 3% LES-V mean and RMS velocity profiles . . . . . . . . . 88 6.4 Spectral 1% LES-V C p and C f . . . . . . . . . . . . . . . . . . . . . 91 6.5 Spectral 1% LES-V boundary layer thicknesses . . . . . . . . . . . . 92 6.6 Spectral 1% LES-V mean velocity profiles in wall units . . . . . . . 93 6.7 Spectral 1% LES-V mean and RMS velocity profiles . . . . . . . . . 94 6.8 Spectral 1% LES-V maximum RMS velocity . . . . . . . . . . . . . 95 vii Abstract The flow over blades and airfoils at moderate angles of attack and Reynolds numbers ranging from 10 4 to 10 5 undergoes separation due to the adverse pressure gradient generated by surface curvature. In many cases, the separated shear layer then transitions to turbulence and reattaches, closing off a recirculation region - the laminar separation bubble. To avoid body-fitted mesh generation problems and numerical issues, an equivalent problem for flow over a flat plate is formu- lated by imposing boundary conditions that lead to a pressure distribution and Reynolds number that are similar to those on airfoils. Spalart & Strelets (2000) tested a number of Reynolds-averaged Navier-Stokes (RANS) turbulence models for a laminar separation bubble flow over a flat plate. Although results with the Spalart-Allmaras turbulence model were encouraging, none of the turbulence models tested reliably recovered time-averaged direct numerical simulation (DNS) results. The purpose of this work is to assess whether large eddy simulation (LES) can more accurately and reliably recover DNS results using drastically reduced resolution – on the order of 1% of DNS resolution which is commonly achiev- able for LES of turbulent channel flows. LES of a laminar separation bubble flow over a flat plate are performed using a compressible sixth-order finite-difference code and two incompressible pseudo-spectral Navier-Stokes solvers at resolutions corresponding to approximately 3% and 1% of the chosen DNS benchmark by viii Spalart & Strelets (2000). The finite-difference solver is found to be dissipative due to the use of a stability-enhancing filter. Its numerical dissipation is quan- tified and found to be comparable to the average eddy viscosity of the dynamic Smagorinsky model, making it difficult to separate the effects of filtering versus thoseofexplicitsubgrid-scalemodeling. Thenegligiblenumericaldissipationofthe pseudo-spectral solvers allows an unambiguous assessment of the performance of subgrid-scale models. Three explicit subgrid-scale models – dynamic Smagorinsky, σ, and truncated Navier-Stokes (TNS) – are compared to a no-model simulation (under-resolved DNS) andevaluated against the benchmark DNSdata focusing on two quantities of critical importance to airfoil and blade designers: time-averaged pressure (C p ) and skin friction (C f ) predictions used in lift and drag calculations. Results obtained with these explicit subgrid-scale models confirm that accurate LES of laminar separation bubble flows are attainable with as low as 1% of DNS resolution, and the poor performance of the no-model simulation underscores the necessity of subgrid-scale modeling in coarse LES with low numerical dissipation. ix Chapter 1 Introduction Laminar flow separation, transition to turbulence, and reattachment – coined laminarseparationbubblebecauseoftherecirculationregionthephenomenacreate together – directly affect the performance of increasingly important applications ranging from small unmanned aerial vehicles to low pressure turbines found in jet engines and gas generators. The presence of a laminar separation bubble on a small unmanned aerial vehicle’s wing can increase the lift over drag ratio L/D and thus the efficiency of the craft. Airfoil shapes that maximize this effect are sought at the design stage. On the other hand, laminar separation bubbles cause detrimental unsteadiness and transient structural loads on low pressure turbine blades, so blade shapes and flow control schemes to mitigate this phenomenon are sought at the design stage. Faster, more accurate computational tools are necessary to enable the design optimization of technologies to better control the onset of laminar separation bubbles on blades and airfoils. 1.1 Laminar Separation Bubble Flows The Reynolds number (Re) is a dimensionless ratio of inertial to viscous forces characterizing the flow regime. Reynolds numbers for small unmanned aerial vehi- cles and micro-air vehicles are low to moderate. Based on wing chord length, they aretypically less than2×10 6 andareinsomecases onlyontheorder of 10 4 to 10 5 . By comparison, civilian airplanes are characterized by Reynolds numbers ranging 1 from a few million to 80× 10 6 for the Boeing 747 at cruising velocity. Recent experimental investigations of low Reynolds number aerodynamics by Hu et al. (2007);Hain et al.(2009);Spedding & McArthur (2010) reveal that low tomoder- ate Reynolds number flows over airfoils are often dominated by the effects of flow separation and reattachment. The appearance of unsteady recirculation regions due to separation and reattachment greatly influence the aerodynamic forces the wing is subjected to. They change the lift and drag characteristics and thus the flight stability of small unmanned aircraft and micro aerial vehicles. Separation- induced transition is also seen in higher Reynolds number flows. For example, flow separation and turbulent reattachment is sometimes observed on the blades of wind turbines as well as on low pressure gas turbine blades due to the influence of the interior wall. This phenomenon causes unsteadiness in the flow, which is a determining factor in high cycle fatigue of turbomachinery components. Although most wind and gas turbines’ operating Re are higher than 10 5 , the physics gov- erning separation-induced transition and the stability of the recirculation region it creates are akin to those seen in laminar separation bubbles that occur on lower Re applications like micro aerial vehicles. A number of experiments have been performed to elucidate laminar sep- aration bubble flows. The seminal works by Gaster (1963) and Horton (1968) established the foundation for the physical understanding of the lam- inar separation bubble and advanced some semi-empirical predictions for the location of separation and reattachment. Since then, a number of experi- ments have been conducted to study in more detail the structure, stability, and dynamics of laminar separation bubble flows Alving & Fernholz (1996); Häggmark(2000);Marxen et al.(2003);Burgmann et al.(2006);Yarusevych et al. (2006); Burgmann et al. (2007); Burgmann & Schröder (2008); Hu et al. (2007); 2 Figure 1.1: Features of the flow field in the vicinity of a transitional separation bubble Horton (1968); Lin & L.Pauley (1996); Castiglioni et al. (2014). Hain et al. (2009); Spedding & McArthur (2010). As a result, the physical origin of laminar and transitional flow separation is now qualitatively well understood. Asfigure1.1 illustrates, theattached laminar boundarylayer developing onawing or blade is subjected to an adverse pressure gradient due to the airfoil’s curvature, which causes it to separate. Immediately behind the separation point there is an effectively stagnant flow region, the so-called dead air region. A reverse flow vor- tex develops downstream of the dead air region, leading to an inflectional mean velocity profile in the boundary layer. This triggers the growth of convective and secondary instabilities which quickly break down to turbulence. Kelvin-Helmholtz rollsareindeedvisibleinasnapshotoftheflowfieldofalaminarseparationbubble flow over a NACA0012 airfoil shown in figure 1.2. As the separated shear layer transitions to turbulence, its interaction with the reverse flow vortex causes it to 3 Figure 1.2: Instantaneous iso-surfaces of spanwise vorticity on a NACA0012 airfoil at 5 degrees angle of attack and Re c = 5×10 4 taken from a STAR-CCM+ large eddy simulation Castiglioni (2015). reattach, thereby closing off the recirculation region. Aft of the laminar separa- tion bubble, unlike in the two dimensional case, the flow is fully turbulent and three-dimensional. Clear spanwise vortices are not shed, but the size and shape of the bubble changes in time due to the complex balance between the effects of the pressure gradient, convective instabilities, and viscous dissipation. This pic- ture emerges from these experimental investigations as well as from direct numer- ical simulations (DNS) results Lin & L.Pauley (1996); Spalart & Strelets (2000); Alam & Sandham (2000); Marxen & Rist (2010); Jones et al. (2008, 2010). 4 Many DNS were carried out to match and augment experimental data. These DNSshedlightonthemechanisms ofenergytransferatworkinlaminarseparation bubble flows Skote et al. (1998); Spalart & Strelets (2000); Skote & Henningson (2002), the process of transition Na & Moin (1998); Wu & Moin (2010); Alam & Sandham (2000); Marxen & Rist (2004), the stability characteristics of such flows Marxen et al. (2003); Marxen & Rist (2010); Jones et al. (2010), and theeffectofdisturbancesandforcingonitsdynamicsHerbst & Henningson(2006); Marxen & Henningson (2011); Jones et al. (2008). Agreement with experiments was found to be very favorable in most cases. Key to getting agreement with experimental results was resolving the reverse flow region near the wall, and the shear layer which transitions to turbulence above the separation bubble. Although three-dimensional laminar separation bubble flows have been shown to differ substantially from their two-dimensional analog, there is some disagree- ment on the degree of importance of spanwise structures, their dominant wave- length, and the role they play in laminar separation bubble flows. Recent simula- tionresults show that airfoil sections may require one chord length in the spanwise directiontoensurethefullthree-dimensionalnatureofthelaminarseparationbub- ble can be captured Eisenbach & Friedrich (2008). Experimental studies confirm that the dynamics and turbulent statistics of laminar separation bubble flows are sensitive to incoming levels of turbulence, boundary conditions, and even acoustic vibrations. This impedes the direct comparison of different experiments, and the search for universal empirical parameters and predictors for its impact on lift and dragaswellasthelevelofturbulenceitgeneratesintheflow. Researchers devising technologies to control or mitigate the onset of separation bubbles face the same issues. A predictive tool fast enough to reliably explore the sensitivity to these parameters and flow control schemes is needed. 5 1.2 Motivation In order to produce more efficient airfoil or blade designs, to create control schemes to reduce separation effects, and to better predict high cycle fatigue, numerical prediction tools for laminar separation bubble flows are needed. Airfoil and blade designers are primarily interested in obtaining accurate time-averaged quantities relating to the aerodynamic forces on the airfoil or blade, namely the coefficient of pressure to measure the lift and drag, and the skin friction, another important component of drag. To enable optimization in a realistic industrial setting, thesequantities must becalculatedinamatterofafewhoursorless. Such computationallyaffordable,accurateandreliablenumericalpredictionsforlaminar separation bubble flows had not been obtained until this study. As a proof of conceptandtoavoidnumericalissuesassociatedwithbody-fittedmeshgeneration, anequivalent problemforflowoveraflatplateisformulatedbyimposingboundary conditions that lead to a pressure distribution and Reynolds number based on bubble length similar to those observed on the airfoils of small unmanned aircraft and low-pressure turbine blades. The objective of this research is thus to test the accuracy of large eddy simulation at resolutions drastically reduced compared to a benchmark DNS for a laminar separation bubble flow over a flat plate. Corollary objectives are to identify key factors in obtaining accurate predictions, and to investigate the performance of different subgrid-scale models in this highly under- resolved environment. 6 Chapter 2 Background 2.1 Governing Equations The relevant governing equations are the incompressible Navier-Stokes equa- tions because the fluid density and temperature may be considered constant for laminar separation bubble flows at moderate Reynolds numbers (Re = 10 4 to Re = 2×10 5 ). ∂u i ∂t + ∂ ∂x j (u i u j ) =− ∂p ∂x i +ν ∂ ∂x j ∂u i ∂x j + ∂u j ∂x i ! +F i , i = 1,2,3 (2.1) ∂u i ∂x i = 0, (2.2) The velocity fieldu = (u 1 ,u 2 ,u 3 ) expressed in reference to a Cartesian coordinate system x = (x 1 ,x 2 ,x 3 ) is a solution of the momentum and continuity equations, (2.1) and (2.2). p = P/ρ is the static pressure and ν is the kinematic viscosity which is assumed to be constant and uniform in space. No accepted analytical closed form solution to the Navier-Stokes equations for laminar separation bubble flows is known. A numerical solution to the Navier- Stokes equations is thus sought on a discrete set of grid points x ijk or discrete cell volumes. A number of approaches exist to discretize the equations such as the finite difference, finite volume, and finite element methods to name only a few. Despite fundamental differences in their approach, each discretization method’s conservation properties are directly dependent on the accuracy of the schemes 7 used to approximate each of the terms in the equations solved. The accuracy of the solution depends both on the exact formulation of the numerical scheme, the quality (level ofnon-orthogonality)ofthegrid,andthelevel ofresolutiononwhich the discretized equations are solved. 2.2 Direct Numerical Simulation Once discretized, the Navier-Stokes equations are integrated in time. Doing so without further approximations is called direct numerical simulation (DNS). DNS is the most accurate and reliable fluid dynamics simulation tool available. However, to obtain an accurate solution using DNS requires that the mesh or cell size captures all relevant scales of motion in the problem. If the flow to be simulated is turbulent, then the mesh or cell size must be within one order of magnitude from the Kolmogorov length scale, the smallest length scale at which turbulenceexists. Dominatedbyviscosity,thesesmallscalesarelargelyresponsible for the dissipation of turbulent motion into heat. Capturing Kolmogorov scales imply solving the Navier-Stokes equations on a very large number of cells or a very fine mesh. The computational work required may take a prohibitively large amount of time even on the best available supercomputer. Indeed, DNS requires substantial computational resources, long wall-clock runs, and long analysis times; e.g.a relatively simple 3-D airfoil configuration at a Reynolds number of 5×10 4 required over 170 million grid points Jones et al. (2008). For laminar separation bubble flow over a flat plate at a Reynolds number of 10 5 , DNS results required over 16000 processor-hours Cadieux et al. (2012). A number of 3-D configurations and angles of attack need to be quickly investigated to allow for the optimization 8 of airfoil and turbine blade designs. For this case, a DNS approach is generally impractical and other simulation approaches must be considered. 2.3 Reynolds-averaged Navier-Stokes A widely-used computationalfluid dynamics simulation approach for moderate to high Reynolds number flows is to solve for the mean flow quantities directly instead of calculating the primary quantities at each instant in time as is done in DNS. The primary quantities are split into fluctuating and mean flow components u i =U i +u ′ i (2.3) where U i =hu i i is the ensemble average of u i satisfying Reynolds’ averaging con- ditions. Substituting eq.(2.3) into the Navier-Stokes equations and noting that hu ′ i i = 0 gives the Reynolds-averaged Navier-Stokes (RANS) equations for incom- pressible flow: ∂U i ∂t + ∂ ∂x j (U i U j ) =− 1 ρ ∂P ∂x i + ∂ ∂x j ν ∂U j ∂x i + ∂U i ∂x j ! −τ rans ij ! , (2.4) ∂U i ∂x i = 0 (2.5) where τ rans ij =hu ′ i u ′ j i. (2.6) All information about the instantaneous velocity fields is lost through this averag- ing in favor of obtaining an estimate of the mean flow. Unfortunately, averaging procedures used in practice, namely time and spatial averaging, do not strictly satisfy Reynolds last condition that hhf(x)ig(x)i=hf(x)ihg(x)i for any flow that is not fully developed with a large separation of scales, limiting its applicability to 9 such flows Wilcox (2006). In other words, if the mean and fluctuating components are correlated, then the time or spatial average of their product does not vanish and the RANS equations are no longer valid Wilcox (2006). As seen in eq. (2.4), the RANS equations have an additional unknown term, the Reynolds stress tensor τ rans ij , which must be modeled using prior knowledge about the flow being simu- lated and/or additional equations. The majority of RANS models forτ rans ij rely on the Boussinesq eddy viscosity approximation. The Boussinesq hypothesis states that by analogy to momentum transfer in the molecular motion of a gas which can be described by a molecular viscosity ν, the Reynolds stress tensor should be proportional to the mean strain rate tensor using a turbulent eddy viscosity ν rans : τ rans ij ∝ν rans ∂U j ∂x i + ∂U i ∂x j ! (2.7) This reduces the number of extra unknowns from six to one, the eddy viscosity. Most RANS models use dimensional arguments and analogy to other physical processes to set the value for ν rans . For example, Prandtl’s mixing length model uses the Boussinesq analogy to molecular momentum transport and assumes that there exists a turbulent mixing length analogous to the mean free path Wilcox (2006). Using dimensional analysis, the mixing length model for a boundary layer simplifies to τ rans =ν rans ∂U ∂y , (2.8) ν rans =ℓ 2 mix | ∂U ∂y |, (2.9) ℓ mix =δ ∗ = Z ∞ 0 1− U(y) U ∞ ! dy ≈ 1.72x/ q Re x , (2.10) 10 where the mixing length is assumed to be proportional to the displacement thick- ness of the boundary layer. Such a simple turbulence model provides results in agreement with experiments for boundary layers Wilcox (2006). However, it is incomplete because the appropriate mixing length for the flow being simulated must be known a priori. Moreover, the mixing length model is only valid for sim- ple flows with slow-varying properties (so-called equilibrium turbulent flows) due to the assumptions used in deriving it. Modern models are far more complex and perform better in a wider range of flows, but still struggle with accurately pre- dicting transition from laminar to turbulent flow and the lengths of recirculation zones. 2.3.1 RANS for Laminar Separation Bubble Flows Laminar separation bubble flows present a challenge to RANS models because boundary layer separation and reattachment involves subtle interactions between viscous, advective, andpressureeffectsandisinherentlyanon-equilibriumprocess, especially when driven by an adverse pressure gradient instead of geometry (e.g. backward facing step). Although the rate of strain changes rapidly as the flow separates, the turbulence adjusts to changes in the flow on an unrelated, longer time scale. A perturbed turbulent boundary layer was experimentally shown not to return to equilibrium for at least 10 boundary-layer thicknesses downstream of the perturbation Wilcox (2006). Since the Reynolds stresses modeled in typical two-equation RANS models are adjusted based solely on the mean rate of strain, they preclude any such transient effects from flow history. Attempts were made to adjust forthiseffect byrelaxing theeddyviscosity intheregionbehindseparation. However, this method requires prior knowledge of the separation point - and thus 11 precludes strong predictive capabilities even for unsteady RANS Howard et al. (2000). Spalart & Strelets (2000) tested a number of typical RANS turbulence model, namely Spalart-Allmaras, Menter’s shear stress transport, modified shear stress transport, and Secundov’s ν 92 t , for a simple laminar separation bubble flow on a flat plate driven by suction from the top. All RANS turbulence models except Spalart-Allmaraspredictearliertransitionandreattachmentthanobservedintheir spectral DNS results. Obvious disagreement in the location and magnitude of the peak negative skin friction is observed between the results of different models. All RANS results were found to under-predict the level of skin friction down- stream of reattachment. Without modification, Menter’s shear stress transport model transitions even before the expected separation point, predicting attached flowthroughout thedomain. AlthoughSpalart-Allmarasresults wereencouraging, none of the RANS turbulence models Spalart & Strelets (2000) tested recovered all the important features of the time-averaged skin friction DNS curve: location of the separation point, location and magnitude of the peak negative skin friction, reattachment point, and turbulent skin friction levels immediately downstream of reattachment. The large variations observed in skin friction predictions depending on the RANS turbulence model chosen may be symptomatic of a lack of robust- nessofthetypicalRANSturbulencemodelingapproachwhenappliedtoinherently unsteady phenomena. To better predict transition, Howard et al. (2000) sensitized two-equation RANS model coefficients of the Launder andSharma andk−g models to the local turbulent Reynolds number as proposed by Wilcox (2006), removing the need for a priori knowledge of the separation point. With the local turbulent Reynolds 12 number sensitivity modification, two-equation models were shown to predict tran- sition more reliably in an unsteady RANS solver and showed improved agreement for separation bubble length. However, Launder and Sharma and k − g model results were still inconsistent amongst each other and amongst different sensitiza- tion approaches. None matched the DNS results for peak negative skin friction and its level immediately downstream of reattachment Howard et al. (2000). RANS turbulence models using second moment closures show improvement over single and two-equation models without the need for artificially triggering transition, but require tuning of their closure coefficients for lower Reynolds num- ber flows with separation and transition Hadžić & Hanjalić (2000). While results show better agreement for the reattachment point and turbulent region at coarser resolution than the Spalart-Allmaras model, second-moment closure RANS results still fail to recover the peak negative skin friction. Since the authors mention but do not present results for a finer grid, it can be surmised that results did not improve significantly Hadžić & Hanjalić (2000). That the method does not con- verge to DNS results with increased resolution is further evidence that it is not well-suited for the problem. Similarly, results for a laminar separation bubble on a flat plate with a semi- circular leadingedgeusingatwo-layer modelmatched experimental resultsreason- ablywell, butrequired anumber ofempiricalcorrelationalongwithspecialmodifi- cations to the model to capture transition Papanicolaou & Rodi (1999). Although RANS can predict separation and transition reliably with models optimized for such flows (second-moment closure, two-layer model), it still struggles to recover wall skin friction accurately for this proof-of-concept laminar separation bubble flow over a flat plate. 13 2.4 Large Eddy Simulation Anotheroptionistoemploylargeeddysimulation(LES)techniques. LEStech- niques were developed based on the observations made by Kolmogorov that the smallest scales of turbulence are dominated by viscosity, behave mostly isotropi- cally and account for the majority of turbulent dissipation. The understanding of the energy cascade from large to small scales led to the idea of modeling the more universal small scales of turbulent motion while resolving the energy containing eddies directly affected by the flow boundary conditions. This greatly relaxes the DNSrequirementthatthemeshorcellsizebeonthesameorderastheKolmogorov length scale. LES hinges on the use of subgrid-scale models to predict the correct small scale dissipation rate based solely on information from the larger resolved scales. 2.4.1 Filtered Navier-Stokes Equations A low-pass filter operation is used to separate small scales from larger ones and derive the filtered Navier-Stokes equations. This filtering operation can be described by a convolution integral with the filter function or kernel G, u i (x,t) =G∗u i = Z +∞ −∞ G(x−x ′ ;∆) u i (x ′ ,t) d 3 x ′ (2.11) where the resolved scale or filtered velocity is ¯ u i , and the subgrid-scale velocity is defined as u ′ i = u i − ¯ u i . ∆ is the filter width which is generally taken to 14 be proportional to the grid or cell size ∆ = (∆ x∆ y∆ z) 1 3 Sagaut (2006). The incompressible filtered Navier-Stokes equations can then be written as follows: ∂¯ u i ∂t + ∂ ∂x j (¯ u i ¯ u j +τ sgs ij ) =− ∂¯ p ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j + ∂¯ u j ∂x i ! +F i , (2.12) ∂¯ u i ∂x i = 0, (2.13) where τ sgs ij =u i u j − ¯ u i ¯ u j (2.14) isthesubgrid-scalestresstensor. Itcontainsthetermu i u j whichisanewunknown. Just as turbulence models are required to close the RANS equations, the fil- tered Navier-Stokes equations require a subgrid-scale model for τ sgs ij . However, the derivation of subgrid-scale models generally rely only on the assumption that the effect of the scales of motion smaller than the filter width on the larger scales are small and mostly dissipative. This assumption is more robust than those made in deriving the RANS equations and most RANS turbulence models because it is satisfied for a much wider array of flow conditions and relative filter widths. 2.4.2 Subgrid-scale Models There exists avarietyof different subgrid-scale models toclose equations (2.12) and (2.13) that generally belong to one of two distinct categories: structural and functional modeling. Structural modeling directly approximates the subgrid-scale stress tensor or subgrid-scale velocities based on the resolved velocities or a formal series expansion. This approach assumes that the structure of the small scales is universal and the energy contained in the subgrid-scales are a function of the 15 resolved scales Sagaut (2006). Examples include models based on the deconvolu- tion procedure, stress transport models, and subgrid-scale velocity reconstruction models. Arecentexampleisthevelocity estimationmodelbyDubois et al.(2002). Instead of directly approximating τ sgs ij , functional modeling seeks to approxi- mate the effects of inter-scale energy transfer on the resolved scales. Rather than assuming that the structure of the small scales is universal, this approach posits thattheeffectsofthesmallscalesonthelargerresolved scalesareuniversalSagaut (2006). Knowledge of the turbulent energy cascade and the concepts of forward and backscatter are used to justify the assumption that the effect of the small scales on the large is universal and depends only on the energy of the large scales driving the flow. Recent examples of such subgrid-scale models are the σ and the interscale energy transfer models Nicoud et al. (2011); Anderson & Domaradzki (2012). The increasingly popularimplicit LES(ILES) approachalso generally falls into this category, where the numerical scheme is adjusted such that its truncation errors and associated dissipative and dispersive effects have the desired impact on the resolved scales. The most famous functional modeling example remains the Smagorinsky model. The Smagorinsky Model The Boussinesq approximation used in many RANS models is invoked again, but this time to calculate a turbulent eddy viscosity ν sgs describing the subgrid- scale dissipation based solely on the resolved scales of motion. τ sgs ij − 1 3 τ sgs kk δ ij = 2ν sgs ( ¯ S ij − 1 3 ¯ S kk δ ij ), (2.15) 16 where ν sgs = (C S ∆) 2 | ¯ S| ¯ S ij , (2.16) | ¯ S| = (2 ¯ S ij ¯ S ij ) 1 2 (2.17) ¯ S ij = 1 2 ∂¯ u i ∂x j + ∂¯ u j ∂x i ! , (2.18) ∆ = (∆ x∆ y∆ z) 1 3 . (2.19) C S is a closure coefficients determined ahead of simulation by matching avail- able experimental or DNS data. This is a limitation to the model’s applicability because foreknowledge of the flow being simulated is required to set the closure coefficient. It also cannot account for the local reduction of eddy viscosity near physical boundaries without the use of explicit Van Driest damping functions. A modification of the Smagorinsky model which eliminates these issues and makes it universal consists of letting the closure coefficients be functions of time and space (e.g. C 2 S = C d (x,t)). In this dynamic version of the Smagorinsky model, these coefficients are computed dynamically using the Germano identity. The Dynamic Procedure The dynamic procedure used to compute the local instantaneous closure coeffi- cient C d (x,t) for the Smagorinsky model is used here as an example, but it can be applied to other models formulated in a similar fashion. Its purpose is to provide this coefficient based on local resolved strain rate without prior knowledge of the flow. It effectively reduces the model contribution in laminar shear flows where 17 the stress tensor is not zero, but no turbulence exists, removing the issues that the static Smagorinsky model faces with physical boundaries. ν sgs =C d ∆ 2 | ¯ S| ¯ S ij (2.20) C d = min 0.2 2 , max " hL ij M ij i hM ij M ij i , 0 #! , (2.21) L ij = d ¯ u i ¯ u j − b ¯ u i b ¯ u j , (2.22) M ij = 2 \ ∆ 2 | ¯ S| ¯ S ij −2 b ∆ 2 | b ¯ S| b ¯ S ij , (2.23) b ∆ = 2∆ = 2(∆ x∆ y∆ z) 1 3 . (2.24) The overbar is used to represent grid-filtered terms (filter width ∆) and the hat is used to indicate test-filtered quantities using Simpson’s rule (filter width 2∆): b f(x)≈ 1 6 f(x−∆ x)+ 2 3 f(x)+ 1 6 f(x+∆ x). (2.25) The weights of the filter are adapted to the non-uniform vertical grid using quadratic interpolation. The h i symbol used in (2.21) denotes averaging in any uniformdirection(ifoneexists) orlocalaveraginginalldirections. Althoughsmall negative values of C d may be physically justified, mimicking backscatter phenom- ena, such negative values can quickly lead to numerical instability. So in practice, averaging is used in tandem with clipping to avoid negative as well as rapidly oscillating values of C d . The constant computed through this dynamic procedure is local in space and time. It effectively reduces the model contribution in lami- nar shear flows where the stress tensor is not zero, but no turbulence exists. The dynamic procedure is particularly computationally expensive due to the applica- tionofatest-filter inthree-dimensions onat least two tensorialquantities d ¯ u i ¯ u j and \ ∆ 2 | ¯ S| ¯ S ij as well as local spatial averaging, each filtering and averaging operation 18 often requiring communication among different processes or blocks in parallelized implementations. Despite this shortcoming, the dynamic Smagorinsky model has become the benchmark against which other subgrid-scale model are tested due to its success in academia and its universality. The σ Model The σ model follows definitions set forth in (2.15) and (2.18) but computes ν sgs using the singular values σ i of the velocity derivative tensor g ij . This choice is motivated by the desire to improve on the dynamic Smagorinsky model by provid- ing more appropriate near-wall behavior, as well as providing zero contributions in pure two-dimensional shear or pure rotation cases Nicoud et al. (2011). ν sgs = (C σ ∆) 2 σ 3 (σ 1 −σ 2 )(σ 2 −σ 3 ) σ 2 1 is the subgrid-scale eddy viscosity, (2.26) where σ 1 ≥σ 2 ≥σ 3 ≥ 0, are the singular values of g ij = ∂¯ u i ∂x j , and (2.27) C σ = 1.35 is the closure coefficient. (2.28) C σ is determined from homogeneous turbulence and validated using channel flow simulations Nicoud et al. (2011). The singular values σ i are obtained using the invariants of G ij = g ki g kj and their angles to avoid the overhead of linear algebra library calls to an eigenvalue solver for each cell in the domain at each time step. 19 The truncated Navier-Stokes Approach The truncated Navier-Stokes approach follows the method developed by Domaradzki et al. (2002) in which periodic filtering is used as a substitute for a subgrid-scale model. Periodic filtering is used to remove energy from the small- est resolved scales by the use of a low-pass approximate deconvolution method (ADM) filter Stolz et al. (2001). The filtering operation is implemented using the product of an approximate deconvolution filter Q N ≈ G −1 described in Tantikul & Domaradzki (2010) with filter G: Q N G =I−(I−G) N+1 . (2.29) The orderN = 5 is chosen such that the filter only affects scales smaller thanfilter width ∆ = ∆ x when using a simple three point filter in physical space G(∆ x): G(∆ x)∗f(x)≈ 1 8 f(x−∆ x)+ 3 4 f(x)+ 1 8 f(x+∆ x). (2.30) The filter weights are adjusted for the non-uniform vertical grid using quadratic interpolation. Sincetheenergycascadesfromlargetosmallscalesandaccumulates there slowly in a high order under-resolved simulation, it is only necessary to filter afterafractionofapercent ofthelargeeddyturnover time. Hence, filteringisonly applied when the kinetic energy at high wave numbers reaches unphysical levels in the truncated Navier-Stokes approach. This is fundamentally different from implicit LES, where the inherent approximation errors can be said to act as a low- pass filter at each time step. Excessive energy accumulating in the small scales is detected using a criterion based on the ratio of energy removed I(∆ x)/I(2∆ x) by 20 two ADMtest-filters with filter widths ∆ = ∆ x denoted bythe tildeand∆ = 2∆ x denoted by the hat as follows I(∆ x) I(2∆ x) = Z V E− ˜ E E− b E dV (2.31) ≈ Z Y 0 * P 3 i=1 1 2 (u i − ˜ u i )(u i − ˜ u i ) P 3 i=1 1 2 (u i −b u i )(u i −b u i ) + (y)dy (2.32) ˜ u i = (Q 5 G(∆ x))∗u i (2.33) b u i = (Q 5 G(2∆ x))∗u i (2.34) G(2∆ x)∗f ≈ 1 4 f(x−∆ x)+ 1 2 f(x)+ 1 4 f(x+∆ x). (2.35) The ratio I(∆ x)/I(2∆ x) represents the energy contained in the small scales com- pared to the larger scales. When it reaches values in excess of those obtained for a typical dissipation, inertial, or Batchelor energy spectrum – determined to be 0.007 to 0.009 from theory by Tantikul & Domaradzki (2010, 2011) – primary variables are filtered in physical space with filter Q 5 G(∆ x). Using this criterion, thefilter isappliedatvarying intervalscentered around200hundred timestepsfor the coarsest resolution simulation presented here, corresponding to approximately 0.5% of one non-dimensional time unit t =t ∗U 0 Lx . 2.4.3 LES Results for Laminar Separation Bubble Flows A number of LES of laminar separation bubble flows over flat plates and airfoils have been completed recently by Wilson & Pauley (1998); Yang & Voke (2001); Roberts & Yaras (2005); Eisenbach & Friedrich (2008); Xu et al. (2010); Kojima et al.(2013). Forinstance, LESresultsYang & Voke(2001) obtainedwith thedynamicSmagorinskymodelwerereportedtobeingoodagreementwithexper- iments for boundary-layer separation and transition caused by surface curvature 21 at Re = 3,450. Yet even for this relatively low Reynolds number, the two critical issues ingettingagreement were anumericalresolution (472×72×64 meshpoints) comparable to DNS of the same flow, and a high order numerical method. Such strict requirements are difficult to satisfy in simulations of practical flows often performed with low order finite difference or finite volume methods (e.g. commer- cial codes). Similarly, LES of flow separation on anairfoilat a high angleof attack was performed at Re = 10 5 using Cartesian grids Eisenbach & Friedrich (2008). This case also required very high resolutions between 50 and 100 million mesh points to obtain good agreement. Using LES with such high resolution and higher order methods implies a time-to-solution on the same order as DNS. Therefore, the question remains: can LES produce sufficiently accurate results for laminar separation bubble flows with drastically reduced resolution, around 1% of DNS resolution, commonly achievable for fully turbulent flows? While a handful of other investigators also performed low resolution LES, the effect of different subgrid-scale models on the quality of important results such as time-averaged skin friction and pressure coefficient remains largely unknown for laminar separation bubble flows. Only the constant Smagorinsky model has been compared to its dynamic counterpart and a no-model case in Wilson & Pauley (1998). Other investigators like Eisenbach & Friedrich (2008); Yang & Voke (2001); Xu et al. (2010) relied entirely on the dynamic Smagorinsky model, or on the numerical dissipation of their chosen scheme as in implicit LES or ILES in Kojima et al. (2013) and even without any prior knowledge of the dissipative scheme’seffectsontheresolvedscalesinthecaseofRoberts & Yaras(2005). With- out benchmark DNS data or a baseline case with no subgrid-scale model active to compare directly to, the performance of their models or ILES results could not be evaluated quantitatively. 22 2.5 Wall-modeled LES and Hybrid RANS-LES For LES to be accurate, it requires a mesh nearly as fine as a DNS near physi- cal boundary where boundary layers develop. To mitigate this stringent resolution requirement, wall models are developed to give approximate boundary conditions to the LES solver away from the surface. A wide variety of models with differ- ent assumptions have been proposed and reviewed by Piomelli & Balaras (2002); Piomelli (2008) and Sagaut & Deck (2009). The use of wall models generally mitigates noise generated by poor approximation of curved surfaces or highly non- orthogonalbody-fittedmeshes inlow resolution settings andpermit simulations to reach much higher Reynolds numbers that are closer to operating conditions for mostturbomachinerybladesMcMullan & Page(2012). However, wallmodelshave historicallyhadapoortrackrecordinpredictingseparationandreattachment, and introduce another source of error into LES due to further approximations made at the wall Bose & Moin (2014). But more importantly, when wall models are used it becomes difficult to distinguish between the performance of the wall model and thatofthesubgrid-scalemodelbecauseobtainingcorrectamountofturbulent con- tent near the wall is key to the overall accuracy of LES. Detached eddy simulation (DES) solves the RANS equations near the wall and smoothly transition to LES on a single grid using a single hybrid RANS-LES turbulence model developed by Spalart(2006). Thisremovestheneedforwall-layermodelingwhilestilldrastically reducing near-surface resolution requirements, allowing the simulation of higher Reynols number flow. Despite many successes, difficult issues such as modeled stress depletion in the log law region, and non-monotonic grid convergence arise in DES as well as its derivatives delayed DES and zonal DES as acknowledged by Spalart (2006, 2009); Deck et al. (2011); Deck (2012). Partially averaged Navier- Stokes(PANS)developedbyGirimaji & Abdol-Hamid(2005);Basara et al.(2011) 23 and other variable resolution approaches such as the scale-adaptive simulation (SAS) proposed by Menter & Egorov (2010); Egorov et al. (2010) and turbulence- resolving RANS (TRANS) put forth by Shur et al. (2008) all avoid these issues by tyingtheeddyviscositytophysicalquantitieslikeenergygriddensityoranintegral length scale. However, their turbulence modeling approaches and by consequence their results are often closer to unsteady RANS than LES and are thus limited in their ability to capture transient and unsteady effects accurately as pointed out in Menter & Egorov (2010). In fact, neither hybrid approaches nor wall models have been validated to the same extent as LES, even when in pure RANS or pure LES mode(whenthemethodallowsit)Sagaut & Deck(2009). Toinvestigatetheeffects of different subgrid-scale models without the unknown influence of wall models or hybrid RANS-LES approaches, wall-resolved LES is chosen for this work. 24 Chapter 3 Methodology 3.1 Approach Laminar separation bubble flows occur on blades and airfoils at low to mod- erate Reynolds numbers ranging from 10 4 to 10 5 due the curvature of the airfoils and blades. Simulating flow over blades and airfoils requires the creation of non- orthogonal body fitted meshes, unstructured grids, or the use of immersed bound- ary methods to properly represent the airfoil or blade’s surface. Grid creation not only presents its own challenges, but also often limits numerical solution methods tosecondorderaccuracyinspaceandtime,withtheexceptionoffiniteelementand discontinuous Galerkin methods. The choice of meshing technique may also have inherent approximations that in turn affect stability and accuracy of the numeri- cal methods at low resolution. For example, approximating a curved surface by a series ofconnected straightlines asopposedtobezier curves orsplines cangiverise to inaccurate results in second order codes, and catastrophic numerical instability when using higher order methods. In order to investigate the capacity of LES to reduce the resolution require- ment for laminar separation bubble simulations and the performance of different subgrid-scale models free from the numerical issues associated with geometry, flow over a flat plate with an adverse pressure gradient strong enough to cause separa- tionas described insection 3.2 is studied. This approach has been used sucessfully tostudylaminarseparationbubbleflowsbothinexperimentsbyHäggmark(2000); 25 Marxen et al. (2003); Sohn et al. (1998) and in simulations by Spalart & Strelets (2000); Alam & Sandham (2000); Wilson & Pauley (1998); Herbst & Henningson (2006); Na & Moin (1998); Skote et al. (1998); Wu & Moin (2010). The resulting pressure distribution is qualitatively comparable to what is seen on blades and airfoils as is shown in figure 3.1: a smooth increase in pressure is followed by a plateau over the separation bubble. The plateau ends with a sharp rise in pressure indicating the flow’s transition to turbulence and reattachment. Downstream of the sharp rise, the pressure plateaus again over the developing attached turbulent boundary layer. The only difference to note between the airfoil and the flat plate laminar separation bubble pressure distribution is that the sharp peak in pressure at the stagnation point of the airfoil is not observed on the flat plate. Despite this difference, the Reynolds number based on bubble length (Re ≈ 67,000) is similar to those found on airfoils and blades indicating a degree of physical equiva- lence. In fact, flat plate laminar separation bubble flows display the same physical features as those seen on airfoils as evidenced by the similarity of the mean veloc- ity contours and mean velocity profiles in figure 3.2 to those in figure 1.1. The mechanisms for separation-induced transition are also the same as in the airfoil case. Notice that the Kelvin-Helmholtz rolls visible on a flat plate with suction from the top shown in figure 3.3 near x = 3.5 closely resemble those seen on an airfoil shown in figure 1.2, each figure displaying iso-surfaces of spanwise vorticity. As such, general conclusions reached from investigating flat plate separation and reattachment should also be applicable to blades and airfoils. 26 3.2 Flow Specification ThecomputationalsetupusedbySpalart & Strelets(2000)tostudyseparation- induced transition flow over a flat plate is followed. The physical domain is a rectangular box with height Y, length 7.5Y, and width 0.6Y (see figure 3.4). At the inflow a laminar Blasius boundary layer velocity profile is imposed with the freestreamvelocityU 0 . Atthetopboundary, averticalsuctionvelocity isimposed in a narrow slot oriented perpendicular to the mean flow direction. The suction produces an adverse pressure gradient that causes flow separation. The flow then transitionstoturbulence andreattachesdownstream. The verticalsuctionvelocity is specified as V(x) =aexp(−[(x−x s )/(0.24Y)] 2 ), (3.1) where a is the peak velocity and x s is its streamwise location Spalart & Strelets (2000). The resulting separation bubble is sensitive only to the upper-wall bound- ary conditions through the nominal flow deceleration parameter S, S = 1 YU 0 Z V(x)dx. (3.2) Using the height Y to non-dimensionalize all relevant lengths the parameters in the equations above are set such that x s = 3, S = 0.3 and the Reynolds number at x s is Re xs = 10 5 , giving a ≈ 0.7U 0 and Re Y = Re xs /3, matching those in Spalart & Strelets (2000). These choices are driven by the requirement that the flow separates naturally, without additional forcing mechanisms like those used in Alam & Sandham (2000). 27 3.3 Numerical Methods The numerical schemes used to approximate the equations in space and inte- grate them in time have an impact on three important quantities: the rate of con- vergence, numerical dissipation and dispersion. The choice of numerical scheme determines the rate of convergence to the true solution. For example, a second order scheme in space implies that doubling the number of mesh points over a given solution domain should decrease the error of the simulation by a factor of 4 over a constant time of integration. The concepts of numerical dissipation and dispersion are linked to the exact form of the truncation error terms of any approximation made by the scheme. For example, a spatialderivative approximated by acentraldifference hasa truncation error E calculated from its Taylor series expansion as follows. ∂ x f(x) = f(x+h)−f(x−h) 2h +h 2 ∂ 3 x f(x)+O(h 3 ) (3.3) ≈ f(x+h)−f(x−h) 2h (3.4) E =h 2 ∂ 3 x f(x)+O(h 3 ) (3.5) The omission of these higher order terms in the simulation have effects that are unknown apriorianddependonthegoverning equations, theparticularflowsimu- lated and the degree of under-resolution. Ascheme is described as dissipative if its totalenergykineticdecreasesfasterthantheexactsolution. Itsvisibleeffectcanbe likened to artificially increasing viscosity. The classic example is the reduction in amplitude of a half sine wave, and its increasing wave length during linear convec- tion in space using a simple finite-difference upwind scheme as seen in figure 3.5a reproduced from Fletcher (1991). Numerical dispersion is linked to the scheme 28 amplifying and attenuating different Fourier modes of a derivative approximation causing oscillations that travel with different wave speeds. For example, the same half-sine wave linear convection problem solved using a Crank-Nicolson scheme results in spurious oscillations as seen in figure 3.5b reproduced from Fletcher (1991). Such spurious oscillations are characteristic of a scheme with a dominat- ing dispersive term. Dissipative schemes have been preferred historically for their robustness: their ability to artificially smooth out any sharp changes or disconti- nuities insolutionswherelessdissipative ormoredispersive schemes might become numerically unstable and not provide any results. This is problematic for three reasons. Solutions for flows with shocks where discontinuities are physical will be increasingly inaccurate over time. For flows where transition to turbulence occurs naturally, as in laminar separation bubble flows, excessive dissipation may inhibit transition to turbulence and reattachement completely. Finally, any highly under- resolved simulation witha dissipative scheme will likely be particularly inaccurate. Under-resolution already implies that not all length scales of motion relevant to the problem will be captured. The amount of numerical dissipation is generally inversely proportional to the resolution – the more under-resolved, the higher the numerical dissipation. Combined, under-resolution and dissipative schemes may even preclude thedevelopment ofhighwave number content inprimaryquantities. SincemostLESmodelsarepredicatedontheabilitytopredict thecorrectsubgrid- scale dissipation rate based solely on coarser resolved scales, energy conservation in the numerical methods used is paramount as evidenced in Kravchenko & Moin (1997). Second-order methods often damp and deform high wave number con- tent of the primary variables. Subgrid-scale models are generally not designed to account for these effects, and given such flawed input are unlikely to compute the 29 correct subgrid-scale dissipation. As such, understanding, controlling, or remov- ing numerical dissipation at low resolution is of critical importance in predicting laminar separation bubble flows accurately and quickly. For these reasons, two Navier-Stokes solvers that employ high-order numerical methods were chosen to investigatethecapabilityofLEStoreducetheresolutionrequirements foraccurate laminar separation bubble predictions. 3.3.1 Center for Turbulence Research Code Developed by graduate students at the NASA Center for Turbulence Research (CTR) at Stanford, this code solves the compressible LES equations for a perfect gas Nagarajan et al. (2007). Henceforth this solver will be referred to as the CTR code. Derivatives are computed using a sixth-order finite difference approximation similar to a Padé scheme. The free parameters are chosen such that the resulting derivative approximations resolve higher waves numbers than otherwise possible Lele (1992). These high wave numbers are generally not well approximated or even severely damped in standard finite difference schemes of the same order. To maintain the spectral-like efficiency and high order of convergence of these derivative approximations, the horizontal directions are treated as periodic and sponge regions are used to simulate non-periodic flow. A split implicit-explicit time integration is used. For explicit time advancement in the freestream flow, a third-order Runge-Kutta scheme (RK3) is employed. A second-order A-stable implicit scheme is used near the wall to allow for larger stable integration time steps. Compact filtering as described in Lele (1992) is employed at each time step, bothinthefreestream andwall-normaldirections. Filteringisnecessary toremove aliasing errors, for overall stability and to ensure a smooth continuous solution at theinterface of theimplicit andexplicit computationaldomainsNagarajan(2004). 30 The numerical scheme is constructed on a structured curvilinear grid, and the variables are staggered in space. The freestream Mach number is chosen to be 0.2. Forcompressibleflowssubgrid-scalemodelsaredevelopedintermsofFavre-filtered quantities, as addressed for the first time in detail by Erlebacher et al. (1992), and in this work the dynamic Smagorinsky model is used in a form described in Sayadi & Moin (2012). Duetotheperiodicboundaryconditionsinthestreamwise direction, numerical spongesarenecessarytosimulatespatiallyevolvingflow. Thisspongeregionallows the flow to be recycled from outlet to inlet by forcing a return to the desired inlet boundary layer profile. Due to the code being compressible, a numerical sponge at the top boundary is also necessary to ensure sound and vortical waves are not reflectedbackintothecomputationaldomainMani(2012). Theinletspongeregion spansfromx = 0.03tox = 0.5,whereas theoutletspongestartsatx = 8andends at x = 9.2. The top sponge extends the domain from y = 1 to y = 1.8. Sponge regions account for one third of the total number of mesh points. These sponge layers relax the computed Navier-Stokes solution to the scale-similar compressible boundary layer case obtained a priori as a reference solution. From y = 1 to y = 1.4, the reference solution’s wall-normal velocity is changed to the suction profile specified in eq. (3.1). It is then smoothly brought back to its precomputed scale-similar value from y = 1.4 to y = 1.8. The sponge relaxation parameter increases from zero at y = 1, the end of the physical domain, and reaches its maximum close to the end of the sponge region at y = 1.8. The suction velocity is thus enforced indirectly through the influence of the forced solution above the top of the physical domain. 31 3.3.2 Spectral code in vorticity form A pseudo-spectral incompressible Navier-Stokes solver originally developed by Domaradzki & Metcalfe (1987) was modified and parallelized to perform LES of laminar separation bubble flow. The derivatives in the horizontal directions are approximatedusingFourierexpansions, whereastheverticalisapproximatedusing collocated Chebyshev polynomials. The flow variables are integrated in time using a fractional time step method described in Orszag & Kells (1980). It solves the filtered Navier-Stokes equations in rotational form: ∂¯ u i ∂t = [ǫ ijk ¯ u j ¯ ω k ]− ∂Π ∂x i +ν ∂ ∂x j ∂¯ u i ∂x j + ∂¯ u j ∂x i ! +F i − ∂τ sgs ij ∂x j , i = 1, 2, 3 (3.6) ∂¯ u i ∂x i = 0, (3.7) where ¯ ω i =ǫ ijk ∂¯ u k ∂x j is the vorticity, (3.8) Π =P/ρ+ 1 / 2 u i u i is the pressure head, (3.9) F i =f i (x,t), is a body force, and (3.10) τ sgs ij =u i u j − ¯ u i ¯ u j is the subgrid-scale stress tensor. (3.11) StrongalgebraicgridstretchingisusedintheverticaltoredistributetheChebyshev gridpointssuchthatapproximately 2 / 3 ofthepointslieinsidetheboundarylayerat its thickest while still resolving the inflow Blasius boundary layer. The streamwise velocityisthendecomposedintoabaseinflowprofileandvelocitydefect asfollows: u total = ˜ u+u, where ˜ u = [˜ u(y),0,0] is known and computed to be the Blasius 32 solution at the inlet x =x 0 . Hence,u is the departure away fromthe inflow profile at each location in space that is solved for at each time step. The non-linear term in u is computed in rotational form as in eq. (3.6). It is advancedintimeusingAdams-Bashforthscheme, whereastheadvectionduetothe base flow ˜ u is computed separately using Crank-Nicolson to obtainbetter stability and kinetic-energy-preserving characteristics as seen in eq. (3.12) and (3.13). ¯ u ∗ i = ¯ u n i + 3 2 ∆ t " N i +F i − ∂τ sgs ij ∂x j # n − 1 2 ∆ t " N i +F i − ∂τ sgs ij ∂x j # n−1 (3.12) ¯ u ∗∗ i = ¯ u ∗ i + 1 2 ˜ u i ∆ t ∂¯ u ∗∗ i ∂x j + ∂¯ u ∗ i ∂x j ! (3.13) where (3.14) N i =ǫ ijk ¯ u j ¯ ω k − ¯ u i ∂˜ u i ∂x j represent the non-linear term contributions in ¯ u i . (3.15) Thebaseflowadvection,pressureandviscousstepsareperformedusingtheFourier representation of the flow variables ¯ u i (x,y,z,t) = X |m|<M X |n|<N ˆ ¯ u i (k x ,y,k z ,t)exp(ik x x)exp(ik z z), (3.16) where k x = 2πm/L x and k z = 2πn/L z are the horizontal wavenumbers in the corresponding streamwise x and spanwise directions z, and L x and L z are the periodicity lengths. The base flow advection step described in eq. (3.13) then simplifies to eq. (3.17). ˆ ¯ u ∗∗ i = (1+b(y)k x ) (1−b(y)k x ) ˆ ¯ u ∗ i , where b(y) = 1 2 ˜ u(y)∆ t (3.17) 33 The pressure step consists of a correction to ensure the resulting flow field is divergence-free where the pressure is treated implicitly. Used together with the continuity equation (3.7), this results in a Poisson equation for vertical velocity with known Dirichlet boundary conditions as shown in eq. (3.18). The scheme of Orszag & Kells (1980) removes the need to impose explicit pressure boundary conditionsattheexpense ofa smallnumericalerror nearthewallwhich isO(ν∆ t). If a consistent boundary condition for the vertical velocity is used as in eq. (3.22), O(ν∆ t 3 2 ) may be reached Guermond et al. (2006). D 2 −k 2 ˆ ¯ v ∗∗∗ =−D[ik x ˆ ¯ u+ik z ˆ ¯ w] ∗∗ −k 2 ˆ ¯ v ∗∗ , (3.18) where (3.19) D = ∂ ∂y , (3.20) k 2 =k 2 x +k 2 z , (3.21) The above Poisson equation (3.18) is solved with Dirichlet boundary conditions ˆ ¯ v ∗∗∗ | b = ˆ ¯ v ∗∗ | b −ν∆ t h D(ik x (u ∗∗ −u n )+ik z (w ∗∗ −w n ))+k 2 (v ∗∗ −v n ) i | b (3.22) such that the pressure head and horizontal velocities can then be obtained alge- braically as follows ˆ Π =−[D ˆ ¯ v ∗∗∗ +ik x ˆ ¯ u ∗∗ +ik z ˆ ¯ w ∗∗ ]/(k 2 ∆ t) (3.23) ˆ ¯ u ∗∗∗ = ˆ ¯ u ∗∗ −ik x ∆ t ˆ Π (3.24) ˆ ¯ w ∗∗∗ = ˆ ¯ w ∗∗ −ik z ∆ t ˆ Π (3.25) 34 The viscous terms are then treated implicitly, which results in one Helmholtz equation for each direction in spectral space as shown in eq.(3.26). D 2 −k 2 − 2 ν∆ t ˆ ¯ u n+1 i =− 1 ν∆ t ( ˆ ¯ u ∗∗∗ i + ˜ u i ) (3.26) Numerical sponge regions are implemented in the streamwise direction using the fringe method formulation of Spalart & Watmuff (1993) to damp all turbulence and return the outflow to the desired inflow Blasius profile ˜ u(y) using the body force term (3.27): F i =λ(x)(˜ u i − ¯ u i ), (3.27) where λ(x) =κ f e −((x−x 0 )/σ l ) 2 +e −((x−xmax)/σr) 2 , tends to zero outside inflow and outflow, (3.28) κ f = 10 U 0 L f , is the forcing strength, (3.29) σ l = 1 10 L f , σ r = 3 8 L f , are the widths of sponge regions, (3.30) L f = 15 100 (x max −x 0 ) is the total size of the sponge regions. (3.31) This forcing formulation enables the simulation of a spatially evolving boundary layer while maintaining periodic boundary conditions and spectral accuracy. The sponge region at the inflow extends fromx = 0.25 tox = 0.45, whereas the sponge at the outflow spans x = 8.7 to x = 10. Fringe method forcing terms are active 35 over15%ofthetotaldomain. Theirinfluencedoesextendsomewhatbeyondwhere the terms are active, so to be safe the region x = 0.5 to x = 7.5 is considered to be the physically representative region of our computational domain. A ceiling suction boundary condition that matches Spalart’s and is designed to cause a 30% free stream flow deceleration are enforced. The resulting top streamwisevelocityiscomputedduringthepressurestepbyapplyingthecontinuity equationandimposingzerospanwisevelocity. Inpractice,theeffectivedeceleration is closer to 25% due to viscous effects as reported by Spalart & Strelets (2000). To ensure mass is conserved in the complete domain, the amount of fluid removed by suction from the top of the domain is injected in a narrow slot outside of the region of interest (from x = 8.5 to x = 9.5). The streamwise and vertical velocity boundary conditions are shown in figure 3.6. Blowing is applied at the wall to avoid any stability problems associated with disturbances near the top of thedomainwhere theresolutionisvery coarse. Blowing throughtheplateseverely limits the CFL number used in our simulations, increasing overall computational time. The upshot of this strict CFL restriction is that it removes any possibility of implicit filtering due to large time steps, a phenomenon observed in previous results using implicit schemes. 3.3.3 Spectral code in skew-symmetric form The vorticity form of the discretized Navier-Stokes equations was shown to be prone to aliasing errors and more sensitive to Gibbs phenomena by Kravchenko & Moin (1997), which can lead to persistent numerical oscillations at the smallest resolved scales in a spectral code due to the discretization method’s inherently negligible numerical dissipation. Concern that these sources of error 36 might be exacerbated at low resolution and affect the ability of explicit subgrid- scalemodelstoprovide thecorrect eddy viscosity ledtothedevelopment ofaspec- tral code with a higher order time integration scheme where the non-linear term is computed in skew-symmetric form. This form is known to be less prone to alias- ing errors and has better energy-conservation characteristics Kravchenko & Moin (1997). The solver is based a ona classic two-step pressure-correction scheme with a third-order backward difference formula (BDF3) outlined in Guermond et al. (2006), with boundary conditions used in the Zang-Hussaini algorithm outlined in Canuto et al. (2007) that should lead to O(∆ t 3 ) accuracy. The filtered Navier- Stokes equations are integrated in time using a combined non-linear viscous step 1 ∆ t   γ 0 ¯ u ∗ i − J−1 X q=0 α q ¯ u n−q i   −ν ∂ ∂x j ∂¯ u ∗ i ∂x j + ∂¯ u ∗ j ∂x i ! =− J−1 X q=0 β q N n−q i (3.32) with boundary conditions ¯ u ∗ | y=0 = ∆ t γ 0 " 2 ∂p ∂x n − ∂p ∂x n−1 # y=0 (3.33) D 2 ¯ u ∗ | y=Y = 0 (3.34) ¯ v ∗ | y=0 = ν∆ t γ 0 D 2 ¯ v n | y=0 (3.35) ¯ v ∗ | y=Y =V(x)+ ν∆ t γ 0 D 2 ¯ v n | y=Y (3.36) ¯ w ∗ | y=0 = ∆ t γ 0 " 2 ∂p ∂z n − ∂p ∂z n−1 # y=0 (3.37) D 2 ¯ w ∗ | y=Y = 0 (3.38) where γ 0 ,α q ,β q are backward difference formula of order J coefficients, and where N i = 1 2 ¯ u j ∂¯ u i ∂x j + ∂ ∂x j ( 1 2 ¯ u i ¯ u j +τ sgs ij )−F i . (3.39) 37 The resulting Helmholtz equations for ¯ u ∗ i are solved in spectral space. A pressure- correction step γ 0 ∆ t ¯ u n+1 i − ¯ u ∗ i + ∂p ∂x i n+1 = 0 (3.40) ∂¯ u i ∂x i = 0 (3.41) is then performed to obtain divergence free fields. To avoid applying explicit Neumann boundary conditions to the pressure, the continuity equation is used to derive a Poisson equation for ¯ v n+1 in spectral space that is equivalent to eq. (3.18), and subsequently obtaining the pressure and horizontal velocities algebraically as in equations (3.23), (3.24), and (3.25). The resulting velocity fields are then divergence free in the entire computational domain including any active sponge regions and satisfy the boundary conditions to O(ν∆ t 5 2 ) Guermond et al. (2006). This time integration scheme is stifly-stable for higher CFL numbers that the previous vorticity-based scheme Guermond et al. (2006). This property was only fully utilized in the initialization phase of simulations, whereas the time step was kept to the same order of magnitude as the vorticity-form spectral code during the time-averaging phase to ensure direct comparisons could be made. 3.4 Validation The CTR code was validated against analytical predictions for the growth rate of the inviscid mixing layer instability, the growth rate of Tollmien-Schlichting (T- S) waves in a non-parallel boundary layer, the level of distortion in representing a convecting Taylor vortex in time, and the pressure fluctuations due to sound scattered bya circular cylinder by Nagarajan(2004). Ithasbeen usedina number 38 of DNS and LES of bypass transition in Nagarajan et al. (2007) which reported good agreement with experiments. The pseudo-spectral incompressible code in vorticity form was previously vali- datedagainstlinearstabilitytheoryforparallelboundarylayers(u bl i ={U(z),0,0} and u bc i = 0) where errors of less 2% were reported Domaradzki & Metcalfe (1987) for the T-S wavelengths tested. This validation test was performed again after modifications were made to parallelize the code and replace obscure fast Fourier transformandlinearalgebrasubroutineswithcalltostandardopensourcelibraries like LAPACK and FFTW. Results are shown in table 3.1 were consistent and obtained numerical growth rates matching growth rates predicted by linear stabil- ity theory within 5% for a wide range of Reynolds numbers, domain sizes and T-S wavelengths, discarding 4outliers(outof25tests) wherethegrowthratepredicted was too close to zero. The same validation test was also performed for the spectral code in skew-symmetric form, which performed slightly better, with errors less than 2% with the exception of a few outliers as shown in table 3.2. 39 Table 3.1: Numerical (h ∂E ∂t i) vs theoretical (ω i ) growth rates of most unstable Orr-Sommerfeld modes in vorticity-form spectral code. Reynolds number α ∗ ω i h ∂E ∂t i h % error i 1805.00 0.500000 -0.115938E-01 -0.115268E-01 0.577699 1805.00 0.675000 -0.658913E-02 -0.653586E-02 0.808385 1805.00 0.850000 -0.226416E-02 -0.223016E-02 1.50171 1805.00 1.02500 -0.159557E-03 -0.154442E-03 3.20600 1805.00 1.20000 -0.108623E-02 -0.111583E-02 2.72447 2555.00 0.500000 -0.485721E-02 -0.481590E-02 0.850413 2555.00 0.675000 -0.970083E-03 -0.935548E-03 3.56002 2555.00 0.850000 0.185558E-02 0.187327E-02 0.953377 2555.00 1.02500 0.267053E-02 0.266628E-02 0.158898 2555.00 1.20000 0.920218E-03 0.888659E-03 3.42944 3305.00 0.500000 -0.211583E-02 -0.208443E-02 1.48432 3305.00 0.675000 0.921599E-03 0.945822E-03 2.62837 3305.00 0.850000 0.283207E-02 0.284275E-02 0.377046 3305.00 1.02500 0.293234E-02 0.292242E-02 0.338410 3305.00 1.20000 0.760838E-03 0.729025E-03 4.18128 4055.00 0.500000 -0.793901E-03 -0.769478E-03 3.07645 4055.00 0.675000 0.162942E-02 0.164762E-02 1.11693 4055.00 0.850000 0.294678E-02 0.295130E-02 0.153382 4055.00 1.02500 0.261108E-02 0.259848E-02 0.482493 4055.00 1.20000 0.192827E-03 0.164394E-03 14.7455 4805.00 0.500000 -0.915698E-04 -0.715104E-04 21.9061 4805.00 0.675000 0.188049E-02 0.189410E-02 0.723902 4805.00 0.850000 0.279080E-02 0.279260E-02 0.645979E-01 4805.00 1.02500 0.216841E-02 0.215651E-02 0.548716 4805.00 1.20000 -0.414861E-03 -0.435599E-03 4.99871 40 x Cp 0 0.2 0.4 0.6 0.8 1 -1.5 -1 -0.5 0 0.5 1 Jones et al. UDNS Dyn. Smag. WALE (a) NACA0012 airfoil at 5 degree angle of attack at Re c = 0.5×10 5 Castiglioni et al. (2014). 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x −C p (b) Flat plate with suction from the top applied at Re x = 10 5 . Figure 3.1: Time-average coefficients of pressure for airfoils is approximated using suction boundary conditions for a flat plate. 41 x y 1 2 3 4 5 6 7 0 0.1 0.2 0.3 0 0.5 1 (a) Contour plot of mean streamwise velocity. 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 (x=2) y 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 (x=3) y 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 (x=4) y 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 (x=6) y (b) Profiles of mean streamwise velocity before separation (x = 2), after separation (x = 3), attherecirculating vortex(x = 4), andin theturbulentboundarylayer (x = 6). Figure 3.2: The anatomy of a typical laminar separation bubble over a flat plate taken from large eddy simulation (LES) results using the truncated Navier-Stokes approach at 3% of DNS resolution. 42 Figure 3.3: Instantaneous iso-surfaces of spanwise vorticity on a flat plate due to suction from the top at Re x = 10 5 from spectral LES results with the σ-model at 3% of Spalart & Strelets (2000) DNS resolution. Figure 3.4: Physical domain, boundary and inlet conditions used to investigate laminar separation bubble flow over a flat plate. 43 (a) Upwind scheme solution (b) Crank-Nicolson scheme solution Figure 3.5: Solutions for the convection equation ∂ t T +u∂ x T = 0 after 40 time steps with CFL=0.8 reproduce from Fletcher (1991) 1 2 3 4 5 6 7 8 9 10 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.8 x velocity Figure3.6: SpectralLESboundaryconditions. Line: topverticalvelocityv(x,y = Y); dashed line: top streamwise velocity minus unity u(x,y =Y)−1; dash-dotted line: wall vertical velocity v(x,y = 0). 44 Table 3.2: Numerical (h ∂E ∂t i) vs theoretical (ω i ) growth rates of most unstable Orr-Sommerfeld modes in skew-symmetric spectral code. Reynolds number α ∗ ω i h ∂E ∂t i h % error i 1805.00 0.500000 -0.115938E-01 -0.115823E-01 0.994404E-01 1805.00 0.675000 -0.658913E-02 -0.658011E-02 0.136815 1805.00 0.850000 -0.226416E-02 -0.225150E-02 0.558928 1805.00 1.02500 -0.159557E-03 -0.143196E-03 10.2543 1805.00 1.20000 -0.108623E-02 -0.106467E-02 1.98480 2555.00 0.500000 -0.485721E-02 -0.485301E-02 0.864110E-01 2555.00 0.675000 -0.970083E-03 -0.963559E-03 0.672520 2555.00 0.850000 0.185558E-02 0.186313E-02 0.406872 2555.00 1.02500 0.267053E-02 0.268117E-02 0.398365 2555.00 1.20000 0.920218E-03 0.933770E-03 1.47269 3305.00 0.500000 -0.211583E-02 -0.211210E-02 0.176624 3305.00 0.675000 0.921599E-03 0.925912E-03 0.467992 3305.00 0.850000 0.283207E-02 0.283767E-02 0.197892 3305.00 1.02500 0.293234E-02 0.293782E-02 0.186724 3305.00 1.20000 0.760838E-03 0.768801E-03 1.04665 4055.00 0.500000 -0.793901E-03 -0.791426E-03 0.311795 4055.00 0.675000 0.162942E-02 0.163252E-02 0.190116 4055.00 0.850000 0.294678E-02 0.294896E-02 0.737562E-01 4055.00 1.02500 0.261108E-02 0.261351E-02 0.930939E-01 4055.00 1.20000 0.192827E-03 0.199857E-03 3.64557 4805.00 0.500000 -0.915698E-04 -0.896177E-04 2.13182 4805.00 0.675000 0.188049E-02 0.188215E-02 0.883248E-01 4805.00 0.850000 0.279080E-02 0.279187E-02 0.382952E-01 4805.00 1.02500 0.216841E-02 0.217091E-02 0.115136 4805.00 1.20000 -0.414861E-03 -0.403731E-03 2.68299 45 Chapter 4 Center for Turbulence Research Results Results for three cases computed using the CTR code are reported here: a benchmark DNS case (DNS), a wall-resolved LES with the dynamic Smagorin- sky model (LES), and an under-resolved DNS. Parameters for these simulations are summarized in Table 4.1. Both the DNS and LES were set up and run by a collaborator at CTR, Dr Taraneh Sayadi, whereas I set up and performed the UDNS cases. The DNS by Spalart & Strelet (2000) was initially intended to be the benchmark case. However, it was run using an incompressible spectral code with animposed vorticity-free boundary conditionat thetopboundary. These top boundary conditions could not be matched exactly in simulations with the CTR Spectral DNS CTR DNS CTR LES CTR UDNS N x 1022 1536 512 240 N y 120 300 140 90 N z 120 128 32 32 N total ×10 6 14.7 59.0 2.3 0.7 % of spectral DNS 100 401 15.6 4.7 % of CTR DNS 25 100 3.9 1.2 ∆ x + 20 9.7 26.4 57.0 ∆ y + at X = 7Y 1 0.5 1.0 1.6 ∆ z + 6.7 7.6 27.5 29.6 S effective 0.25 0.21 0.21 0.20 Table 4.1: Resolution and parameters for all cases run with the CTR code com- pared to the spectral DNS by Spalart & Strelets (2000). 46 code due to the fringe layer formulation used. The effective top boundary condi- tion is compared with the spectral DNS boundary condition of Spalart & Strelets (2000) in figure 4.1; the nominal deceleration parameter S effective = 0.21 is less than for the spectral DNS case (see Table 4.1). Results for additional simulations performed with a different numerical code are reported in Cadieux et al. (2012). For those additional cases, however, only LES and under-resolved DNS were per- formed because direct comparisons to the benchmark DNS by Spalart & Strelets (2000) were made. CTR DNS, LES, and UDNS were run until the separation bubble stabilized and turbulent flow was well established downstream of reattachment as illustrated in figure 4.2 and 4.3. Results were then averaged over multiple bubble “breathing” periods. All time-averaged results relating to pressure and friction coefficients obtained are in good qualitative agreement with the DNS benchmark (see fig- ure4.4a and 4.4b). The wall pressure coefficients C pw = (P w −P ∞ )/( 1 2 ρU 2 0 ) shown in figure 4.4a for the UDNS and LES cases are both in good quantitative agreement with the DNS benchmark with the exception of a slight difference in bubble length. The downward slope in C p in figure 4.4a after x = 5 indicates the existence of a slight favorable pressure gradient which extends to the end of the physical domain. This favorable pressure gradient is caused by blowing at the top boundary seen in figure 4.1. This presents a limitation in the applicability of results to the suction side of airfoils in MAVs and blades in turbo-machinery where such persistent favorable pressure gradients are seldom encountered downstream of the separation bubble Jones et al. (2010). Although weak, the favorable pressure gradient may also artificially improve agreement of LES and UDNS 47 1 2 3 4 5 6 7 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x Normalized velocity Figure 4.1: Normalized wall-normal velocity top boundary condition (V/U 0 at Y = 1): S & S 2000 Spalart & Strelets (2000) (circles), and UDNS (dashed line). Normalized mean streamwise difference from freestream velocity ((U −U 0 )/U 0 at Y = 1): UDNS (line). resultswiththeDNSbenchmarkbecauseofitseffectonthereattachment location. At resolutions onthe order of 1% of their respective benchmark DNS,andeven without models, all simulations predict the separation point seen in DNS bench- marks exactly. This can be observed in the first zero-crossing on the wall skin friction C f = µ ∂U ∂y | y=0 /( 1 2 ρU 2 0 ) plots in figure 4.4b. The UDNS predicts the same shape and maximum value of the peak negative skin friction as the benchmark DNS. Wall-resolved LES with dynamic Smagorinsky modeling performs slightly worse than the UDNS run, but still reaches within 15% of the DNS peak negative skin friction coefficient value. UDNS and LES predict the location of the reattach- ment point with less than 5% difference with the DNS. UDNS recovers benchmark 48 Figure 4.2: DNS snapshot of iso-surfaces of vorticity: Kelvin-Helmholtz rolls are visible over the separated shear layer leading to transition to turbulence and subsequent turbulent flow reattachment, closing of the separation bubble. DNS results almost exactly for the turbulent C f in the region downstream of the bubble whereas LES results underpredict the skin friction in that region. 4.1 Numerical Dissipation The good quantitative agreement between the no-model highly under-resolved DNS and benchmark DNS results suggests that the code used may belong to a category of implicit LES (ILES) where the numerical dissipation plays the role of subgrid-scalemodels. Asisevidentintheresultspresentedinfigures4.4aand4.4b, the addition of a subgrid-scale model, even when coupled with higher resolution, visibly worsens agreement with the DNS benchmark compared to the no-model case. Such behavior is expected for codes that already provide enough dissipation 49 x y 1 2 3 4 5 6 7 0 0.2 0.4 0 0.5 1 Figure 4.3: Contour plot of normalized average streamwise velocity U/U 0 from the UDNS case. Notice the laminar boundary layer growth followed by a clear separation bubble spanning from x = 2.8 to x≈ 4.6. through their numerics so that additional explicit subgrid-scale dissipation is not required. The code used has two primary sources of numerical dissipation: truncation error in derivative approximations, and explicit filtering. Since the code uses sixth order compact finite differences and is claimed to be conservative Nagarajan et al. (2003), focus is placed on quantifying the amount of effective numerical viscosity introducedbytheexplicitfilteringateachtimestepNagarajan(2004). Theexplicit high wave number filtering is based on the formulation of a sixth order compact filter Lele (1992). It is used to remove “spurious” and unstable high frequency oscillationsthatmaydevelopattheinterfaceoftheimplicitlyandexplicitlytreated regionsduetothecode’suseofhighorderfinitedifferencesNagarajan(2004). This is done to stabilize the code, and to ensure the implicit and explicit grid solutions matchattheir interface. It replaces theuseofartificialviscosity ornewer weighted essentially non-oscillatory(WENO) typeschemes, aswellaspenalty-type methods used to stabilize physical or numerical interfaces. Numerical dissipation from this filtering operation is quantified using two different methods. 50 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x C p w (a) Coefficient of pressure at the wall. 1 2 3 4 5 6 7 −2 −1 0 1 2 3 4 5 6 x 10 −3 x C f (b) Wall coefficient of friction. Figure 4.4: Time-averaged C p and C f . CTR DNS (circles), CTR LES with dynamic Smagorinsky model (line), and CTR UDNS (dashed line). 51 4.1.1 Estimating Numerical Dissipation due to Filtering The amount of effective viscosity the filtering operation imparts to the sim- ulation is estimated by comparing the energy decay rates of runs with filtering and without filtering. The number of time steps in such an analysis is limited to 10 to ensure that no numerical instabilities develop. First, two runs are per- formed using the same value of molecular viscosity, one with filtering and the other without. Second, the run without filtering is then repeated several times with larger values of the molecular viscosity until its energy decay curve matches that of the filtered case. The excess of the molecular viscosity in a run for which the best match is achieved provides an estimate of the effective viscosity that can be attributed to the filtering operation Diamessis et al. (2008). Since this approach was developed and validated for wakes and isotropic turbulence by Domaradzki et al. (2003); Domaradzki & Radhakrishnan (2005); Bogey & Bailly (2006); Diamessis et al. (2008), the procedure is carried through for a limited por- tion of the domain where a turbulent boundary layer is established after the sep- aration bubble. Incidentally, this is also where the dynamic Smagorinsky model should be particularly active and where LES predictions were the poorest match to DNS benchmark data. To confirm this, the same numerical procedure was per- formed with subgrid-scale modeling turned on. Average eddy viscosities calculated in the code were output to compare against assessed numerical viscosity. Total energy curves were obtained for the domain x = 5 to x = 7.5 and y = 0 to y = 0.5 from a set of restart files starting after 30,000 time steps. At such time, a statistically steady laminar separation bubble flow is well established. A typical total energy decay curve can be seen in figure 4.5a. Similar results were obtained using the same procedure with different starting iterations (e.g. 36,000, and 40,000). Since the domain selected for the analysis is only a portion of the 52 2 4 6 8 10 221,640.5 221,640.7 221,640.9 221,641.1 Iteration n, where n total = 30000 + n Total Energy (a) Total energy decay case 0 2 4 6 8 10 221,643.6 221,643.7 221,643.8 221,643.9 221,644 221,644.1 221,644.2 221,644.3 221,644.4 221,644.5 Iteration n, where n total = 36000 + n Total Energy (b) Total energy growth case Figure 4.5: Total energy in turbulent boundary layer following the laminar separa- tion bubble as a function of time. UDNS (squares), UDNS without filtering (line), UDNS without filtering and 18% larger molecular viscosity (dashed line), UDNS without filtering and 33% larger molecular viscosity (dash-dotted line). 53 entire computational domain and the solution is at a statistically steady state, the energy can either decay or grow over our short analysis times of a few time steps. However, thedescribedprocedureforestimatingtheeffectiveviscosityisapplicable also to cases where the total energy increases (see figure 4.5b). An over-arching conclusion fromall cases analyzed is that the effect of the filtering operation in the turbulent boundary layer downstream of reattachment is readily comparable to increasing the molecular viscosity by 18%. Given the non-dimensional molecular viscosity in the simulations is 2×10 −6 , the non-dimensional numerical viscosity is estimated at 0.35× 10 −6 . The average non-dimensional eddy viscosity in the same domain calculated using the dynamic Smagorinsky model from the same starting iteration (30,000) as in figure 4.5a was 0.71× 10 −6 . This corresponds to approximately 36% of molecular viscosity. The filtering operation’s numerical viscosity is thus of the same order of magnitude as the turbulent eddy viscosity provided by the dynamic Smagorinsky model, making it difficult to distinguish between effects of the model and that of numerical dissipation. 4.1.2 Quantifying Numerical Dissipation The amount of effective numerical viscosity is quantified more accurately using a recently developed technique Schranner et al. (2015). An equation for the evolu- tionofkineticenergyisobtainedbymultiplying theNavier-Stokes equationsbyu i . TheresidualE n ofthekinetic energy equationiscomputed over agiven subdomain from consecutive snapshots of the velocity fields at times t n−1 ,t n ,t n+1 . −E n = Z V ∂ρE kin ∂t dV + Z V ∂ρE kin u j ∂x j +u i ∂p ∂x i +ρνu i ∂τ ij ∂x j ! dV (4.1) 54 where E kin = 1 2 u i u i (4.2) τ ij = ∂u i ∂x j + ∂u j ∂x i − 2 3 ∂u k ∂x k δ ij (4.3) ∂ρE kin ∂t ≈ 1 2∆ t (ρE kin ) n+1 −(ρE kin ) n−1 (4.4) ∂ρE kin u j ∂x j +u i ∂p ∂x i +ρνu i ∂τ ij ∂x j ! ≈ ∂ρE kin u j ∂x j +u i ∂p ∂x i +ρνu i ∂τ ij ∂x j ! n HO (4.5) Eqn (4.1) represents the residual of the kinetic energy transport equation as the volumetric dissipation rate due to all numerical approximations made in the code Schranner et al. (2015). It is computed by comparing the rate of change of kinetic energyusingacentraldifferenceformula(4.4)wherethekineticenergyiscomputed directlyfromthevelocityfieldsattimet n−1 andt n+1 outputfromtheNavier-Stokes solver used – in this case the sixth-order finite difference CTR code. The exact transport terms are approximated using the velocity fields at time t n using the highest order discretization available, denoted by the subscript HO in eqn (4.5). Theresidualcomputedthiswaycanreadilybecomparedtothephysicaldissipation rate function E ν = Z V ρντ ij ∂u i ∂x j dV (4.6) computed at time t n using the same high order discretization, from which an effective relative numerical viscosity ν n ν = E n E ν (4.7) 55 canbecalculatedprovidedthephysicaldissipationratedoesnottendtozero. This analysis is carried out on the same 10 snapshots the previous estimates were com- puted from for direct comparison. The goalis to provide a more accurate estimate of numerical viscosity and better explain the role of filtering in this simulation by repeating the procedure for the 10 snapshots where there is no filtering active. The numerical dissipation rate and effective numerical viscosity of the CTR UDNS is computed using finite differences in the vertical and Fourier expansions in the horizontal for four regions: in the laminar boundary layer before separation for x ∈ [1.5,3.0], inside the separated region for x ∈ [3.0,4.5] immediately after reattachment for x ∈ [4.5,5.5], and further aft of the bubble in the turbulent boundary layer forx∈ [5.5,7.0]. These subdomains are also limited in the vertical to the maximum 99% boundary layer thickness in that region. These quantities along with the physical dissipation rate for comparison are shown in table 4.2 and 4.3. As expected, the effective numerical viscosity is nearly zero where the flow is laminar and well resolved - the numerical dissipation rate is much smaller than the physical dissipation rate due to the presence of the shear layer. The results for the filtered case and the unfiltered case are consistent in the laminar and separated regions, giving less than 4% difference. The filter can be said to have no effect on these regions, and the contribution of the truncation error of the scheme makes up most of the effective viscosity. In the region near the reattachment point, the flow experiences strong instabilities and transitions to turbulence. The Kelvin-Helmholtz rolls and secondary instabilities create sharp streamwise and vertical gradients. As the shear layer breaks down to turbulence, a number of small eddies are also created as shown in figure 3.3 and 4.2. These features form the basis of high wavenumber content which is then promptly removed by the 56 region <E ν >×10 −7 <E n >×10 −7 <ν n /ν >(%) laminar 8.40 1.20 5.98 separated 4.94 3.30 16.49 transitional 10.52 20.66 103.30 turbulent 21.43 6.38 31.91 Table 4.2: Time-averaged physical and numerical dissipation in 10 CTR UDNS snapshots. filter, generating an effective numerical viscosity as large as physical viscosity. In the turbulent region, all velocities oscillate randomly in space and time at the smallest scale resolved by the mesh. However, the vertical gradients are now better resolved due to the larger boundary layer thickness and smaller vertical fluctuations. Filtering explicitly damps the small scales in the streamwise and spanwise directions, removing energy from the simulation in the turbulent region. The numerical dissipation rate reaches 31.91% of the physical dissipation rate in the turbulent boundary layer. This confirms previous estimates that placed the effective numerical viscosity in the turbulent boundary layer on the same order of magnitude as the molecular viscosity. For comparison, the average ν sgs given by the dynamic Smagorinsky model in the turbulent region, where it should be at its maximum, is 36% of molecular viscosity ν. The numerical dissipation inherent in the CTR code due to filtering is thus approximately equal to the eddy viscosity predicted by the model. Adding explicit subgrid-scale modeling to the explicitly filtered simulation thus produces excessivedissipation,explainingwhyLESconsistentlypredictsbubblelength,peak negative skin friction, and turbulent skin friction less accurately than UDNS for this problem. When comparing the no-filter case to the filtered case, it becomes obvious that the filter provides the most dissipation in the transitional and turbulent regions, 57 region <E ν >×10 −7 <E n >×10 −7 <ν n /ν > (%) laminar (no filter) 8.40 1.18 5.90 separated (no filter) 4.95 2.53 12.63 transitional (no filter) 10.60 2.44 12.19 turbulent (no filter) 12.46 0.05 0.27 Table 4.3: Time-averaged physical and numerical dissipation in 10 CTR UDNS snapshots with filtering turned off. which fortuitously corresponds to where small scales are expected to be generated by the physics of the problem, and the same regions where a subgrid-scale model would be active. The no-filter case has nearly zero effective numerical viscosity in the turbulent region in comparison to the laminar and separated regions. This can be explained partly by the absence of sharp vertical gradients observed in the laminar and separated shear layer that may be under-resolved by the finite difference scheme. It could also indicate that the code is dispersive or nearly unstable in situations where high wavenumber content is present at the numerical interface between the explicit and implicit grids and nothing is done to smooth any mismatch. By filtering high wave number components in all directions regardless of whether they develop as a result of the physics or the numerics, the explicit filter effectively removes energy fromthesimulation andactsasasmall-scaleenergy dis- sipation mechanism. This filtering operation can thus be regarded as an implicit subgrid-scale model. For this laminar separation bubble flow and particular simu- lation parameters, it performed well and produced numerical viscosity comparable totheaverageeddyviscosity prescribedbythedynamicSmagorinskysubgrid-scale model in the correct regions. However, unlike an explicit subgrid-scale model, the results are unlikely to be repeatable with a different code or a different grid. 58 Chapter 5 Spectral Results I The skew-symmetric spectral large eddy simulation tool lsbflow detailed in section 3.3.3 is used to unambiguously assess the performance of different subgrid- scale models. The following simulations of a laminar separation bubble over a flat plate were performed at 3% and 1% of the corresponding DNS benchmark reso- lution: a no-model LES or under-resolved DNS to serve as a baseline (UDNS), an LES with the dynamic Smagorinsky model (DSM), another with the σ-model (SIG), and one using the truncated Navier-Stokes approach with automatic fil- tering (TNS). Each simulation is started from the same initial condition before the flow has transitioned to turbulence. All simulations were run until the sep- aration bubble stabilized and turbulent flow was well established downstream of reattachment beforestarting tocollect statistics. Time-averaging wasdoneover 40 non-dimensional time units t =t ∗U 0 Lx , well beyond the point of reaching averaging- period independence. Parameters for these simulations are summarized in Table 5.1. Results obtained sharply contrast previous LES results using sixth-order finite differences with explicit filtering where the no-model simulation provided the best agreement and adding the dynamic Smagorinsky model deteriorated agreement with the benchmark DNS Cadieux et al. (2014). Spectral under-resolved DNS results instead demonstrate poor agreement with the benchmark DNS at 3% of DNS resolution and even lose the ability to predict the separation point at 1%. Such disagreement is expected because the smallest relevant physical scales of the 59 Table 5.1: Resolution and parameters for LES cases and DNS by Spalart & Strelets (2000). DNS 3% LES 1% LES N x 1022 256 240 N y 120 56 32 N z 120 32 24 N total ×10 6 14.7 0.46 0.18 % of DNS 100 3.1 1.25 ∆ x + at x = 7 12 54 54 ∆ y + at x = 7 ≤1 0.2 0.27 ∆ z + at x = 7 6.3 26 34 laminarseparationbubbleproblemarenotcapturedwithsuchcoarsegridsandthe incompressible spectral code provides little dissipation due to approximations in the discretization of the Navier-Stokes equations in space and time. The no-model simulations instead suffer from a known problem common to all poorly resolved spectral simulations: slow accumulation of energy in the largest resolved wave numbers. This can in turn cause numerical oscillations and artificial backscatter which contaminates the larger scales, creating a feedback loop. Due to this phe- nomenon, UDNS predictions at 1% are no longer representative of the physics of the problem. The addition of a subgrid-scale model to provide dissipation successfully stops energy from accumulating in the smallest resolved scales and markedly improves results in all cases. Explicitly-modeled LES capture separation, transition, and turbulent reattachment accurately on both grids. The instantaneous spanwise vorticity and vertical streamwise velocity gradient for the truncated Navier-Stokes simulationat3%ofthebenchmarkDNSresolutionarepresentedinfigures5.1aand 5.1b to illustrate the laminar, transitional, and turbulent regions and the scale of theirstructures. Thecapabilitytopredictaccuratelyatlowcomputationalcostthe 60 x y 1 2 3 4 5 6 7 0 0.1 0.2 0.3 0.4 (a) Contours of spanwise vorticity for a slice along the x-y plane at z = 0. x z 1 2 3 4 5 6 7 0 0.2 0.4 0.6 (b) Contours of streamwise velocity vertical gradient at the wall ∂u ∂y | y=0 . Figure 5.1: Instantaneous LES results with TNS model at 3% of DNS resolution showing shear layer break up and flow reattachment. averageskinfrictionC f = 2µU y /(ρU 2 ∞ ),pressurecoefficientC p = 2(p−p ∞ )/(ρU 2 ∞ ), and the location of separation and reattachment is of particular interest to air- foil and blade designers. The performance of the LES performed is evaluated using the above-mentioned quantities by comparing to the spectral DNS data by Spalart & Strelets (2000). Other quantities such as boundary layer thicknesses, mean and RMS velocity profiles are also investigated to further characterize the performance of the different subgrid-scale models. 5.1 LES at 3% of DNS Resolution The coefficient of skin friction and coefficient of pressure predictions from each simulationisdisplayedinfigures5.2aand5.2balongwiththebenchmarkDNSdata 61 by Spalart & Strelets (2000). All 3% LES cases capture the location of the separa- tionpointexactly, visibleinthefirstzerocrossingoftheskinfrictioncoefficientC f . Beyond thepeaknegative frictionlocation, UDNSresults stronglydepart fromthe benchmark DNS. A larger separation bubble is observed. The delayed rise in the coefficient of pressure, as well as the skin friction going from negative to positive indicate a later reattachment point. The skin friction downstream ofreattachment is also poorlypredicted by the UDNS. In contrast, the simulations with anexplicit subgrid-scale model recover both C p and C f curves of the benchmark DNS almost exactly. A small delay in the onset of the negative peak skin friction is observed in all models. Similarly, a small discrepancy in nominal pressure in the turbulent boundary layer downstream of reattachment is observed. These discrepancies are likely related to the rapid coarsening of the grid due to the algebraic mapping. Figure 5.3 displays similar agreement for all LES with explicit models for the boundary layer thicknesses calculated using a pseudo-velocity U p based on the vertical integral of spanwise vorticity as follows U p (x,y)=− Z y 0 ∂v(x,y ′ ) ∂x − ∂u(x,y ′ ) ∂y ! dy ′ . (5.1) Thedisplacementthicknessisslightlyover-predictedintheseparatedregion,which could be due in part to the rapid coarsening of the grid at that vertical location. It is very difficult to distinguish between models in LES at 3% of DNS resolution, each recovering the pressure, skin friction and boundary layer thicknesses very accurately, discounting small grid-related discrepancies. Mean velocity profiles shown in figure 5.4 show minor disagreements with the benchmark DNS in terms of the vertical location of the separated shear layer as well as the amount of reverse flow present at x = 3.5. RMS fluctuation velocities 62 are found in good agreement for the dynamic Smagorinsky model and σ-model, whereas the truncated Navier-Stokes predictions agree perfectly at x = 2.5 but depart slightly from the benchmark at x = 3.5. 5.2 LES at 1% of DNS Resolution Furtherreducing thenumber ofcollocationpointsineachdirectiontoreach1% of the DNS resolution changes the performance of LES with explicit models only slightly. In contrast, the no-model simulation now converges to a flow field where the flow separates later, but reattaches and transitions quickly to turbulence with- out forming a proper laminar separation bubble. This is made evident in figures 5.5b and 5.5a by the later first C f zero crossing of the UDNS, its earlier second zero crossing, and its C p curve’s smooth rise without a characteristic plateau. Explicitly-modeled LES show consistent improvement in agreement with the benchmark DNS when compared to the UDNS for each key feature of the flow as shownintable5.2andinfigures5.5band5.5a. Thisimprovement isdirectlylinked with their subgrid-scale dissipation contribution −τ sgs ij ¯ S ij . This is shown in figure 5.6wherethemodels’subgrid-scaledissipationisplotted. Noticethatthestrongest contributionsforallmodelsarenearx = 3.5wheretheseparatedshearlayerbreaks upasisvisualizedinfigure3.3usingiso-surfacesofspanwisevorticity. Theeffective turbulent eddy viscosity reaches ν sgs ≈ 10ν there for both dynamic Smagorinsky and σ models. ν sgs then decreases until it is O(ν) at reattachment, and further decreases downstream in the turbulent boundary layer. As expected, the models arenotactiveinthelaminarandfreestreamregion,butdohavesomeimpactinside 63 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x C p (a) Pressure coefficients. 1 2 3 4 5 6 7 −4 −3 −2 −1 0 1 2 3 4 5 6 x 10 −3 x C f (b) Skin-friction coefficients. Figure 5.2: Pressure and friction coefficients for LES at 3% of DNS resolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; line: truncated Navier- Stokes. 64 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x δ * (a) Displacement thickness δ ∗ . 1 2 3 4 5 6 7 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x θ (b) Momentum thickness θ. Figure 5.3: Boundary layer thicknesses for LES results at 3% of DNS resolu- tion. Squares: Blasius solution; Circles: Spalart& Strelets(2000)DNS;diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorin- sky; line: truncated Navier-Stokes. 65 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (a) under-resolved DNS. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (b) dynamic Smagorinsky. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (c) σ-model. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (d) truncated Navier-Stokes. Figure 5.4: Profiles in LES at 3% of DNS resolution. Line: U; dashed line: u ′ /u ′ max . From left to right, plots spaced by 1.5: x = 2.5, x = 3.5. 66 the separated shear layer as early as x = 2.5, corresponding to the earliest non- zero turbulent fluctuations in v ′ as plotted in figure 5.10b. The average subgrid- scale dissipation provided by the truncated Navier-Stokes approach is estimated by summing the energy removed each time the automatic filtering operation is triggered as denoted by t m over one non-dimensional time unit t = t ∗U 0 Lx = N t ∆ t for which the filter is triggered M ≈N t /200 times: −hτ sgs ij ¯ S ij i≈ 1 N t ∆ t M X m=1 h E− ˜ E i t=tm (5.2) = 1 N t ∆ t M X m=1 " 3 X i=1 1 2 (u i − ˜ u i )(u i − ˜ u i ) # t=tm (5.3) This estimate is not representative of the dynamics of the truncated Navier-Stokes approach, but it would seem to indicate that its average dissipation contribution is one order of magnitude smaller than the dynamic and σ models. The maximum effective viscosity of the truncated Navier-Stokes approach is O(ν) and is located aft of x = 3.5. Its downstream location and negligible contributions before x = 3 may explain why truncated Navier-Stokes results better capture the onset of the peak negative skin friction. Truncated Navier-Stokes predicitons reach remarkable quantitative agreement with the benchmark DNS for C f and C p with an overall RMS error lower than the other two models. Its only departure is seen from x = 4.5 to x = 5 where it over-predicts the skin friction just downstream of reattachment by approximately 10%. LES with the σ-model obtains the second best agreement, followed closely by the dynamic Smagorinsky model. Both eddy viscosity models do not recover the turbulent skin friction coefficient as well as their 3% counterparts, especially inside the bubble from x = 3 to x = 3.5. This is explained in part by the more rapidly coarsening grid in the vertical, but also because of their dissipative effect 67 on all scales of motion (as opposed to only the smallest resolved scales). The dynamic Smagorinsky also under-predicts the peak negative skin friction and the reattachment location. Displacement and momentum thickness results seen infigures 5.7a and5.7bfor allexplicitlymodeledLESareconsistentwiththeircorrespondinghigherresolution results and show good agreement with the DNS benchmark. A clear departure is again observed for UDNS. Its displacement thickness confirms the lack of a typical separation bubble, whereas its momentum thickness indicates a quicker transition and relaxation to a near-equilibrium turbulent boundary layer. This is further evidenced in the agreement of the mean UDNS velocity in wall units with the log law of the wall near the outflow at x = 7 shown in figure 5.8. Perfect agreement with the log law of the wall was not reported in the DNS by Spalart & Strelets (2000). It is not seen in explicitly-modeled LES. The turbulent boundary layer they predict has yet to return to equilibrium at x = 7, more in line with the DNS benchmark. All simulations capture the viscous sublayer accurately. This is expected as the algebraic mapping is designed to be strong enough to resolve the inlet Blasius boundary layer and gives us a ∆ y + of 0.27 at x = 7. The performance of the models is further delineated in the mean velocity and RMS fluctuation velocity profiles shown in 5.9. RMS fluctuation velocities are computed from 40 flow fields taken approximately one non-dimensional time unit apart substracted from the mean flow which is averaged over a total of 40 non- dimensional time units. The mean and RMS velocity profiles are averaged in the spanwise direction. The departures of the dynamic Smagorinsky and σ models seen in the skin friction curve between x = 3 and x = 3.5 are mirrored in their mean velocity profiles and associated vertical derivatives. A stronger reverse flow 68 Table 5.2: Performance of subgrid-scale models at 1% of DNS resolution. UDNS DSM SIG TNS C p average % error 51 7.8 4.2 1.6 Peak −C f % error 52 17 3.4 3.9 Reattachment point % error 29 4.4 1.4 0.5 Turbulent C f average % error 2.5 2.9 4.7 5.1 Reattachment shape factor % error 9.1 10 8.4 2.2 Rel. CPU cost 1.0 2.5 1.4 1.2 thanfoundintheDNSispredicted bybothmodels, whereas thetruncatedNavier- Stokes results compare well with the benchmark. On the other hand, the RMS fluctuation velocity profiles of all explicitly modeled LES generally display satis- factory agreement with the benchmark. UDNS profiles confirm what is visible in its C f curve, that the flow transitions and reattaches between x = 3 and x = 3.5. The maximum RMS fluctuation of the streamwise and vertical velocity as a function of streamwise distance is shown in figures 5.10 along with the benchmark DNS data. Only the truncated Navier-Stokes approach recovers accurately the peak fluctuation velocities, whereas all models predict the fluctuation intensities downstream of reattachment beyond x = 5 satisfactorily. The accuracy with which each subgrid-scale model recovers key features of the flow of critical importance to airfoil and blade designers such as the pressure and skin friction, as well as the shape factor at reattachment are quantified by computing the relative percent error with respect to the benchmark DNS and shownintable5.2. Thecomputationalcostofeachsubgrid-scalemodelismeasured intermsofaveragenumberofiterationspersecond(performedonthesamenumber of MPI processes on the same machine). It is displayed in table 5.2 as the ratio of UDNS speed (iterations/second) to that of the modeled simulations to illustrate the extra computational time required for each model. The computational cost 69 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x C p (a) Pressure coefficients. 1 2 3 4 5 6 7 −4 −3 −2 −1 0 1 2 3 4 5 6 x 10 −3 x C f (b) Skin-friction coefficients. Figure5.5: Pressure andskinfrictioncoefficients forLESat1%ofDNSresolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; line: truncated Navier- Stokes. 70 of the dynamic procedure in the dynamic Smagorinsky model is higher than that of finding the singular values of the derivative tensor in the σ-model and nearly twice that of the periodic filtering in TNS. This is due to the communication cost incurred in the multiple test-filtering operations in the dynamic procedure. It should be noted that the dynamic Smagorinsky model’s relative computa- tional cost appears higher than previously reported by Piomelli for a comparable spectral code Piomelli (1999). The difference is not in model implementation, but in the algorithm chosen for the non-linear step. A three-step Runge-Kutta scheme (RK3) requires three times as many Fourier transforms as the scheme used here, and thus takes roughly three times as long to compute. Piomelli’s implementation uses RK3 but performs the dynamic procedure only once Piomelli (1999). This is computationally equivalent to performing the dynamic procedure every third iter- ationinthe solver presented here. It isthus reasonablethat therelative cost of the dynamic Smagorinsky model to the no-model simulation reported here is nearly three times the one reported by Piomelli. Such differences are to be expected in different Navier-Stokes solvers because the relative cost corresponds to the ratio of the total number of operations, including all the extra operations added by the model, to the number of operations required by the solver alone. Overall, the truncated Navier-Stokes approach and σ-model both outperform dynamic Smagorinsky models due to their better agreement with the DNS bench- mark and their lower computational cost. The accuracy of the TNS results at 1% of DNS resolution demonstrates that the use of periodic filtering alone can provide enough dissipation to act as a substitute for an explicit subgrid-scale model. It also lends weight to the applicability of the TNS approach beyond isotropic turbulence and channel flow previously tested Domaradzki et al. (2002); Tantikul & Domaradzki (2010, 2011). Its better performance is attributed to the 71 filtering operation only affecting the smallest resolved scales, whereas turbulent eddy viscosity models often have a dissipative effect on a wider range of scales. Excluding the drastic changes in the baseline no-model simulation (UDNS) results in going to the coarser mesh, LES results at 3% and 1% are consistent in their ability to predict C f , C p , and boundary layer thicknesses. Consistency of results for each respective model at these different resolutions indicates a healthy degree of grid independence, which is seldom obtained at such low resolutions. The σ model results match those of the dynamic Smagorinsky model but improve upon it in the separated shear layer where intense shear is present, vindicating its stated purpose Nicoud et al. (2011). 5.3 Spanwise Structures To investigate the importance of spanwise structures in laminar separation bubble flow, the velocity auto-correlation functions were computed with respect to spanwisedistanceatdifferentlocationsinthestreamwisedirectionataheightequal to the displacement thickness δ ∗ for the 1% LES case using the truncated Navier- Stokes approach with automatic filtering. They were then ensemble-averaged over 35 realizations separated by approximately one non-dimensional time unit. The plots of these auto-correlationfunctions in figure 5.11have a number of interesting features. figure 5.11a shows that the streamwise velocity fluctuations are strongly correlated from x = 2 to x = 4 where the effects of the suction boundary layer are strongest, whereas elsewhere in the domain, R uu tends to zero with distance r. Similarly, R vv shows a strong correlation nearest to the suction location at x = 3, but otherwise tends to zero. Spanwise fluctuations R w w all tend to zero more rapidly except for locations where the suction is strongest at x = 3 and where the 72 presence of Kelvin-Helmholtz rolls have been observed around x = 4. That the spanwise fluctuations do not tend to zero near the bubble likely indicates one of two things: a lack of spanwise resolution or that the spanwise structures are larger than the spanwise domain size. AtruncatedNavier-Stokessimulationwithautomaticfilteringwheretheextent of the domain and the number of collocation point in the spanwise direction were doubled was performed to test whether the size of the spanwise structures was not well captured inprevious simulations. The auto-correlationfunctions arecom- puted andensemble-averaged inthe sameway over 35realizations. Thefigure 5.12 displays each velocity correlation function at different streamwise locations. From comparison with figure 5.11, it is immediately visible that most velocity auto- correlations are more closely clustered around zero at r = 20 than their counter- parts on the smaller spanwise domain. One exception remains across all velocities at x = 3 where the boundary conditions impose a structure on the flow. Another exception is seen in the streamwise velocity correlation at x = 2 which is the location where separation occurs, again driven by the suction boundary condition. However, previous non-zero correlations at x = 4 are no longer visible, indicating that the turbulent structures are now smaller than the domain extent. For future very coarse LESoflaminar separationbubble over aflat plate, thespanwise extent should be larger than that used in Spalart & Strelets (2000). Although doubling the spanwise size may not be necessary, a minimum of L z = 1Y corresponding to L z = 0.5L b where L b is the bubble length should be used to ensure velocity auto-correlations fall to zero. Resolution also plays an important role in the size of spanwise structures. To visualize the difference, the velocity auto-correlation functions were computed and ensemble-averaged in the same fashion for a truncated Navier-Stokes simulation 73 with automatic filtering at 3% of the DNS resolution. The number of points in the spanwise direction is increased from the 1% LES-S by a factor of 4 3 , going from a ∆ z + = 34 to ∆ z + = 26. This results in an improvement in the tendency of all correlationstogotozerowithdistancerasseeninfigure5.13,buttheimprovement is not as dramatic as for a larger spanwise domain size. Since increasing resolution alone may not be enough, the earlier conclusion that the domain size should be increased in future coarse LES of laminar separation bubble flows over flat plates stands. 74 x y 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 0 0.1 0.2 0.3 0.4 −5 −4.5 −4 −3.5 −3 (a) Dynamic Smagorinsky model. x y 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 0 0.1 0.2 0.3 0.4 −5 −4.5 −4 −3.5 −3 (b) σ model. x y 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 0 0.1 0.2 0.3 0.4 −8 −7.5 −7 −6.5 (c) truncated Navier-Stokes. Figure5.6: Contoursofthelogarithmoftimeandspanwise-averaged subgrid-scale dissipation log 10 (h−τ sgs ij ¯ S ij i). Each contour represents a jump of 0.25. 75 1 2 3 4 5 6 7 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 x δ * (a) Displacement thickness δ ∗ . 1 2 3 4 5 6 7 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x θ (b) Momentum thickness θ. Figure 5.7: Boundary layer thicknesses for LES at 1% of DNS resolution. Squares: Blasiussolution;Circles: Spalart&Strelets(2000)DNS;diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; full line: truncated Navier-Stokes. 76 10 0 10 1 10 2 10 3 0 5 10 15 20 25 y+ u+ Figure 5.8: Turbulent boundary layer mean streamwise velocity in wall units at x = 7 for LES at 1% of DNS resolution. Squares: log law of the wall; circles: u + =y + ; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; full line: truncated Navier-Stokes. 77 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (a) under-resolved DNS. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (b) dynamic Smagorinsky. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (c) σ-model. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (d) truncated Navier-Stokes. Figure5.9: ProfilesinLESat1%ofDNSresolution(lines), andinDNS(symbols). Line, circles: U; dashed line, diamonds: u ′ /u ′ max . From left to right, plots spaced by 1.5: x = 2.5, x = 3.5. 78 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x u′ (a) u ′ max (x). 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x v′ (b) v ′ max (x). Figure 5.10: Maximum RMS fluctuation velocity in LES at 1% of DNS resolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dottedline: dynamicSmagorinsky; fullline: truncatedNavier- Stokes. 79 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz uu (a) R uu = 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz ww (b) R vv = <v ′ (x,y,z)v ′ (x,y,z+r)> <v ′ (x,y,z)v ′ (x,y,z)> 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz vv (c) R ww = <w ′ (x,y,z)w ′ (x,y,z+r)> <w ′ (x,y,z)w ′ (x,y,z)> Figure 5.11: Ensemble-averaged spanwise auto-correlation functions for TNS at 1% of DNS resolution. Blue squares: x = 1; light green circles: x = 2; red crosses: x = 3; black diamonds: x = 4; green triangles: x = 5; fuschia dash-dotted line: x = 6; teal dashed line: x = 7. 80 0 0.1 0.2 0.3 0.4 0.5 0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔ z uu (a) R uu = 0 0.1 0.2 0.3 0.4 0.5 0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔ z ww (b) R vv = <v ′ (x,y,z)v ′ (x,y,z+r)> <v ′ (x,y,z)v ′ (x,y,z)> 0 0.1 0.2 0.3 0.4 0.5 0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔ z vv (c) R ww = <w ′ (x,y,z)w ′ (x,y,z+r)> <w ′ (x,y,z)w ′ (x,y,z)> Figure 5.12: Ensemble-averaged spanwise auto-correlation functions for TNS with doubled spanwise domain size and spanwise collocation points. Blue squares: x = 1; light green circles: x = 2; red crosses: x = 3; black diamonds: x = 4; green triangles: x = 5; fuschia dash-dotted line: x = 6; teal dashed line: x = 7. 81 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz uu (a) R uu = 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz ww (b) R vv = <v ′ (x,y,z)v ′ (x,y,z+r)> <v ′ (x,y,z)v ′ (x,y,z)> 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r=nΔz vv (c) R ww = <w ′ (x,y,z)w ′ (x,y,z+r)> <w ′ (x,y,z)w ′ (x,y,z)> Figure 5.13: Ensemble-averaged spanwise auto-correlation functions for TNS at 3% of DNS resolution. Blue squares: x = 1; light green circles: x = 2; red crosses: x = 3; black diamonds: x = 4; green triangles: x = 5; fuschia dash-dotted line: x = 6; teal dashed line: x = 7. 82 Chapter 6 Spectral Results II Priortothecompletionofthedevelopmentofthespectrallargeeddysimulation tool lsbflow, another spectral solver described in section 3.3.2 was used to perform comparable simulations. LES at 3% of DNS resolution were performed with the truncated Navier-Stokes approach with a constant filtering period corresponding to 0.5% of one non-dimensional time unit and compared to a no-model simulation. LES at 1% of DNS resolution were done with three subgrid-scale models: the dynamicSmagorinskymodel(DSM),theσ model(SIG),andthetruncatedNavier- Stokes approach with a constant filtering period (TNS). An under-resolved DNS (UDNS) was also performed to serve as a base line comparison. Parameters for thesesimulationsaresummarizedintable6.1anddenotedbyLES-Vtodistinguish them from previous LES obtained with lsbflow. Simulations were run until the separation bubble stabilized and turbulent flow was well established downstream ofreattachment. Both3%and1%simulationresultsweretime-averagedovermore than 25 non-dimensional time units t =t ∗U 0 Lx . 6.1 LES-V at 3% of DNS Resolution At 3% resolution, no-model results display a larger separation bubble and do not recover the expected skin friction levels downstream of reattachment. This is consistent with previous spectral results, but the UDNS predictions are in even greater disagreement with the DNS benchmark in this case as is evidenced in 83 DNS 3% LES-V 1% LES-V N x 1022 240 224 N y 120 56 32 N z 120 32 24 N total ×10 6 14.7 0.43 0.17 % of DNS 100 2.9 1.1 ∆ x + at x = 7 12 54 58 ∆ y + at x = 7 <1 0.07 0.27 ∆ z + at x = 7 6.3 26 34 Table 6.1: Resolution and parameters for LES-V cases and DNS by Spalart & Strelets (2000). figures 6.1a, 6.1b, and 6.2. Discrepancies in C f and C p around x = 3 seen in the UDNS could be an indication of the no-model simulation reaching a different steady state with stronger and perhaps slower bubble ‘’breathing” cycles which were not smoothed out by the time averaging. Truncated Navier-Stokes results improve on the no-model simulation drasti- cally. The C p curve is well predicted with the exception of a later rise in pressure, indicating a slightly larger bubble. This is confirmed by its downstream second zero crossing oftheC f curve anditslarger displacement thickness intheseparated flow region visibile in figure 6.2a. TNS overpredicts the peak negative skin friction by 12%, predicts the reattachment point accurately with 2.4% relative error, and under-predicts the turbulent skin friction downstream of reattachment by 5.3%. The level of mean reverse flow it predicts at x = 3.5 is in good agreement with the benchmark DNS as shown in figure 6.3. Departures from the benchmark DNS in its mean velocity and RMS profiles at that location occur where the resolution coarsens rapidly and the boundary layer thickness is at its largest. Notice that the ∆ y + for the 3% LES is much smaller than for typical wall- resolved LES in which ∆ y + ≈ 1. It is in fact significantly smaller than required 84 for DNS ∆ y + ≈ 0.5. This means that many points are clustered very close to the wall and that very few points are used to resolve vertical locations beyond y = 0.2, which is where the separated shear layer is found to transition. Given thatpreviousinvestigatorshaveidentifiedresolvingthespearatedshearlayerbreak up and transition as the key to obtaining agreement, it is both surprising and encouraging that the truncated Navier-Stokes approach managed to recover the C p and C f curves this accurately. Conversely, this difference explains most of the inconsistencies found between the previous 3% LES results and the 3% LES-V results. 6.2 LES-V at 1% of DNS Resolution The grid is significantly coarsened in the vertical direction to reach 1% of the DNS resolution. The stretching parameter is also changed to ensure a larger ∆ y + which matches the previous 1% LES. These changes gravely deteriorate the no-model simulation predictions. As seen in figures 6.4a and 6.4b, the UDNS now predicts a later separation point but a quicker transition to turbulence and reattachment. This is largely consistent with previous 1% UDNS results. One exception is that the under-resolved DNS turbulent C f levels downstream of reat- tachment are under-predicted by nearly 40%. Another is that its mean velocity profileinwallunits doesnotagreewiththe loglawofthewallatx = 7 asexempli- fied in figure 6.6, whereas it did in previous results (see figure 5.8). The reason for this discrepancy between spectral UDNS results can only be numerical in nature, because all other parameters were kept constant. In any case, both UDNS suffer from excessive energy accumulation in the small scales which eventually starts to 85 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x C p (a) Pressure coefficients. 1 2 3 4 5 6 7 −6 −4 −2 0 2 4 6 x 10 −3 x C f (b) Skin-friction coefficients. Figure 6.1: Pressure and friction coefficients for LES-V at 3% of DNS resolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; line: truncated Navier- Stokes. 86 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x δ * (a) Displacement thickness δ ∗ . 1 2 3 4 5 6 7 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x θ (b) Momentum thickness θ. Figure 6.2: Boundary layer thicknesses for LES-V results at 3% of DNS resolu- tion. Squares: Blasius solution; Circles: Spalart& Strelets(2000)DNS;diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorin- sky; line: truncated Navier-Stokes. 87 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (a) under-resolved DNS. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (b) truncated Navier-Stokes. Figure 6.3: Profiles in LES-V at 3% of DNS resolution. Line: U; dashed line: u ′ /u ′ max . From left to right, plots spaced by 1.5: x = 2.5, x = 3.5. drive motion at the large scales and pollute results in ways that are not represen- tative of the underlying physics. By contrast, the explicitly-modeled LES predictions show improved agreement with the benchmark DNS. All models predict the pressure distribution and reat- tachment location very accurately. The truncated Navier-Stokes approach with a fixed filtering period performs overall better on the coarser grid than it did on the 3% grid. This is attributed to the better distribution of points away from the viscous sublayer toward the midplane so as to better resolve the separated shear layer at its thickest. It is the only model to the recover the turbulent skin friction downstream of reattachment accurately as seen in figure 6.4b. It has the lowest 88 RMS error for C f and C p . However, it overpredicts the peak negative skin friction as noted in table 6.2 and the size of the separated displacement thickness by 5.9% as shown in figure 6.5. It shows the best agreement for the mean velocity profile at x = 3.5 seen in figure 6.7, and for the streamwise distribution of the maximum RMS velocity u ′ max (x). The dynamic Smagorinsky model displays the second best agreement with DNS data. It recovers the skin friction, pressure, displacement and momentum thickness curves of the DNS benchmark almost exactly. How- ever, its predicted skin frictionlevels areconsistently below the DNScurve beyond the reattachment point. This is reflected in the larger maximum RMS velocity u ′ max (x) aft of x = 4.5 in figure 6.8a. The same behavior is observable in the σ model results for C f and u ′ max (x). Both the dynamic Smagorinsky and σ models display stronger reverse flow in figure 6.7 than TNS and the DNS benchmark at x = 3.5. Theσ modelsignificantly overshoots thepeaknegative skinfriction. This is surprising as its subgrid-scale dissipation contribution was previously shown to be almost identical to the dynamic Smagorinsky model’s as shown in figure 5.6. All models otherwise agree fairly well, but not exactly with the log law of the wall as displayed in figure 6.6, in agreement with Spalart & Strelets (2000) who did not report perfect agreement. Results obtained with this vorticity-form solver largely corroborateconclusions drawn from lsbflow results. Both sets of simulations demonstrate that accurate LES are possible at 1% of DNS resolution. The only notable discrepancy is the LES-V models’ consistent under-prediction of the pressure and skin friction down- streamofreattachment. Allotherparametershavingbeenkeptconstant,thismust be due to implementation differences. Namely, the exact way the top boundary conditions areenforcedforthestreamwise velocity differs between thetwo spectral solvers, and the vorticity-form is known to be more prone to aliasing errors and 89 Table6.2: Performanceofsubgrid-scalemodelsinLES-Vat1%ofDNSresolution. UDNS DSM SIG TNS C p average % error 46 3.2 3.4 4.0 Peak −C f % error 37 5.6 35 13 Reattachment point % error 20 1.1 2.0 1.0 Turbulent C f average % error 38 14 19 8.9 Reattachment shape factor % error 9.7 7.3 18 8.3 Gibbs oscillations when highly under-resolved, which could impede the models’ ability to predict the correct subgrid-scale dissipation. 90 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x C p (a) Pressure coefficients. 1 2 3 4 5 6 7 −4 −3 −2 −1 0 1 2 3 4 5 6 x 10 −3 x C f (b) Skin-friction coefficients. Figure 6.4: Pressure and friction coefficients for LES-V at 1% of DNS resolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; line: truncated Navier- Stokes. 91 1 2 3 4 5 6 7 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 x δ * (a) Displacement thickness δ ∗ . 1 2 3 4 5 6 7 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 x θ (b) Momentum thickness θ. Figure 6.5: Boundary layer thicknesses for LES-V at 1% of DNS resolution. Squares: Blasius solution; Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorin- sky; full line: truncated Navier-Stokes. 92 10 0 10 1 10 2 10 3 0 5 10 15 20 25 y+ u+ Figure 6.6: Turbulent boundary layer mean streamwise velocity in wall units at x = 7 for LES-V at 1% of DNS resolution. Squares: log law of the wall; circles: u + =y + ; diamonds: under-resolved DNS; dashed line: σ-model; dash-dotted line: dynamic Smagorinsky; full line: truncated Navier-Stokes. 93 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (a) under-resolved DNS. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (b) dynamic Smagorinsky. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (c) σ-model. 0 0.5 1 1.5 2 2.5 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 y (d) truncated Navier-Stokes. Figure 6.7: Profiles in LES-V at 1% of DNS resolution (lines), and in DNS (sym- bols). Line, circles: U; dashed line, diamonds: u ′ /u ′ max . From left to right, plots spaced by 1.5: x = 2.5, x = 3.5. 94 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x u′ (a) u ′ max (x). 1 2 3 4 5 6 7 0 0.05 0.1 0.15 0.2 0.25 x v′ (b) v ′ max (x). Figure6.8: MaximumRMSfluctuationvelocityinLES-Vat1%ofDNSresolution. Circles: Spalart & Strelets (2000) DNS; diamonds: under-resolved DNS; dashed line: σ-model; dash-dottedline: dynamicSmagorinsky; fullline: truncatedNavier- Stokes. 95 Chapter 7 Conclusions 7.1 Summary and Conclusions Large eddy simulations of a laminar separation bubble flow over a flat plate with a Reynolds number (10 5 ) and pressure distribution comparable to those seen on airfoils of small unmanned aircraft and on blades of low pressure turbines were performed using three high-order codes: a sixth-order compact finite dif- ference compressible large eddy simulation solver provided by the NASA Stanford Center for Turbulence Research (CTR code) and two Fourier-Chebyshev pseudo- spectral incompressible Navier-Stokes solvers. The performance of each simulation was evaluated against benchmark DNS data focusing on two quantities of crit- ical importance to airfoil and blade designers: time-averaged pressure (C p ) and skin friction (C f ) predictions used in lift and drag calculations. Boundary layer thicknesses, mean velocity profiles, and RMS fluctuation velocity profiles were also investigated. Agreement obtained with two spectral largeeddy simulation toolswith explicit subgrid-scale models demonstrates that accurate large eddy simulation of laminar separation bubble flows are attainable with as low as 1% of DNS resolution. The level of agreement obtained with explicitly-modeled LES is far superior to RANS results to date and has not been previously demonstrated for laminar separation bubble flows at such low resolution. Strong departures from the DNS bench- mark are observed in the no-model spectral simulations, indicating that explicit 96 subgrid-scale models are required in low-dissipation low-resolution cases. In con- trast, no-model predictions given by the CTR code were in better agreement with the CTR DNS due to its use of an explicit stability-enhancing filter providing enough dissipation without the use of a subgrid-scale model. In fact, the addition of an explicit subgrid-scale model in the CTR code dete- riorated agreement with the CTR DNS benchmark. Numerical dissipation due to filtering was first estimated by comparing energy decay rates in the turbulent boundary layer region downstream of reattachment to new simulations with no filtering active, but higher physical viscosity. The equivalent numerical viscosity was estimated to be of the same order of magnitude as physical viscosity and com- parableto theaverage eddy viscosity provided by thedynamic Smagorinsky model in the same region. Numerical dissipation due to filtering was then confirmed to play an important role equivalent to an implicit subgrid-scale model using a novel methodtoquantifyitbasedsolelyonpost-processingavailableconsecutivevelocity fields. To assess the performance of different subgrid-scale models without the effects of numerical dissipation seen in the CTR code results, a spectral solver with negli- gible numerical dissipation was used to perform large eddy simulation at two reso- lutions corresponding to 3% and 1% of the benchmark DNS by Spalart & Strelets (2000). Three explicit subgrid-scale models were compared to a no-model baseline case or under-resolved direct numerical simulation (UDNS) for these two reso- lutions. Simulations without an explicit subgrid-scale model (UDNS) failed to capture accurately the benchmark C f and C p at 3% resolution, and even failed to capture qualitative features of the expected flow on the coarsest mesh where its results were contaminated by excessive energy accumulating in the smallest resolved scales. The addition of explicit subgrid-scale models markedly improved 97 agreement with the benchmark DNS. At 3% of DNS resolution, all explicitly- modeled LES tested recover all important features of the flow almost exactly, with only minor discrepancies. The truncated Navier-Stokes approach performed the best on the coarsest grid, recovering the overall skin friction, pressure, boundary layerthicknesses, andvelocityprofilesoftheDNSbenchmarkmostaccurately. The σ model performed the second best, but like the dynamic Smagorinsky model, it did not capture the onset of the peak negative skin friction or the reattachment point as accurately at 1% of DNS resolution. However, both recovered the turbu- lent skin friction downstream of reattachment almost exactly. All results obtained with explicit subgrid-scale models are self-consistent as well as consistent amongst each other for the two resolutions and codes tested. The observed grid indepen- dence strengthens the validity of these subgrid-scale modeling approaches. Results obtained using a different spectral solver corroborate the above-mentioned con- clusions with only minor discrepancies attributed to differences in implementation and boundary conditions. That the physics of the problem were well captured at this low resolution is encouraging for applications with body-fitted meshes and higher operating Reynoldsnumbers. However, dissipative schemes maydeterioratelargeeddysimu- lation predictions as illustrated in the CTR largeeddy simulation results, meaning that higher order methods like finite element or discontinuous Galerkin may be required. Large eddy simulation coupled with such high-order methods should be able to converge to DNS results at resolutions two orders of magnitude lower than would otherwise be required for DNS, providing significant computational savings that could soon enable design optimization for applications that are dominated by the unsteady effects of separation and transition. For applications where such highordermethodsarenotavailableorarenotpractical, effortsshouldbemadeto 98 remove sources of numerical dissipation where possible. For example, the highest order central difference flux approximations available should be chosen instead of other more dissipative schemes like upwinding or essentially non-oscillatory schemes. Once a statistically-steady state is reached on a coarse grid, numerical dissipation should be quantified in different subvolumes using the method devel- oped by Schranner et al. (2015), especially in transitional and turbulent regions of the flow. If necessary, the resolution should be increased locally in these subvol- umes in order to decrease the effective numerical viscosity until it is an order of magnitude smaller than the subgrid-scale viscosity of an established subgrid-scale model. This approach should help the user converge to a solution more quickly. It isalso less costly thantypical grid-refinement studies because longaveraging times are not necessary. However, it does not provide the same invaluable convergence rate information. 7.2 Future Work The consistently good agreement obtained with the truncated Navier-Stokes at 1% of the DNS resolution may be evidence that this approach more accurately represents the underlying physics than the traditional approach to large eddy sim- ulation using the filtered Navier-Stokes equations with a modeled subgrid-scale stress tensor. The truncated Navier-Stokes approach thus warrants further study. Namely, it should be extended to and tested on other more flexible high order codes for both canonical flows as well as more complex geometries where DNS are available to assess its accuracy. To better understand its underlying dynamics, the evolution oftheenergy spectra andtransfer termsshouldbeplottedandcompared to a no-model simulation as well as a DNS of the same flow if available. 99 The simple method to quantify numerical dissipation proposed by Schranner et al. (2015) should be investigated as a tool to gauge whether explicit subgrid-scale models will yield worthwhile improvements in results for a given grid and a nominally under-resolved solution with no model for commercial package users. Moreover, this method could beused on-the-flyto help determine when and where subgrid-scale models should be active, or when filtering is necessary in the case of the truncated Navier-Stokes approach. Another application could be found in solution-adaptive mesh refinement in which areas where numerical dissipation is computed and found to be higher than an acceptable threshold are refined. The capability of large eddy simulation to succesfully predict time-averaged quantities at resolutions on the order of 1% of DNS should be tested on blades andairfoilsat ornear their operating Reynolds number inconfigurations forwhich DNS or experimental data is available. The same simulation setups should also be used to compare low resolution large eddy simulation results with different wall models active, and with hybrid RANS-LES formulations so their impact on the quality of results can be assessed. 100 Reference List AlamM,SandhamN(2000) Direct numericalsimulation of‘short’laminarsepara- tionbubbleswithturbulent reattachment. Journal of Fluid Mechanics410:1–28. Alving AE, Fernholz HH (1996) Turbulence measurements around a mild sep- aration bubble and downstream of reattachment. Journal of Fluid Mechan- ics 332:297–328. Anderson BW, Domaradzki JA (2012) A subgrid-scale model for large-eddy sim- ulation based on the physics of interscale energy transfer in turbulence. Physics of Fluids 24:065104. Basara B, Krajnovic S, Girimaji S, Pavlovic Z (2011) Near-wall formulation of the partially averaged navier stokes turbulence model. AIAA Journal49:2627–2636. Bogey C, Bailly C (2006) Large eddy simulations of transitional round jets: Influence of the Reynolds number on flow development and energy dissipation. Physics of Fluids 18:065101. Bose ST, Moin P (2014) A dynamic slip boundary condition for wall-modeled large-eddy simulation. Physics of Fluids (1994-present) 26:–. Burgmann S, Bücker C, Schröder W (2006) Scanning PIV measurements of a laminar separation bubble. Experiments in Fluids 41:319–326. Burgmann S, Dannemann J, Schröder W (2007) Time-resolved and volumetric PIV measurements of a transitional separation bubble on an SD7003 airfoil. Experiments in Fluids 44:609–622. Burgmann S, Schröder W (2008) Investigation of the vortex induced unsteadi- ness of a separation bubble via time-resolved and scanning PIV measurements. Experiments in Fluids 45:675–691. Cadieux F, Domaradzki JA, Sayadi T, Bose T (2014) DNS and LES of laminar separation bubbles at moderate Reynolds numbers. Journal of Fluids Engineer- ing 136. 101 Cadieux F, Domaradzki JA, Sayadi T, Bose T, Duchaine F (2012) DNS and LES of separated flows at moderate Reynolds numbers In Proceedings of the 2012 Summer Program, pp. 77–86. Center for Turbulence Research. Canuto C, Hussaini MY, Quarteroni A, Zang TA (2007) Spectral Methods : Fun- damentals in Single Domains Springer, Dordrecht. Castiglioni G (2015) Numerical Modeling of Separated Flows at Moderate Reynolds Numbers Appropriate for Turbine Blades and Unmanned Aero Vehi- cles Ph.D. diss., University of Southern California. Castiglioni G, Domaradzki J, Pasquariello V, Hickel S, Grilli M (2014) Numerical simulations of separated flows at moderate reynolds numbers appropriate for turbine blades and unmanned aero vehicles. International Journal of Heat and Fluid Flow 49:91 – 99 8th Symposium on Turbulence & Shear Flow Phenomena (TSFP8). Deck S (2012) Recent improvements in the zonal detached eddy simulation (zdes) formulation. Theoretical and Computational Fluid Dynamics 26:523–550. Deck S, ÃĽlie Weiss P, PamiÃĺs M, Garnier E (2011) Zonal detached eddy simu- lation of a spatially developing flat plate turbulent boundary layer. Computers & Fluids 48:1 – 15. DiamessisP,LinY,DomaradzkiJA(2008) Effectivenumericalviscosityinspectral multidomain penalty method-based simulations of localized turbulence. Journal of Computational Physics 227:8145–8164. Domaradzki JA, Metcalfe RW (1987) Stabilization of laminar boundary layers by compliant membranes. Physics of Fluids 30:695–705. Domaradzki JA, Radhakrishnan S (2005) Effective eddy viscosities in implicit modeling of decaying high Reynolds number turbulence with and without rota- tion. Fluid Dynamics Research 36:385–406. Domaradzki JA, Xiao Z, Smolarkiewicz PK (2003) Effective eddy viscosi- ties in implicit large eddy simulations of turbulent flows. Physics of Flu- ids 15:3890–3893. DomaradzkiJA,LohKC,YeePP(2002) Largeeddysimulationsusingthesubgrid- scale estimation model and truncated navier-stokes dynamics. Theoretical and Computational Fluid Dynamics 15:421–450. Dubois T, Domaradzki J, Honein A (2002) The subgrid-scale estimation model applied to large eddy simulations of compressible turbulence. Physics of Flu- ids 14:1781–1801. 102 Egorov Y, Menter F, Lechner R, Cokljat D (2010) The scale-adaptive simulation method for unsteady turbulent flow predictions. part 2: Application to complex flows. Flow, Turbulence and Combustion 85:139–165. Eisenbach S, Friedrich R (2008) Large-eddy simulation of flow separation on an airfoil at a high angle of attack and Re = 10 5 using Cartesian grids. Theor. Comp. Fluid Dyn. 22:213–225 10.1007/s00162-007-0072-z. Erlebacher G,HussainiM,Speziale C,ZangT(1992) Toward thelarge-eddysimu- lation of compressible turbulent flows. Journal of Fluid Mechanics 238:155–185. Fletcher CAJ (1991) Computational techniques for fluid dynamics,Vol. 1 Springer Verlag. Gaster M (1963) On the stability of parallel flows and the behaviour of separation bubbles Ph.D. diss., University of London. Girimaji S, Abdol-Hamid K (2005) chapter Partially-Averaged Navier Stokes ModelforTurbulence: ImplementationandValidation AerospaceSciencesMeet- ings. American Institute of Aeronautics and Astronautics 0. Guermond J, Minev P, Shen J (2006) An overview of projection methods for incompressible flows. Computer Methods in Applied Mechanics and Engineer- ing 195:6011 – 6045. Hadžić I, Hanjalić K (2000) Separation-induced transition to turbulence: Second- moment closure modelling. Flow, Turbulence and Combustion 63:153–173. Häggmark C (2000) Investigation of disturbances developing in a laminar separa- tion bubble flow Ph.D. diss., Royal Institute of Technology, Sweden. Hain R, Kähler CJ, Radespiel R (2009) Dynamics of laminar separation bubbles at low-Reynolds-number aerofoils. Journal of Fluid Mechanics 630:129–153. Herbst A,HenningsonD(2006) Theinfluence ofperiodicexcitation onaturbulent separation bubble. Flow, Turbulence and Combustion 76:1–21. HortonH (1968) Laminar separation bubbles in two and three dimensional incom- pressible flow Ph.D. diss., University of London. Howard R, Alam M, Sandham N (2000) Two-equation turbulence modelling of a transitional separation bubble. Flow, Turbulence and Combustion 63:175–191. HuH,YangZ,IgarashiH(2007)Aerodynamichysteresisofalow-Reynolds-number airfoil. Journal of Aircraft 44:2083–2085. 103 Jones LE, Sandberg RD, Sandham ND (2008) Direct numerical simulations of forced and unforced separation bubbles on an airfoil at incidence. Journal of Fluid Mechanics 602:175–207. Jones L, Sandberg R, Sandham N (2010) Stability and receptivity characteris- tics of a laminar separation bubble on an aerofoil. Journal of Fluid Mechan- ics 648:257–296. Kojima R, Nonomura T, Oyama A, Fujii K (2013) Large-eddy simulation of Low-Reynolds-Number flow over thick and thin NACA airfoils. Journal of Air- craft 50:187–196. Kravchenko A, Moin P (1997) On the effect of numerical errors in large eddy simulations of turbulent flows. Journal of Computational Physics 131:310 – 322. Lele SK (1992) Compact finite difference schemes with spectral like resolution. Journal of Computational Physics 103:16–42. LinJCM,L.PauleyL(1996) Low-Reynolds-number separationonanairfoil. AIAA Journal 34:1570–1577. Mani A (2012) Analysis and optimization of numerical sponge layers as a nonre- flective boundary treatment. Journal of Computational Physics 231:704 – 716. Marxen O, Henningson D (2011) The effect of small-amplitude convective distur- bances on thesize and bursting of a laminar separationbubble. Journal of Fluid Mechanics 671:1–33. Marxen O, Lang M, Wagner S (2003) A combined experimental/numerical study of unsteady phenomena in a laminar separation bubble. Flow, Turbulence and Combustion 71:133–146. MarxenO,RistU(2004) Directandlargeeddysimulationofthetransitionprocess in a laminar separation bubble In Proceedings of the 5th Workshop on Direct and Large-Eddy Simulation, pp. 231–240, Munich, Germany. Kluwer Academic Publishers. Marxen O, Rist U (2010) Mean flow deformation in a laminar separation bubble: separation and stability characteristics. Journal of Fluid Mechanics 660:37–54. McMullan W, Page G (2012) Towards large eddy simulation of gas turbine com- pressors. Progress in Aerospace Sciences 52:30 – 47 Applied Computational Aerodynamics and High Performance Computing in the {UK}. 104 Menter F, Egorov Y (2010) The scale-adaptive simulation method for unsteady turbulent flow predictions. part 1: Theory and model description. Flow, Turbu- lence and Combustion 85:113–138. Na Y, Moin P (1998) Direct numerical simulation of a separated turbulent bound- ary layer. Journal of Fluid Mechanics 374:379–405. NagarajanS (2004) Leading edge effects in bypass transition Ph.D. diss., Stanford University. Nagarajan S, Lele SK, Ferziger JH (2003) A robust high-order compact method for large eddy simulation. Journal of Computational Physics 191:392–419. NagarajanS,LeleSK,FerzigerJH(2007) Leading-edgeeffectsinbypasstransition. Journal of Fluid Mechanics 572:471–504. Nicoud F, Toda HB, Cabrit O, Bose S, Lee J (2011) Using singular values to build a subgrid-scale model for large eddy simulations. Physics of Fluids (1994- present) 23:–. Orszag SA, Kells LC (1980) Transition to turbulence in plane poiseuille and plane couette flow. Journal of Fluid Mechanics 96:159–205. Papanicolaou EL, Rodi W (1999) Computation of separated-flow transition using a two-layer model of turbulence. Journal of Turbomachinery 121:78–87. Piomelli U (1999) Large-eddy simulation: achievements and challenges. Progress in Aerospace Sciences 35:335 – 362. Piomelli U (2008) Wall-layer models for large-eddy simulations. Progress in Aerospace Sciences 44:437 – 446 Large Eddy Simulation - Current Capabilities and Areas of Needed Research. Piomelli U,BalarasE(2002) Wall-layer models forlarge-eddysimulations. Annual Review of Fluid Mechanics 34:349–374. Roberts S, Yaras M (2005) Large-eddy simulation of transition in a separation bubble. Journal of Fluids Engineering 128:232–238. SagautP(2006) Large Eddy Simulation for Incompressible Flows Springer, Berlin, Germany. Sagaut P, Deck S (2009) Large eddy simulation for aerodynamics: status and per- spectives. Philosophical Transactions of the Royal Society of London A: Mathe- matical, Physical and Engineering Sciences 367:2849–2860. 105 Sayadi T, Moin P (2012) Large eddy simulation of controlled transition to turbu- lence. Physics of Fluids 24:114103. Schranner FS, Domaradzki JA, Hickel S, Adams NA (2015) Assessing the numer- ical dissipation rate and viscosity in numerical simulations of fluid flows. Com- puters & Fluids 114:84 – 97. Shur ML, Spalart PR, Strelets MK, Travin AK (2008) A hybrid RANS-LES approachwithdelayed-DESandwall-modelled{LES}capabilities. International Journal of Heat and Fluid Flow 29:1638 – 1649. Skote M, Henningson D (2002) Direct numerical simulation of a separated turbu- lent boundary layer. Journal of Fluid Mechanics 471:107–136. Skote M, Henningson D, Henkes R (1998) Direct numerical simulation of self- similar turbulent boundary layers in adverse pressure gradients. Flow, Turbu- lence and Combustion 60:47–85. Sohn KH, Shyne RJ, DeWitt KJ (1998) Experimental investigation of boundary layer behavior in a simulated low pressure turbine National Aeronautics and Space Administration, Lewis Research Center. Spalart PR (2009) Detached-eddy simulation. Annual Review of Fluid Mechan- ics 41:181–202. SpalartP(2006) Topicsindetached-eddy simulation InGrothC,ZinggD,editors, Computational Fluid Dynamics 2004, pp. 3–12. Springer Berlin Heidelberg. Spalart P, Strelets M (2000) Mechanisms of transition and heat transfer in a separation bubble. Journal of Fluid Mechanics 403:329–349. Spalart P, Watmuff J (1993) Experimental and numerical study of a turbulent boundarylayerwithpressuregradients. JournalofFluidMechanics249:337–371. Spedding G, McArthur J (2010) Span efficiencies of wings at low Reynolds num- bers. Journal of Aircraft 47:120–128. Stolz S, Adams NA, Kleiser L (2001) An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows. Physics of Fluids 13:997–1015. Tantikul T, Domaradzki J (2010) Large eddy simulations using truncated Navier- Stokes equations with the automatic filtering criterion. Journal of Turbu- lence 11:1–24. 106 Tantikul T, Domaradzki J (2011) Large eddy simulations using truncated Navier- Stokes equations with the automatic filtering criterion: Reynolds stress and energy budgets. Journal of Turbulence 12:1–25. Wilcox DC (2006) Turbulent Modeling for CFD DCW Industries, La Canada, California. Wilson P, Pauley L (1998) Two-and three-dimensional large-eddy simulations of a transitional separation bubble. Physics of Fluids 10:2932–2940. WuX,MoinP(2010)Transitionalandturbulentboundarylayerwithheattransfer. Physics of Fluids 22:85–105. Xu T, Sullivan P, Paraschivoiu M (2010) Fast large-eddy simulation of low Reynolds number flows over a NACA0025. Journal of Aircraft 47:328–333. Yang Z, Voke P (2001) Large-eddy simulation of boundary-layer separation and transition at a change of surface curvature. Journal of Fluid Mechan- ics 439:305–333. Yarusevych S, Sullivan PE, Kawall JG (2006) Coherent structures in an airfoil boundarylayer andwake atlowReynoldsnumbers. Physics of Fluids18:044101. 107

Abstract (if available)

Abstract The flow over blades and airfoils at moderate angles of attack and Reynolds numbers ranging from 10⁴ to 10⁵ undergoes separation due to the adverse pressure gradient generated by surface curvature. In many cases, the separated shear layer then transitions to turbulence and reattaches, closing off a recirculation region—the laminar separation bubble. To avoid body‐fitted mesh generation problems and numerical issues, an equivalent problem for flow over a flat plate is formulated by imposing boundary conditions that lead to a pressure distribution and Reynolds number that are similar to those on airfoils. Spalart & Strelets (2000) tested a number of Reynolds‐averaged Navier‐Stokes (RANS) turbulence models for a laminar separation bubble flow over a flat plate. Although results with the Spalart‐Allmaras turbulence model were encouraging, none of the turbulence models tested reliably recovered time‐averaged direct numerical simulation (DNS) results. The purpose of this work is to assess whether large eddy simulation (LES) can more accurately and reliably recover DNS results using drastically reduced resolution—on the order of 1% of DNS resolution which is commonly achievable for LES of turbulent channel flows. LES of a laminar separation bubble flow over a flat plate are performed using a compressible sixth‐order finite‐difference code and two incompressible pseudo‐spectral Navier‐Stokes solvers at resolutions corresponding to approximately 3% and 1% of the chosen DNS benchmark by Spalart & Strelets (2000). The finite‐difference solver is found to be dissipative due to the use of a stability‐enhancing filter. Its numerical dissipation is quantified and found to be comparable to the average eddy viscosity of the dynamic Smagorinsky model, making it difficult to separate the effects of filtering versus those of explicit subgrid‐scale modeling. The negligible numerical dissipation of the pseudo‐spectral solvers allows an unambiguous assessment of the performance of subgrid‐scale models. Three explicit subgrid‐scale models—dynamic Smagorinsky, σ, and truncated Navier‐Stokes (TNS)—are compared to a no‐model simulation (under‐resolved DNS) and evaluated against the benchmark DNS data focusing on two quantities of critical importance to airfoil and blade designers: time‐averaged pressure and skin friction predictions used in lift and drag calculations. Results obtained with these explicit subgrid‐scale models confirm that accurate LES of laminar separation bubble flows are attainable with as low as 1% of DNS resolution, and the poor performance of the no‐model simulation underscores the necessity of subgrid‐scale modelling in coarse LES with low numerical dissipation.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

Numerical simulations of linearly stratified flow around a sphere

PDF

On the simulation of stratified turbulent flows

PDF

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

PDF

A Boltzmann model for tracking aerosols in turbulent flows

PDF

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

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

A novel sensor-based blended wall-modeling approach for large-eddy simulation of flows with relaminarization

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Laser manipulation of atomic and molecular flows

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Boundary layer and separation control on wings at low Reynolds numbers

PDF

Pattern generation in stratified wakes

PDF

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

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Modeling and simulation of complex recovery processes

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

Asset Metadata

Creator Cadieux, Francois (author)

Core Title Large eddy simulations of laminar separation bubble flows

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace and Mechanical Engineering (Computational Fluid and Solid Mechanics)

Publication Date 04/20/2015

Defense Date 02/27/2015

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

Tag computational fluid dynamics,direct numerical simulation,fluid mechanics,large eddy simulation,OAI-PMH Harvest,turbomachinery,UAV

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Domaradzki, Julian A. (committee chair), Becker, Thorsten W. (committee member), Lynett, Patrick J. (committee member), Redekopp, Larry G. (committee member), Spedding, Geoffrey R. (committee member)

Creator Email cadieux@usc.edu,francois.cadieux@gmail.com

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

Unique identifier UC11297777

Identifier etd-CadieuxFra-3340.pdf (filename),usctheses-c3-553536 (legacy record id)

Legacy Identifier etd-CadieuxFra-3340.pdf

Dmrecord 553536

Document Type Dissertation

Format application/pdf (imt)

Rights Cadieux, Francois

Type texts

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

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA