TableOfContents Turbulence Guide V4.01
User Manual:
Open the PDF directly: View PDF .
Page Count: 214
Download | ![]() |
Open PDF In Browser | View PDF |
Turbulence Models and Their Application to Complex Flows R. H. Nichols University of Alabama at Birmingham Revision 4.01 CONTENTS 1.0 2.0 3.0 4.0 5.0 Page Introduction 1.1 An Introduction to Turbulent Flow 1.2 Transition to Turbulent Flow 1.3 Statistical Concepts for Turbulent Flow Analysis 1.4 Turbulent Length and Time Scales 1.5 Vorticity and Strain Tensors 1.6 Classification of Turbulence Models References 1-1 1-4 1-9 1-13 1-14 1-15 1-17 Classical Turbulence Closure Methodologies 2.1 Reynolds Averaging 2.2 Favre Averaging 2.3 Boussinesq Approximation References 2-1 2-4 2-5 2-7 Turbulent Boundary Layer and Shear Layer Theory 3.1 Boundary Layer Theory 3.2 Shear Layer Theory References 3-1 3-7 3-10 Algebraic Turbulence Models 4.1 Inner Eddy Viscosity Model 4.2 Cebeci-Smith Model 4.3 Baldwin-Lomax Model 4.4 Wake and Jet Model 4.5 Algebraic Model Shortcomings 4.6 Grid Sensitivity for a Flat Plate with Adiabatic Walls 4.7 Grid Sensitivity for an Axisymmetric Bump 4.8 Grid Sensitivity for a Flat Plate with Heat Transfer 4.9 Grid Sensitivity for a Nozzle with Heat Transfer 4.10 Summary Baldwin-Lomax Application Hints References 4-1 4-1 4-2 4-4 4-4 4-6 4-10 4-14 4-18 4-21 4-22 4-22 One-Equation Turbulence Models 5.1 Theory 5.2 Spalart-Allmaras Model 5.3 Rotation/Streamline Curvature Corrections 5.4 Grid Sensitivity for a Flat Plate with Adiabatic Walls 5.5 Grid Sensitivity for an Axisymmetric Bump 5.6 Grid Sensitivity for a Flat Plate with Heat Transfer 5.7 Grid Sensitivity for a Nozzle with Heat Transfer 5.8 Summary Spalart-Allmaras Application Hints 5-1 5-1 5-3 5-6 5-10 5-13 5-17 5-20 5-20 6.0 7.0 8.0 9.0 References 5-20 Two-Equation Turbulence Models 6.1 Theory 6.2 Traditional k-ε Models 6.3 k-ω Models 6.4 SST Model 6.5 RNG Model 6.6 Numerical Implementation 6.7 Compressibility Correction for Shear Layers 6.8 Initializing Turbulence Values for a Given Profile Shape 6.9 Rotation and Curvature Correction 6.10 Grid Sensitivity for a Flat Plate with Adiabatic Walls 6.11 Grid Sensitivity for an Axisymmetric Bump 6.12 Grid Sensitivity for a Flat Plate with Heat Transfer 6.13 Grid Sensitivity for a Nozzle with Heat Transfer 6.14 Summary Two-Equation Model Application Hints References 6-1 6-3 6-6 6-7 6-9 6-10 6-10 6-13 6-15 6-16 6-19 6-22 6-26 6-29 6-30 6-30 Reynolds and Algebraic Stress Models 7.1 Reynolds Stress Models 7.2 Algebraic Stress Models 7.3 Grid Sensitivity for a Flat Plate with Adiabatic Walls 7.4 Grid Sensitivity for an Axisymmetric Bump RSM and ASM Model Application Hints References 7-1 7-5 7-7 7-12 7-18 7-18 Large Eddy Simulation 8.1 The Filtering Operation 8.2 Derivation of the LES Equations 8.3 Smagorinski Model 8.4 Dynamic Smagorinski Model 8.5 k-Equation Model 8.6 Inflow Turbulence Boundary Condition 8.7 Other LES References 8.8 Spatial Mixing Layer Example References 8-1 8-4 8-7 8-8 8-8 8-9 8-10 8-10 8-13 Hybrid RANS/LES Turbulence Models 9.1 Theory 9.2 Circular Cylinder 9.3 WICS Bay 9.4 Delayed Detached Eddy Simulation (DDES) 9.5 Summary Hybrid RANS/LES Application Hints 9-1 9-4 9-13 9-25 9-30 9-31 10.0 11.0 References 9-31 Wall Function Boundary Conditions 10.1 Theory 10.2 Grid Sensitivity for a Flat Plate with Adiabatic Walls 10.3 Grid Sensitivity for an Axisymmetric Bump 10.4 Grid Sensitivity for a Flat Plate with Heat Transfer 10.5 Grid Sensitivity for a Nozzle with Heat Transfer Wall Function Application Hints References 10-1 10-7 10-10 10-12 10-17 10-20 10-20 Boundary Layer Transition Simulation 11.1 Transition Models Based on Stability Theory 11.2 Transition Models with Specified Transition Onset 11.3 Transition Models with Onset Prediction Capability 11.4 Transition Models with Onset Prediction Capability References 11-2 11-2 11-6 11-10 11-14 1-1 1.0 Introduction Some basic knowledge of turbulence and an understanding of how turbulence models are developed can help provide insight into choosing and applying these models to obtain reasonable engineering simulations of turbulent flows. This effort is directed at production users of Computational Fluid Dynamics (CFD) and attempts to provide the basic information required to choose and use currently available turbulence modeling techniques. The limitations of each modeling technique and application tips are also provided. It is believed that the Navier-Stokes equations can be used to fully describe turbulent flows, but current limitations in computational horsepower have made the direct solution of the Navier-Stokes equations impractical for all but very simple flows at low Reynolds numbers. This is because current computers do not allow for the resolution of the wide range of length and time scales associated with turbulence. Many complex fluid dynamic applications are directed at determining time-averaged quantities, and hence it is desirable to find a means to obtain these mean quantities short of solving the full unsteady Navier-Stokes equations for all of the length and time scales associated with the turbulence. The quest for the ultimate turbulence model has been ongoing for nearly a century now. Early turbulence models were empirically derived algebraic relations. As computers developed and numerical simulation evolved differential equation based transport type turbulence models became the turbulence simulation methodology of choice. It should always be remembered that transport models are empirically calibrated. The use of transport type turbulence models has become standard practice for most engineering applications. Many current researchers are now solving the unsteady Navier-Stokes equations for large-scale, or grid realized, turbulence and modeling the smaller, or subgrid, turbulent scales that cannot be captured on the computational grid. Although the ultimate general-purpose turbulence model has yet to be developed, turbulence modeling has matured to the point that reasonably accurate results can be obtained for a wide range of engineering applications with the current class of computers. As computer technology continues to improve both the role and the form of turbulence models will continue to evolve. 1.1 An introduction to Turbulent Flow An understanding of what constitutes turbulent flow is required before proceeding to discuss turbulence modeling. Turbulence can be parameterized by several nondimensional quantities. The most often used is Reynolds number. Reynolds number represents the ratio of inertial forces to viscous forces. The viscous forces dominate at low Reynolds numbers and disturbances are damped rapidly. These disturbances begin to amplify as Reynolds number is increased and eventually transition into fully turbulent flows. Launder1 gives the following 1-2 definition for turbulent flow: “At moderate Reynolds numbers the restraining effects of viscosity are too weak to prevent small, random disturbances in a shear flow from amplifying. The disturbances grow, become non-linear and interact with neighboring disturbances. This mutual interaction leads to a tangling of vorticity filaments. Eventually the flow reaches a chaotic, nonrepeating form describable only in statistical terms. This is turbulent flow.” Turbulence is a three-dimensional unsteady viscous phenomenon that occurs at high Reynolds number. Turbulence is not a fluid property, but is a property of the flow itself. Turbulent flow can be highly nonlinear and is random in nature. Turbulent disturbances can be thought of as a series of three-dimensional eddies of different sizes that are in constant interaction with each other. This model of turbulence is Lagrangian in nature, with these turbulent flow structures being transported downstream by the mean flow. These structures exist for a limited amount of time before they are dissipated away by molecular viscosity. An example of the fluid structures that appear in a turbulent flow can be seen in the wake of a circular cylinder in Fig. 1.1a and the boundary layer in Fig. 1.1b. Figure 1.1a Turbulent structures in the wake of a circular cylinder. 1-3 Figure 1.1b Turbulent structures in a boundary layer. Turbulence is diffusive in nature. The rates of mass and momentum transfer are much higher in a turbulent flow than in a laminar flow. This diffusiveness can be a highly desirable property for enhancing the mixing of physical and/or chemical properties within a flow. The diffusive nature of turbulence also causes boundary layers and shear layers to become much thicker than their laminar counterparts, which is often an undesirable result of turbulent flow. A generic turbulent kinetic energy spectrum that indicates how the turbulent energy is partitioned among the various size eddies is shown in Fig 1.2. Inertial Range Dissipation Log(Energy) Energy Production Log(Wave number) Figure 1.2 Turbulent energy spectrum. The wave number is inversely proportional to the turbulent length scale. It can be seen in Fig 1.2 that most of the turbulent kinetic energy is in the large turbulent length scales (or low wave numbers). The large scale eddies take kinetic energy from the mean flow in the energy production region of the turbulent spectrum. Smaller eddies feed off of the larger energy producing eddies in the inertial range. Here the turbulence is essentially in equilibrium and the transfer of 1-4 energy by inertial forces is the dominant process. At the small-scale end of the turbulent spectrum the eddies become so small that viscous dissipation converts the kinetic energy into heat. As mentioned earlier, it is generally accepted that the Navier-Stokes equations can be used to directly simulate turbulent flows. This type of fluid modeling is called Direct Numerical Simulation (DNS). Direct numerical simulation of the turbulent energy cascade using the Navier-Stokes equations requires that the grid spacing be smaller than the smallest turbulent length scale. For a typical full-scale aircraft, the smallest length scale of the turbulence may be of the order of 10-6 times smaller than the aircraft reference length. More than 1018 grid points would be required to discretize this example. Since engineering problems of this size are impractical at this time, some model of the turbulence is required so that engineering answers may be obtained for large-scale fluid problems. For many engineering applications the length and time scales of the turbulence are much smaller than the length and time scales of the problem of interest. This tends to simplify the modeling of turbulence and remove the need to treat the unsteady aspects of the turbulent flow since they occur at time scales much smaller than those of interest. When the turbulent and problem scales become of the same size the modeling becomes more difficult. 1.2 Transition to Turbulent Flow Laminar-turbulent transition is an extraordinarily complicated process. The following is a quick summary of the process. More detail can be found in Schlichting2. The process begins by transforming external disturbances into internal instability oscillations in the boundary layer or shear layer. A laminar boundary layer or shear layer is susceptible to disturbances from both the free stream and the body surface. The ability of external disturbances to penetrate the boundary layer or shear layer and then be amplified is defined as the receptivity of the boundary layer. The free stream disturbances may include acoustic waves, particles in the flow, transported vorticity from upstream in the flow, and pressure, density, or temperature fluctuations. Surface disturbances include the roughness of the body or motion of the body. Many of these disturbances eventually get damped by viscosity. The transition process starts when a disturbance in the boundary layer or shear layer is no longer damped but gets amplified. This results in the formation Tollmien-Schlichting waves which are the first mode instability of the flow. While these waves travel downstream, three dimensional waves and vortices begin to develop. At certain points in the boundary layer, small irregularly shaped turbulent spots will occur and will be convected in a wedge shaped region as shown in the Figure.1.3. These spots, known as Emmons spots, appear at random locations on the plates at irregular time intervals. As these spots move downstream, they grow and eventually fuse with each other to encompass the entire boundary layer or shear layer. This process eventually leads to fully developed turbulent flow. 1-5 Laminar Transitional Turbulent Figure 1.3 Growth of turbulent spots in a flat plate boundary layer. The transition mechanism mentioned above is known as natural transition. Flows with strong free stream disturbances can “by-pass” the mechanism of formation of Tollmien-Schlichting waves and form turbulent spots directly from the influence of the free stream disturbances. This is called as by-pass turbulence3. For compressible flows, most notably flows with a Mach number of 2.2 or greater, Mack4,5 showed that multiple modes of instability exist. For flows between M=2.2 and M=4.5, the first mode of the instability, known as the viscous instability, is the most unstable. As the Mach number increases beyond 4.5, the second mode, known as the inviscid instability, becomes the most unstable mode. This second mode includes acoustical disturbances that are characterized by very large fluctuations in pressure and temperature. Boundary layer stability theory seeks to predict the receptivity of the boundary layer in terms the critical Reynolds number for a flow. The critical Reynolds number is defined as the Reynolds number at which the disturbances in a laminar boundary layer begin to amplify and the flow becomes turbulent. According to this theory, the fluid motion is decomposed into a time-averaged mean flow component and a fluctuating component. U = U + u′ (1.1) Here the overbar signifies a time averaged quantity and the prime is used for perturbation quantities. If the quantities in the Navier-Stokes equations are decomposed as in Eq. 1.1 and the mean flow quantities are subtracted out, the two-dimensional incompressible Navier-Stokes equations become 1-6 ∂u ′ ∂v′ + ∂x ∂y ∂u ′ +U ∂t ∂v′ +U ∂t dU 1 ∂p′ ∂u ′ + v′ + = ν (∇ 2u ′) dy ρ ∂x ∂x ∂v′ 1 ∂p′ + = ν∇ 2 v ′ ∂x ρ ∂y (1.2) Eq. 1.2 has three equations and three unknown disturbances u ′, v ′, and p ′ . If a periodic disturbance is assumed then the solution can be expressed as a Fourier series. This solution is linear, so solutions may be summed up. Note this is only valid for very small amplitudes for the disturbances. For two-dimensional flow, the stream-function can be defined as Ψ(x,y,t) = Φ(y)ei(α x-β t) Where α = (1.3) 2π is the wave number and λ is the wavelength of the disturbance λ and β = β r + β i i , where β r = 2πf ( f is the disturbance frequency) and β i i is the amplification factor for the disturbance. The ratio of β/α is c= β β r + βii = = cr + ci i α α (1.4) Here cr is the propagation velocity in the x-direction and is known as the phase velocity. The imaginary part, ci , is the degree of damping if ci is positive and is the degree of amplification if ci is negative. Substituting the expressions from Eq. 1.3 into Eq. 1.2 ∂φ = φ ′( y )e i (αx −βt ) ∂y ∂φ v′ = = −iαφ ′( y )e i (αx − βt ) ∂x u′ = (U − c )(φ ′′ − α 2φ ) − U ′′φ = − (1.5) i (φ ′′′′ − 2α 2φ ′′ + α 4φ ) α Re (1.6) Eq. 1.6 is called the Orr-Sommerfeld equation. The equation is shown in nondimensional form. The terms on the left side represent the inertia terms. The terms on the right side represent the viscous terms from the Navier-Stokes 1-7 equations. Here the Reynolds number (Re) is define as Re= Uδ ν where δ is the boundary layer thickness. Solving the Orr-Sommerfeld equations first requires a steady state fully laminar solution of the Navier-Stokes equations. The pressure and velocity from the laminar solution are used as inputs to the Orr-Sommerfeld equations. An example solution of the Orr-Somerfield is shown in Fig. 1.4. The vertical axis is the ratio of the shear layer thickness to the wavelength of the disturbance. The horizontal axis is the Reynolds number based on the shear layer thickness. Re Figure 1.4 Example solution to the Orr-Sommerfield equation. The curve shown in Fig. 1.4 is the neutral stability curve. Disturbances are amplified in the crosshatched region, while they are damped in the region outside the curve. The Reynolds number at which disturbances are first amplified is called the critical Reynolds number (Rcrit). For Reynolds numbers below Rcrit, disturbances are damped. As the Reynolds number increases above Rcrit, almost all disturbances are amplified and turbulent flow occurs. The effect of pressure gradient is shown in Fig. 1.5. The “0” curve represents no pressure gradient, the “-5” curve represents an adverse pressure gradient, and the “+4” curve represents a favorable pressure gradient. An adverse pressure gradient causes transition to occur more rapidly, while a favorable pressure gradient delays the onset of transition. 1-8 Re Figure 1.5 Effect of pressure gradient on boundary layer stability. Transition locations in boundary layer flows are generally characterized by two parameters: Rex and Reθ. Rex is the Reynolds number based on distance from the leading edge. Reθ is the Reynolds number based on the momentum thickness defined as θ= y =δ ρu ρ e ue y =0 ∫ ⎛ u⎞ ⎜⎜1 − ⎟⎟dy ⎝ ue ⎠ (1.7) Here δ is the boundary layer height and the subscript e denotes the edge of the boundary layer. Fig. 1.6 shows the transition process on the skin friction on a flat plate. The skin friction changes by almost an order of magnitude during the transition process. 1-9 Figure 1.6 Skin friction along a flat plate. 1.3 Statistical Concepts for Turbulent Flow Analysis Several statistical moments are used to describe turbulent flows. This section will review those most often associated with modeling turbulence. These statistics are usually performed on ensembles of discrete observations in both experimental and computational analysis of turbulent flows. The first moment is the mean or average and is defined as u= 1 n ∑ ui n i =1 (1.8) Here the average is taken of n discrete observations. The higher moments are defined relative to the mean. The second moment, or variance, is defined as 2 σ 2 = u ′2 = 1 n ∑ (ui − u ) n i =1 (1.9) The variance provides information on the spread of the data away from the mean. The third moment, or skewness, is given by S= u ′3 σ3 = 1 σ 3n n 3 ∑ (ui − u ) i =1 (1.10) 1-10 Skewness gives the amount of time the signal is above or below the mean. The fourth moment, or kurtosis, is defined by K= u ′4 σ4 = 1 σ 4n 4 n ∑ (u i =1 i − u) (1.11) Kurtosis gives the amount of time the signal is away from the mean. A simple example of these moments can be constructed using sine waves. Four cases are defined as ⎛ 2πt ⎞ Case 1 : sin⎜ ⎟ ⎝ 1000 ⎠ ⎛ 2πt ⎞ Case 2 : sin 2 ⎜ ⎟ ⎝ 1000 ⎠ ⎛ 2πt ⎞ ⎛ 2πt ⎞ Case 3 : sin⎜ ⎟ ⎟ + sin⎜ ⎝ 1000 ⎠ ⎝ 10,000 ⎠ ⎛ 2πt ⎞ 2 ⎛ 2πt ⎞ Case 4 : sin⎜ ⎟ + sin ⎜ ⎟ ⎝ 1000 ⎠ ⎝ 1000 ⎠ The four signals are shown in Fig. 1.7. Figure 1.7 Example signals for statistical analysis. The statistical moments for these four cases are shown in Table 1.1. (1.12) 1-11 Case 1 2 3 4 Mean 0.0 0.5 0.0 0.5 σ2 0.5 0.125 1.0 0.625 S 0.0 0.0 0.0 0.759 K 1.5 1.5 2.25 1.98 Table 1.1 Statistical moments for four test signals. Signals that are symmetric about the mean such as Case 1 and 3 have a zero mean and skewness. Case 4 is the only signal that produces a nonzero skewness. Time ensembles must be statistically stationary for the moments to have any meaning. A statistically stationary flow is defined as a flow in which statistical parameters do not depend on the interval in time used to evaluate them. An example pressure signal is given in Figure 1.8. Window 2 Window 1 Window 4 Window 3 Figure 1.8 Unsteady pressure measurement. It is very difficult to determine if the signal is stationary by visual inspection. Four overlapping data windows of 4096 samples each are used to assess whether this flow is stationary. The overlapped windowing technique is effective if at least two 1-12 of the windows are completely independent as is the case for Window 1 and Window 4. The number of windows and the number of samples in each window is generally determined by the number of total samples available in the signal. Fourier transforms are useful for transforming the signal between the time domain and the frequency domain. The Fourier transform of the four data windows and the average of the four transforms are shown in Figure 1.9. Figure 1.9 Fourier transform of pressure signal. Sound pressure level is proportional to the log of the variance of the pressure signal. The Fourier transforms for the individual windows show a large amount of disagreement. This is due to the discreet nature of the transform and the limited number of samples in each window. The statistical moments and the frequency of the maximum spectral peak for each of the data windows and the average of the windows are shown in Table 1.2. The maximum error in Table 1.2 is defined as the error between the individual data windows and the window average. This signal is reasonably stationary for the moments examined. It should be noted that demonstrating that a signal is stationary with respect to a given moment does not assure that the signal is stationary for higher moments. The windowing technique can be applied to both experimental and computational data. This technique can be used to assess “convergence” for unsteady computations6. 1-13 Average σ Pressure Window 1 1.2980 0.20586 S K 0.33207 2.3974 Frequency of Peak SPL (hz) 433.63 Window 2 1.2957 0.20787 0.36906 2.4891 461.61 Window 3 1.3003 0.20163 0.29161 2.3725 447.62 Window 4 1.2865 0.20218 0.29161 2.5564 461.61 Window 1.2952 0.02044 Average Max Error 0.66% 1.71% 0.32309 2.4539 451.12 14.23% 4.18% 3.88% Table 1.2 Statistics from pressure signal data windows. 1.4 Turbulent Length and Time Scales It is useful to quantify the range of length and time scales associated with a turbulent flow. These estimates of the maximum and minimum eddy length and time scales can be used in selecting the appropriate computational grid spacing and time step for a given problem. Turbulent scales are defined in terms of the turbulent kinetic energy (k), turbulent dissipation (ε), and the kinematic viscosity (ν). The turbulent kinetic energy and turbulent dissipation will be derived in Chap. 2. These two quantities are usually associated with the larger scales of turbulence. From dimensional analysis, the large scale, or turbulence producing, eddies (see Fig. 1.2) have length scales characterized by 3 L= k2 ε (1.13) The large scale eddies have time scales of the order of Τ= k ε (1.14) The smallest scale, or dissipative, eddies (the Kolmogorov scale) have length and time scales given by 1-14 1 ⎛ν 3 ⎞ 4 λ = ⎜⎜ ⎟⎟ ⎝ε ⎠ (1.15) 1 ⎛ν ⎞ 2 τ =⎜ ⎟ ⎝ε ⎠ (1.16) The ratio of the smallest to largest eddy length and time scales is then λ L τ Τ −3 = Rt 4 (1.17) −1 = Rt 2 (1.18) where Rt is the turbulent Reynolds number (k2/(νε)). These relationships indicate that the range of turbulent scales may span orders of magnitude for high Reynolds number flows. Since these turbulent length scales are much smaller than the physical scales (i.e. wing chord, channel height) associated with a flow of interest, it is easy to see that large amounts of grid points would be required to fully simulate a high Reynolds number turbulent flow. 1.5 Vorticity and Strain Tensors Two tensors that are often used in the development of turbulence models are the vorticity and strain tensors. The vorticity tensor is defined as Ω ij = 1 ⎛⎜ ∂ui ∂u j ⎞⎟ − 2 ⎜⎝ ∂x j ∂xi ⎟⎠ (1.19) The trace of the vorticity tensor (the sum of the diagonal elements) is equal to zero since Ω11=Ω22=Ω33=0. The largest magnitudes of vorticity typically occur near the wall in a boundary layer and in the core of a vortex. The strain tensor is defined as S ij = 1 ⎛⎜ ∂ui ∂u j + 2 ⎜⎝ ∂x j ∂xi ⎞ ⎟ ⎟ ⎠ (1.20) r The trace of the strain is ∇ • U , which is zero for steady incompressible flows. The largest magnitudes of strain also occur near the wall in a boundary layer. The magnitude of the strain is almost zero in the core of a vortex. For typical boundary layers and shear layers the magnitude of the vorticity is equal to the 1-15 magnitude of the strain since only one of the cross derivatives will dominate both tensors. 1.6 Classification of Turbulence Models Turbulence models can be classified by what turbulent scales they choose to model and what scales they choose to simulate by solving the unsteady NavierStokes equations on a computational grid. The traditional Reynolds Averaged Navier-Stokes (RANS) use a time averaging process to remove the necessity of simulating all of the scales of the turbulence spectrum. The RANS approach uses one length scale to characterize the entire turbulent spectrum, as shown in Fig 1.10. The use of a single length scale places a tremendous burden on the turbulence modeler since it can be difficult to find one length scale that is appropriate for all cases. When this can be accomplished the flow can be treated as a steady flow since all the unsteadiness is assumed to occur at scales below the computational grid size and are handled by the turbulence model. This allows for the use of time marching numerical algorithms with large numerical dissipation since the object of the calculation is to dissipate all wavelengths (whether they are of physical or numerical origin) of the unsteady flow during the convergence process. Numerical algorithms that converge to a steady state solution rapidly may not be suited for unsteady flow simulations because these algorithms contain large amounts of numerical dissipation and may overdamp the actual flow. The RANS approach is discussed in more detail in Chapter 2. RANS Log(Energy) Hybrid RANS/LES LES DNS Log(Wave number) Figure 1.10 Region of the turbulent energy spectrum modeled by different turbulence model approaches. 1-16 Direct numerical simulation (DNS)7 attempts to simulate all of the scales of turbulence and model nothing. The grid resolution and the maximum allowable time step for a DNS calculation must be small enough to capture the Kolmogorov scales (Eq. 1.15 and 1.16) of the turbulent flow. The inverse relationship of these scales with Reynolds number shown in Eq. 1.17 and 1.18 indicate that the grids must become increasingly fine as Reynolds number is increased. A DNS solution is inherently unsteady, so it must be run for long periods to assure that the resulting solution is statistically stationary and free from the initial conditions provided to start the computation. The numerical algorithm used in the DNS solution process must have very low numerical dissipation in order to allow all the wavelengths of the turbulent flow to naturally persist. The numerical dissipation must be much smaller than the molecular viscosity or it will manifest itself as an increase in the effective molecular viscosity, and subsequently as a decrease in the effective Reynolds number8. Most of the higher order numerical algorithms currently used in DNS are only valid on computational grids that contain uniform spacing. Grid stretching tends to reduce the order of the numerical algorithm and increase the numerical dissipation of the algorithm. The requirement of nearly uniform grids with small cell sizes leads to extremely large numerical grids even for small Reynolds numbers. DNS computations also require a model for the inflow free stream disturbances. For channel flows, the inflow turbulence can be obtained by recycling the outflow turbulence. To date DNS calculations have only been performed for low Reynolds numbers and simple geometries. Large Eddy Simulation (LES) is an attempt to move beyond DNS by modeling only the smallest turbulent scales in a problem9,10,11 as can be seen in Fig. 1.10. LES will be discussed in detail in Chapter 7. The smaller turbulent scales are more nearly isotropic, so they may be modeled with fairly simple turbulence models. Hence a simple turbulence model is used to simulate the “sub-grid” turbulence (turbulent scales that can not be realized on the computational grid) and the Navier-Stokes equations are solved for the remaining scales. This requires a spatial filtering of the turbulence spectrum. Like DNS, LES solutions are unsteady, and care must be taken to choose a numerical algorithm with low numerical dissipation or a decrease in effective Reynolds number will be seen in the calculation. LES solutions must also be run a large number of time steps to eliminate starting transients and to allow the solution to become statistically stationary. While LES allows for simulating much higher Reynolds numbers than DNS for the same computational resources, it is still not capable of simulating flight Reynolds numbers with reasonable engineering turn-around times. LES has been found useful for many lower Reynolds number applications such as combustion calculations where the chemical reactions are driven by the unsteady turbulent mixing of the fuel and the oxidizer. Recently a new class of turbulence models has arisen for unsteady high Reynolds flows. These turbulence models are called hybrid RANS/LES models and will be discussed in detail in Chapter 8. These hybrid turbulence models are extensions of the LES models in which modified RANS turbulence models are 1-17 used as sub-grid models for the Navier-Stokes equations. The goal for most applications is to use RANS to calculate the boundary layer where extremely small turbulent scales are present and use an LES like model to simulate the smaller turbulent scales away from the body. Hence the large scale turbulent structures away from the body are simulated by the unsteady Navier-Stokes equations. These models serve as a bridge between traditional RANS models and LES models as shown in Fig. 1.9. This allows for coarser grids than LES since the RANS-type sub-grid model is valid for non-isotropic scales of turbulence. These models require some sort of spatial filtering to determine the local value for the sub-grid turbulent viscosity. Hybrid RANS/LES models do not need the extremely low numerical dissipation required by DNS and LES and hence can be robust and easy to implement. Chapter 1 References: 1. Launder, B. E., “An Introduction to the Modeling of Turbulence,” VKI Lecture Series 1991-02, March 18-21, 1991. 2. Schlichting, H., Boundary-Layer Theory, McGraw-Hill, Inc., 1979. 3. Morkovin, M., "Bypass - Transition Research Issues and Philosophy," in Instabilities and Turbulence in Engineering Flows, Netherlands: Kluwer Academic Publishers, 1993, pp. 3-30. 4. Mack, L., "Linear Stability Theory and the Problem of Supersonic BoundaryLayer Transition," AIAA Journal, Vol. 13, Mar. 1975, pp. 278-289. 5. Mack, L., "Stability of Axisymmetric Boundary Layers on Sharp Cones at Hypersonic Mach Numbers," AIAA-87-1413, 1987. 6. Nichols, R., “Comparison of Hybrid RANS/LES Turbulence Models for a Circular Cylinder and a Cavity,” AIAA Journal, Vol. 44, No. 6, June 2006, pp. 1207-1219. 7. Mansour, N., “The Use of Direct Numerical Simulation Data in Turbulence Modeling,” AIAA-91-0221, Jan. 1991. 8. Moin, P., “Numerical and Physical Issues in Large Eddy Simulation of Turbulent Flow,” JSME International Journal, 1998, Series B, Vol. 41, No. 2, pp. 454-463. 9. Breuer, M., “Large Eddy Simulation of the Subcritical Flow past a Circular Cylinder: Numerical and Modeling Aspects,” International Journal of Numerical Methods in Fluids, 1998, Vol. 28, pp. 1281-1302. 10. Rodi, W., “Large Eddy Simulation of the Flows past Bluff Bodies,” in: Launder, B. and Sandham, N. (eds) Closure Strategies for Turbulent and Transitional Flows, Cambridge University Press, Cambridge, 2000. 11. Menon, S. “A Numerical Study of Secondary Fuel Injection Techniques for Active Control of Combustion Instability in a Ramjet,” AIAA-92-0777, January 1992. 2-1 2.0 Classical Turbulence Closure Methodologies The direct numerical solution of the Navier-Stokes equations for high Reynolds numbers flow requires a grid fine enough to capture all of the turbulent length scales and a time step small enough to capture the entire turbulent spectrum. The computational requirements for such calculations are well beyond the current state of computers. Obviously, some approximations to the NavierStokes equations are required to allow practical engineering solutions to high Reynolds number flow problems. Two averaging techniques used to simplify the Navier-Stokes equations are Reynolds and Favre averaging. Both of these approaches will be discussed in this chapter. Both of these approaches results in time-averaged correlation terms that must be modeled to close the equation set. Two other methods of closure are large eddy simulation (LES) and hybrid RANS/LES modeling. Both of these techniques use spatial filtering to decide what turbulent scales should be modeled and what scales can be simulated. LES will be discussed in detail in Chapter 8. The hybrid approach will be discussed in Chapter 9. 2.1 Reynolds Averaging Osborne Reynolds1 suggested that for a statistically stationary flow (see Chapter 1) the variables in the Navier-Stokes equations could be decomposed into timeaveraged and turbulent-fluctuation terms. The velocity can then be written as U ( x, t ) = U ( x ) + u ′( x, t ) (2.1) where the average velocity U is defined as U = 1 T Udt 2T ∫−T (2.2) and the time average of the fluctuating velocity u ′ is u′ = 1 2T ∫ T −T u ′dt = 1 2T ∫ (U − U ) dt = 0 T −T (2.3) Note that the time scale of the integration T must be much larger than the turbulent time scales discussed in Chapter 1 or the flow will not be statistically stationary. In the Reynolds averaging process the unsteady behavior of the turbulent flow is replaced by steady-state correlations terms that must then be calibrated with turbulent flow measurements. This has the advantage of significantly reducing 2-2 the computer time requirements for obtaining engineering quality values of steady state variables. It is often difficult to obtain empirical correlations that are valid for a wide range of turbulent flows, and hence the Reynolds Averaged Navier-Stokes (RANS) turbulence models must often be tuned to a particular class of flows. Assuming that fluctuations in density ( ρ ′ ), viscosity ( μ ′ ), and thermal conductivity ( κ T′ ) are negligible (which is often not the case for high speed and reacting flows) and applying Reynolds averaging to the Navier-Stokes equations yields ∂ρ ∂ (ρ U j ) = 0 + ∂t ∂x j (2.4) ( ) ∂ρ U i ∂ (ρU iU j ) + ∂ ρ ui' u 'j = − ∂P + ∂σ ij + ∂t ∂x j ∂x j ∂xi ∂x j ∂ρ E ∂ (ρU j H ) + ∂ ⎛⎜ ρ u 'j e' + ρ U i ui' u 'j + 1 ρ ui' ui' u 'j ⎞⎟ = + ∂t ∂x j ∂x j ⎝ 2 ⎠ ( ) ∂ (U iσ ij ) + ∂ ui'σ ij' − ∂ ∂x j ∂x j ∂x j ⎛ ∂T ⎞ ⎜κ T ⎟ ⎜ ∂x ⎟ j ⎠ ⎝ (2.5) (2.6) where ⎛ ⎝ 1 3 ⎞ ⎠ σ ij = 2μ ⎜ S ij − S kk δ ij ⎟ S ij = 1 ⎛⎜ ∂U i ∂U j ⎞⎟ + ∂xi ⎟⎠ 2 ⎜⎝ ∂x j 1 E = e + U iU i + k 2 H =E+ k= P ρ 1 ' ' ui ui 2 (2.7) (2.8) (2.9) (2.10) (2.11) 2-3 Note that the turbulent kinetic energy (k) has been included in the definition of total energy ( E ). The Reynolds averaging has introduced new unknown correlations k, u i' u 'j , u i' e ' , u i' u i' u 'j , and u i'σ ij' . When the velocity correlation u i' u 'j is multiplied by the density ρ it represents the transport of momentum due to the fluctuating (i.e. turbulent) motion of the fluid. The term ρ ui' u 'j is the transport of xi momentum in the direction of xj (or vice versa). This term acts as a stress on the fluid and is therefore called turbulent or Reynolds stress. In most regions of a turbulent flow, the turbulent stresses are ⎛ ∂U i ⎞ ⎟ . The Reynolds stress is a much larger than their laminar counter parts ⎜ν ⎜ ∂x ⎟ j ⎠ ⎝ symmetric tensor. The trace of the stress tensor is two times the turbulent kinetic energy, and is coordinate system independent. The off diagonal terms of the Reynolds stresses are coordinate system dependent. The velocity-energy correlation ρ ui' e ' can be thought of as a turbulent heat flux. The ρ ui' ui' u 'j triple correlation represents turbulent transport or diffusion. The u i'σ ij' term is the turbulent dissipation. Closure of the equation set requires correlations be developed for these new terms. These terms are often difficult to measure experimentally. The turbulent dissipation, u i'σ ij' , includes correlations of the perturbed velocity and spatial derivatives of the perturbed velocity. This quantity cannot be measured experimentally. Transport equations can be developed for the turbulent stresses and turbulent dissipation to close the system. This approach to closure is called Reynolds Stress Modeling (RSM) and will be discussed in Chapter 6. A simplified approach to modeling the turbulent stresses in which the RSM equations are replaced with algebraic relationships is called Algebraic Stress Modeling (ASM) and will also be discussed in Chapter 6. The Reynolds stresses must meet certain constraints to be physically plausible. Schumann2 introduced the realizability constraint. This constraint specifies that all the component energies of the turbulent kinetic energy (the diagonal terms of the Reynolds stress tensor) remain non-negative and all off-diagonal components of the Reynolds stress tensor satisfy Schwartz’s inequality. This can be written as ui' ui' ≥ 0 (2.12) 2 ui' ui' u 'j u 'j − ui' u 'j ≥ 0 (2.13) 2-4 ( ) 2 u1' u1' ⎛⎜ u2' u2' u3' u3' − u2' u3' ⎞⎟ − u1' u2' u1' u2' u3' u3' − u2' u3' u1' u3' + ⎠ ⎝ ( ) (2.14) u u u u u u −u u u u ≥0 ' ' 1 3 ' ' 1 2 ' ' 2 3 ' ' 2 2 ' ' 1 3 2.2 Favre Averaging Reynolds averaging for compressible flows produces a very complicated set of equations because the fluctuations of density, pressure, and temperature must also be accounted for. Favre3 averaging or Favre3 filtering has been used to simplify the equations set. Favre averaging is a density weighted averaging process in space. The Favre filter can be defined for any variable as ~ U ( x, t ) = U ( x ) + u ′′( x, t ) (2.15) where ∫ ρ (t )U (t )dt = ρU ρ ∫ ρ (t )dt T ~ U= −T (2.16) T −T Note that the over bar signifies a time averaged quantity, and the tilde represents a Favre filtered quantity. Other auxiliary relations include ρu ′′ = 0 (2.17) ~ ~ ρU = ρ U = ρU Thus, using Favre averaging, ~ ~ ρ ( x, t )U ( x, t ) = ρU + ρu ′′ = ρ U (2.18) The equivalent expression using Reynolds averaging would yield ρ ( x, t )U ( x, t ) = ρ U + ρ ′u ′ (2.19) Note that under these definitions the Reynolds average perturbation u ′ = 0 , but the Favre averaged perturbation u ′′ ≠ 0 . Applying this filtering operation to the Navier-Stokes equations, and assuming that the filtering commutes with the derivative operation, yields 2-5 ∂ρ ∂ ~ + ρU i = 0 ∂t ∂xi ( ) (2.20) ~ ∂ρ U i ∂ ∂ ∂P ∂ ~ ~ ~ + + ρ U iU j + ρui′′u ′j′ = − σ ji + σ ′ji′ ∂t ∂x j ∂x j ∂xi ∂x j ( ) ( ) ~ ∂ρ E ∂ ∂ ~ ~ + ρU j H + ∂t ∂x j xj ( ( ) ) 1 ~ ⎛ ′′ ⎞ ⎜ ρu j e + U i ρui′′u ′′ j + ρu ′j′ui′′ui′′ ⎟ = 2 ⎝ ⎠ ~ ∂ ~~ ∂ ∂ ⎛⎜ ∂T ∂T ′′ ⎞⎟ ~ κ U iσ ij + ui′′σ ij + U i σ ij′′ − + κT T ∂x j ∂x j ∂x j ⎜⎝ ∂x j ∂x j ⎟⎠ ( ) (2.21) ( ) (2.22) For almost all flows σ~ij >> σ ij′′ , and hence the σ ij′′ terms are usually neglected. ~ ∂ 2T ∂ 2 T ′′ ∂ 2 T ′′ Similarly, for most flows , and hence the terms are usually >> ∂x 2j ∂x 2j ∂x 2j neglected. The Favre averaged equations are similar in appearance to the incompressible Reynolds averaged Navier-Stokes equations, but the correlations are not the same. Favre averaging should only be used in flows with large density fluctuations such as in combustion or hypersonic flow. In other cases the simpler Reynolds averaging should be used. 2.3 Boussinesq Approximation Boussinesq4 suggested that the turbulent stresses can be treated in an analogous form to the viscous stresses in laminar flows. The Reynolds stresses are modeled as ⎛ ∂U ∂U j − ρ u i' u 'j = μ t ⎜ i + ⎜ ∂x ∂xi ⎝ j ⎞ 2 ⎟ − ρkδ ij ⎟ 3 ⎠ (2.23) The Boussinesq approximation states that the Reynolds stresses are ⎛ ∂U ∂U j ⎞ ⎟ . The eddy viscosity proportional to the local mean flow strain rate ⎜ i + ⎟ ⎜ ∂x x ∂ i ⎠ ⎝ j (μt) is the proportionality factor for the Boussinesq approximation. This is sometimes called a local equilibrium assumption since the turbulent stresses are assumed to be proportional to the mean flow strain rates. The last term in Eq. 2.23 was added so that the normal stresses would sum to the turbulent kinetic energy (k). The Boussinesq approximation reduces the turbulence modeling process from finding the six turbulent stress components as discussed above to determining an appropriate value for the eddy viscosity. 2-6 The eddy viscosity is a property of the flow field and not a property of the fluid. The eddy viscosity can vary significantly form point to point within the flow, and also from flow to flow. The eddy viscosity concept was developed assuming a relationship or analogy exists between molecular motion and turbulent motion. Although turbulent eddies can be thought of as discrete lumps of fluid that collide and exchange energy, they are not rigid bodies and their mean free paths are not small compared to the eddy size, and hence turbulence does not really satisfy the constraints of kinetic gas theory. In spite of the theoretical weakness of the eddy viscosity concept, it does produce reasonable results for a large number of flows. The eddy viscosity concept has been successful in the prediction of flows in which the shear stress ( τ = ρ ui' u 'j ) is the turbulent stress of greatest importance. This class of flows includes simple shear layers and attached boundary layers. The eddy viscosity is introduced as a scalar in Eq. 2.23 so that the eddy viscosity is the same for all stress components. The assumption of isotropic eddy viscosity is a simplification that limits the performance of this class of turbulence models in complex flows. In a direct analogy to the turbulent momentum transport defined by the Boussinesq approximation (Eq. 2.23), the turbulent heat or mass transport is often assumed to be proportional to the gradient of the transported quantity − uiϑ = Γ ∂ϑ ∂xi (2.24) where Γ is the turbulent diffusivity of heat or mass. Like the eddy viscosity, Γ is not a property of the fluid but is a function of the flow. Reynolds analogy between heat or mass transport and momentum transport suggests that Γ is proportional to the eddy viscosity Γ= νt σt (2.25) where σt is the turbulent Prandtl or Schmidt number. Experiments have shown that the σt varies little throughout a flow field and hence it is usually taken as a constant. Using the Boussinesq approximation and the assumption of a constant turbulent Prandtl number, the Reynolds averaged Navier-Stokes equations (Eq. 2.4, 2.5, and 2.6) can be written as ∂ρ ∂ (ρU j ) = 0 + ∂t ∂x j (2.26) 2-7 ∂ρU i ∂ (ρU iU j ) = − ∂P + ∂σ ij + ∂xi ∂x j ∂x j ∂t ∂ ∂ρE (ρU j H ) = ∂ (U iσ ij ) − ∂ + ∂x j ∂x j ∂x j ∂t (2.27) ⎛⎛ μ ⎞ ⎞ ⎜ ⎜ + μ t ⎟ ∂T ⎟ ⎜ ⎜ Pr Pr ⎟ ∂x ⎟ t ⎠ j ⎠ ⎝⎝ (2.28) where ⎛ ⎝ ⎞ ⎠ σ ij = 2(μ + μ t )⎜ Sij − S kk δ ij ⎟ S ij = 1 ⎛⎜ ∂U i ∂U j + ∂xi 2 ⎜⎝ ∂x j 1 3 ⎞ ⎟ ⎟ ⎠ 1 E = e + U iU i + k 2 H =E+ P ρ (2.29) (2.30) (2.31) (2.32) where k is the turbulent kinetic energy. The turbulent kinetic energy is often ignored in the definition of total energy (Eq. 2.31). Under the assumption of a constant turbulent Prandtl number (Prt), the equation set is closed once the eddy viscosity is defined. Some common approaches to defining eddy viscosity models are discussed in Chapters 3, 4, and 5. Chapter 2 References: 1. Reynolds, O., “On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion,” Philosophical Transactions of the Royal Society of London, Series A, Vol. 186, pp. 123-164, 1895. 2. Schumann, U., “Realizability of Reynolds Stress Turbulence Models,” Physics of Fluids, Vol. 20, pp. 721-725, 1977. 3. Favre, A., “Equations Fondanmentales des Fluids a Masse Variable en Ecoulements Turbulents,” Pages 24-78 of Favre, A., Kovasznay, L., Dumas, R., Gaviglio, J., and Coantic, M. (eds.), La Turbulence en Micanique de Fluids, CNRS, Paris, 1962. English edition: The Mechanics of Turbulence, Gordon and Breach, New York, 1964. 4. Boussinesq, J., Theorie de l’ecoulement Turbillonnant et Tumultueux des Liquides Dans les lits Rectilignes a Grande Section, Gaulthier-Villars, 1897. 3-1 3.0 Turbulent Boundary Layer and Shear Layer Theory A brief review of simple theory for turbulent boundary layers and shear layers will aid in the understanding of the development and calibration of turbulence models. Most of these relationships are based on empirical correlations of the shape of the velocity profile across the boundary layer or shear layer. This shape can be used in conjunction with the Navier-Stokes equations to back out relationships for the turbulent stresses or the eddy viscosity. 3.1 Boundary Layer Theory Turbulent boundary layers are usually described in terms of several nondimensional parameters. The boundary layer thickness, δ, is the distance from the wall at which viscous effects become negligible and represents the edge of the boundary layer. Two integral parameters across the velocity profile are the displacement thickness δ * and the momentum thickness θ. ∞ ⎛ ρu ⎞ ⎟ dy δ = ∫ ⎜⎜1 − ρ e ue ⎟⎠ 0⎝ * ∞ θ =∫ 0 ρu ρ e ue ⎛ u⎞ ⎜⎜1 − ⎟⎟ dy ⎝ ue ⎠ (3.1) (3.2) The integration is performed normal to the wall and the subscript e is used to denote the edge of the boundary layer at y=δ. The displacement thickness is a measure of the increased thickness of a body due to the velocity defect of the boundary layer. The momentum thickness is the distance that, when multiplied by the square of the free-stream velocity, equals the integral of the momentum defect across the boundary layer. If a simple power law velocity profile as shown in Eq. 3.3 is assumed and the flow is incompressible, the relationships in Eqs. 3.4-3.7 can be obtained from approximations of the Navier-Stokes equations. u ⎛ y⎞ =⎜ ⎟ ue ⎝ δ ⎠ δ x δ* x = 1/ 7 0.371 Re1x/ 5 = 0.046 Re1x/ 5 (3.3) (3.4) (3.5) 3-2 θ x = 0.036 Re1x/ 5 Cf = 1 τw 0.0577 = 2 Re1x/ 5 2 ρ e ue (3.6) (3.7) Cf is called the skin friction. The distance from the leading edge in the streamwise direction is given by x. Both the boundary layer growth rate and the skin friction decrease as the Reynolds number increases. Reynolds analogy can be used to develop a relationship for heat transfer. Reynolds analogy says that the ratio of the shear stress to the heat transfer is a constant near the wall. Thus the Nusselt number can be defined as Nu x = hx C f Re x Pr = 2 k (3.8) Here h is the heat transfer coefficient, k is the thermal conductivity, and Pr is the Prandtl number. This relationship is independent of the equation used to determine skin friction. The power law relationship is not extremely accurate, but is useful for developing some useful turbulent boundary layer relationships. A more accurate relationship for skin friction for adiabatic incompressible flow on a flat plate is given by White and Christoph1 Cf = 0.42 ln (0.056 Re x ) 2 (3.9) This relationship is only applicable in the turbulent region of a boundary layer, and does not apply in the laminar or transitional regions of the boundary layer. The skin friction is affected by a number of parameters, including pressure gradient, surface roughness, compressibility, and surface heat transfer. Adverse pressure gradients cause the skin friction to be reduced as the boundary layer is pushed toward separation. Boundary layer separation occurs when the skin friction becomes negative. High values of skin friction are an indication of a very stable (i.e. difficult to separate) boundary layer. Unfortunately high values of skin friction also equate to high values of viscous drag. The effect of both compressibility and heat transfer on the flat plate turbulent skin friction are shown in Fig. 3.1. Both a hot wall and compressibility tend to reduce the skin friction over the incompressible adiabatic value. 3-3 Figure 3.1 Effect of compressibility and heat transfer on the skin friction on a flat plate. As mentioned above, the power law relationship in Eq. 3.3 is not very accurate. Better approximations of the velocity profile shape are generally written in terms of the parameters u+ and y+ defined as u+ = u uτ (3.10) y+ = ρ w uτ y μw (3.11) and uτ is the friction velocity uτ = τw ρw (3.12) The subscript w denotes the value at the wall and τw is the wall shear stress defined by 3-4 τw = μ ∂u ∂y (3.13) w The boundary layer velocity profile can be divided into four regions. The incompressible velocity profile in each of the subregions regions of the inner region shown in Fig. 3.2 is given by Laminar sublayer 0 < y+ < 5 u+ = y+ (3.14) Buffer layer 5 < y+ < 30 u+ = 5 ln y+ - 3.05 (3.15) Log layer 30 < y+ < 1000 u+ = 1 κ ln y+ + B (3.16) The values of κ, the von Karmen constant, and B are often debated, but are generally accepted to be 0.4 and 5.5 respectively. The log layer is also called the law of the wall. Inner Region Sublayer Buffer Layer Log Layer Outer Region Wake Figure 3.2 Boundary layer regions. The y+ value where the profile transitions from the inner to the outer profile varies with the Reynolds number and the pressure gradient. The outer region is much more sensitive to pressure gradient. Clauser’s2 equilibrium parameter β is often used to characterize the pressure gradient. 3-5 β= δ * ∂pe τ w ∂x (3.17) Coles3 introduced the wake function W given by ⎛π y ⎞ ⎛ y⎞ W ⎜ ⎟ = 2 sin 2 ⎜ ⎟ ⎝2δ⎠ ⎝δ ⎠ (3.18) The velocity profile in the outer region is given by u+ = 1 κ ln y + + B + Π ⎛ y⎞ W⎜ ⎟ κ ⎝δ ⎠ (3.19) where Π is given by Π = 0.8(β + 0.5) 0.75 (3.20) Spalding4 proposed a composite form for the incompressible velocity profile given by ⎡ + (κu + )2 − (κu + )3 ⎤⎥ y + = u + + e −κB e − ΠW ⎢eκu − 1 − κu + − 2 6 ⎥⎦ ⎢⎣ (3.21) White and Christoph1 give a law of the wall velocity profile that includes the effects of compressibility, heat transfer, and pressure gradient. ⎛ 2Γu0+ − Θ ⎞ Γ1 / 2 ⎛ 2Γu + − Θ ⎞ ⎟⎟ + ⎟⎟ = sin −1 ⎜⎜ sin −1 ⎜⎜ Q Q κ ⎝ ⎠ ⎝ ⎠ Where Q = (Θ 2 + 4Γ ) 1/ 2 ⎡ ⎛ φ − 1 φ0 + 1 ⎞⎤ ⎟⎟⎥ ⎢2(φ − φ0 ) + ln⎜⎜ + − φ 1 φ 1 0 ⎝ ⎠⎦ ⎣ and φ = (1 + αy + ) . 1/ 2 (3.22) The values of y 0+ and u0+ are taken as 6 and 10 respectively. The parameters Θ, Γ, and α represent the effects of heat transfer, compressibility, and pressure gradient respectively. Θ= qw μ w Tw k w ρ w uτ (3.23) Γ= ruτ2 2c p Tw (3.24) 3-6 α= μ w ∂pe ρ wτ w uτ ∂x (3.25) Here qw is the wall heat transfer ( k w ∂T ), μw is the wall viscosity, Tw is the wall ∂y w temperature, ρw is the wall density, r is the recovery factor (normally taken as the Prandtl number to the one third power), kw is the wall thermal conductivity, and cp is the specific heat at constant pressure. It is not obvious what effect each of the parameters defined in Eqs. 3.23-25 has on the velocity profile. Fig. 3.3 shows the effect of an adverse pressure gradient and compressibility has on the velocity profile. Increasing α Increasing Γ Figure 3.3 Effect of adverse pressure gradient and compressibility on boundary layer profile shape. The boundary layer thickens and the skin friction decreases as the pressure gradient is increased. Compressibility also causes the skin friction to decrease. Heat transfer effects are shown in Fig. 3.4. 3-7 Increasing Θ (Cold Wall) Decreasing Θ (Hot Wall) Figure 3.4 Effect of wall heat transfer on boundary layer profile shape. A cold wall (wall temperature less than the adiabatic wall temperature) creates a thinner boundary layer and increases the skin friction. A hot wall (wall temperature greater than the adiabatic wall temperature) thickens the boundary layer and decreases the skin friction. The temperature distribution within the inner part of the boundary layer boundary layer is given by the Crocco-Busemann equation ( T = Tw 1 + Θ u + − Γ(u + ) 2 ) (3.26) where Θ is defined in Eq. 3.23 and Γ is defined in Eq. 3.24. For adiabatic wall cases, β =0 and the Crocco-Busemann equation reduces to (γ − 1)r u 2 ⎞ ⎛ T = Tw ⎜ 1 + ⎟ 2 ⎝ ⎠ (3.27) The pressure is assumed to be constant in the normal direction from the wall in the inner part of the boundary layer. Density distributions can be defined based on the temperature distribution and the equation of state. 3.2 Shear Layer Theory A free shear layer is always initiated from a surface of some kind. The boundary layer profile remains for a short period. If no external pressure gradient is present, the shear layer will eventually become self-similar. A self-similar profile is one in which the profile shape remains unchanged as you move downstream if 3-8 the profile is defined in terms of similarity variables. Chapman and Korst5 suggested the similarity variable for two-dimensional shear layers η= σy (3.28) x where x is the downstream distance from the origin of the shear layer, y is the normal distance across the shear layer (y=0 denotes the center of the shear layer, and σ is the shear layer spread rate parameter. Brown and Roshko6 suggested using y* = y − y 0.5 (3.29) δω where the vorticity thickness δω is defined as δω = U1 − U 2 ⎛ ∂u ⎞ ⎜⎜ ⎟⎟ ⎝ ∂y ⎠ max (3.29) Here U1 is the velocity at the high speed edge of the shear layer and U2 is the velocity at the low speed edge of the shear layer. Chapman and Korst5 suggested that the velocity profile was given by u* = u = 0.5[1 + erf (η )] U1 − U 2 (3.30) where erf is the error function. This profile shape is valid for both laminar and turbulent shear layers. Samimy and Elliot7 obtained Laser Doppler Velocimeter (LDV) data on a shear layer. Measurements were made at several downstream locations between the trailing edge of the splitter plate and a station 210 mm downstream of the splitter plate. The flow parameters are given in Table 3.1. Eq. 3.30 is plotted with data in Fig. 3.5. T0, K P01, kPa M1 M2 Mc U1, m/sec U2/U1 ρ 2/ ρ 1 δ1, mm 291.0 314.0 1.80 0.51 0.52 479.5 0.355 0.638 8.0 Table 3.1 Flow parameters for the spatial mixing-layer case 3-9 Figure 3.5 Shear layer velocity profile. The shear layer thickness is given by b= x σ (3.31) The shear layer spread parameter (σ) is affected by compressibility. The accepted value for subsonic flow issuing into quiescent air is σ=11. This value will increase as compressibility effects become greater, and will cause the shear layer to become thinner. Experimental values8 for σ for two-dimensional jets issuing into quiescent air are shown in Fig. 3.6. Figure 3.6 Shear layer spread parameter for jet issuing into quiescent air. 3-10 Subsonic spread parameter experimental results8 as a function of velocity ratio are shown in Fig. 3.7. Figure 3.7 Shear layer spread parameter for subsonic jets. This indicates that the shear layer will become thinner as the ratio of the velocities of the two streams increases. Chapter 3 References: 1. White, F. M. and Christoph, G. H., “A Simple New Analysis of Compressible Turbulent Two-Dimensional Skin Friction Under Arbitrary Conditions, “ AFFDL-TR-70-133, Feb. 1971. 2. Clauser, F. J., “Turbulent Boundary Layers in Adverse Pressure Gradients,” Journal of the Aeronautical Sciences, Vol. 21, 1954, pp. 91108. 3. Coles, D. “The Law of the Wake in the Turbulent Boundary Layer,” Journal of Fluid Mechanics, Vol. 1, Part 2, July 1956, pp. 191-226. 4. Spalding, D. B., “A Single Formula for the Law of the Wall,” Journal of Applied Mechanics, Vol. 28, No. 3, 1961, pp. 444-458. 5. Chapman, A. and Korst, H., “Free Jet Boundary with Consideration of Initial Boundary Layer,” Proceedings of the Second U. S. National Congress of Applied Mechanics, The American Society of Mechanical Engineers, New York, 1954, pp. 723-731. 6. Brown. G. and Roshko, A., “On Density Effects and Large Structure in Turbulent Mixing Layers,” Journal of Fluid Mechanics, Vol. 64, No. 4, 1974, pp. 775-816. 7. Samimy, M. and Elliot, G., “Effects of Compressibility on the Characteristics of Free Shear Layers,” AIAA Journal, Vol. 28, No. 3, 1990, pp. 439-445. 3-11 8. Birch, S. and Eggers, J., “A Critical Review of Experimental Data for Developed Free Turbulent Shear Layers,” in Free Turbulent Shear Flows, NASA-SP-321, 1972, pp. 11-37. 4-1 4.0 Algebraic Turbulence Models 4.1 Inner Eddy Viscosity Model Most of the earliest turbulence models were based on Prandtl’s mixing length hypothesis. Prandtl1 suggested that the eddy viscosity could be represented by μ t = ρL2m ∂u ∂y (4.1) The Prandtl mixing length relates the eddy viscosity to the local mean velocity gradient. The secret to successful applications of the mixing length hypothesis is to find some general method of defining the mixing length. Most algebraic models divide the boundary layer into an inner and outer region as described in Chapter 3. The inner layer includes the viscous sublayer, the buffer layer, and part of the fully turbulent log region. The outer layer includes the remaining part of the log layer and the wake region. The eddy viscosity in the inner layer follows Prandtl's form and is given by (μ t )inner = ρL2m Ω (4.2) where the mixing length Lm is given by ⎡ ⎛ − y+ Lm = κy ⎢1 − exp⎜⎜ + ⎝ A ⎣ ⎞⎤ ⎟⎟⎥ ⎠⎦ (4.3) where y is the distance normal to the wall, κ is the von Karmen constant, and the vorticity Ω is defined as Ω = ∂u ∂v − ∂y ∂x (4.4) 4.2 Cebeci-Smith Model Cebeci-Smith2 suggested that the outer eddy viscosity be expressed as (μt )outer = ραueδ * (4.5) Here α is usually assigned a value of 0.0168 for flows where the Reynolds number based on momentum thickness (Reθ) is greater than 5000, δ* is the displacement thickness, and ue is the velocity at the edge of the boundary layer. The final eddy viscosity is 4-2 μ t = min[(μ t )inner , (μ t )outer ] (4.6) This model is fairly simple, but it requires knowledge of the conditions at the edge of the boundary layer and the boundary layer thickness. These quantities are not always easy to calculate in complicated flows with a Navier-Stokes code since it is often difficult to define where the boundary layer edge actually occurs. 4.3 Baldwin-Lomax Model Baldwin-Lomax3 developed a form of the outer eddy viscosity that did not require knowledge of the conditions at the edge of the boundary layer. This model has become quite popular for CFD applications. The eddy viscosity in the outer layer is defined as (μ t )outer = ρKC cp Fwake Fkleb (4.7) where Fwake contains the mixing length term and Fkleb is the Klebanoff intermittency factor. These terms are defined as { 2 Fwake = min y max Fmax , C wk y maxU diff / Fmax } (4.8) The quantities Fmax and ymax are taken from the maximum of the F function defined as ⎡ ⎛ − y+ F ( y ) = y Ω ⎢1 − exp⎜⎜ + ⎝ A ⎣ ⎞⎤ ⎟⎟⎥ ⎠⎦ (4.9) and Udiff is given by U diff = (U 2 +V 2 ) −( U max 2 +V 2 ) min (4.10) The F function is calculated along a line normal to the wall. The F function for a typical attached boundary layer is shown in Fig. 4.1. 4-3 (Fmax,ymax) y F Figure 4.1 The F function (Eq. 4.9) for an attached boundary layer. The Klebanoff intermittency factor4 is used to reduct the eddy viscosity to zero at the outer edge of the boundary layer. The Klebanoff intermittency factor is given by 6 ⎡ ⎛ C kleb y ⎞ ⎤ ⎟⎟ ⎥ Fkleb ( y ) = ⎢1 + 5.5⎜⎜ ⎢⎣ ⎝ y max ⎠ ⎥⎦ −1 (4.11) The final eddy viscosity is μ t = min[(μ t )inner , (μ t )outer ] (4.12) The constants used in the Baldwin-Lomax model are A + = 26, C cp = 1.6, C kleb = 0.3, C wk = 1.0, κ = 0.4, K = 0.068 (4.13) Algebraic models work well for flows that can be characterized by a single length scale such as attached boundary layers or simple shear layers. The BaldwinLomax model determines the appropriate mixing length from the location of the peak in the F function. For simple flows there will only be one peak in the F function. In more complicated flows with multiple shear layers such as separated boundary layers or wall jets the F function will contain multiple peaks. When multiple peaks are present it becomes difficult to automatically choose the proper peak to use, if a single peak can be used to model the turbulent flow. Degani and Schiff5 recommended a procedure to automatically select the first significant peak in the F function. This modification to the search procedure for the peak has been shown to improve the performance of the Baldwin-Lomax for high angle-of-attack flows with cross-flow separation. 4-4 4.4 Wake and Jet Model The Baldwin-Lomax model has been modified for use in wakes and jets (Ref. 6). First, Udiff in Eq. 4.11 is redefined as U diff = (U 2 +V 2 ) −( U max 2 +V 2 ) y max (4.14) The exponential term in the definition of F(y) (Eq. 4.9) is set to zero yielding F ( y) = y Ω (4.15) The Fwake function (Eq. 4.8) is redefined to be Fwake = y max Fmax ⎛ U diff =⎜ ⎜Ω ⎝ max 2 ⎞ ⎟ Ω max ⎟ ⎠ (4.16) The ymax location in Eq. 4.16 is defined to be the location where Ω max occurs. Finally the Klebanoff intermittency factor Fkleb (Eq. 4.11) is rewritten as Fkleb ⎛ ⎛ ⎜ ⎜ ⎜ ⎜ = ⎜1 + 5.5⎜ C kleb ⎜ ⎜ ⎜ ⎜⎜ ⎜ ⎝ ⎝ ⎞ ⎟ y − y max ⎟ ⎟ ⎛ U diff ⎞ ⎟ ⎜ ⎟⎟ ⎜Ω ⎟⎟ ⎝ max ⎠ ⎠ 6 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ −1 (4.17) This wake formulation will work for individual simple shear layers. If multiple shear layers are present, the user must partition the flow and apply the turbulence model to each shear layer individually. This can easily become impractical for complicated geometries. 4.5 Algebraic Model Shortcomings The flow in a two-dimensional channel with a circular arc bump contraction can be used to demonstrate two shortcomings of the Baldwin-Lomax turbulence models. These shortcomings are the model’s tendency to switch Fmax peaks and the lack of any transport terms in the turbulence model. Fig. 4.2 shows the geometry and the eddy viscosity contours for both the Baldwin-Lomax and Spalart-Allmaras turbulence models. The eddy viscosity predicted by the Baldwin-Lomax is seen to reduce itself almost to zero at the start of the bump. The Spalart-Allmaras one-equation transport turbulence model (discussed in the next section) predicts a much smoother distribution of eddy viscosity on the bump. The reason for the anomaly in the Baldwin-Lomax model can be seen in 4-5 Fig. 4.3. The location of the Fmax peak moves very close to the wall at the start of the bump (x=0) and rises back to its original level at the bump trailing edge (x=1). The relative magnitude of the Fmax peak is also shown if Fig. 4.3, and is seen to increase in the region of the bump. The eddy viscosity is proportional to the product of the distance of the peak off the wall and the magnitude of the peak, and is seen to decrease in the region above the bump because the Fmax peak has moved very near the wall. The sudden switching of the Fmax peak causes the turbulence level to be greatly reduced in the favorable pressure gradient region near the bump leading edge. If the Baldwin-Lomax model included transport terms, then the higher upstream values of eddy viscosity would be transported downstream at the beginning of the bump, and the Baldwin-Lomax predicted eddy viscosity distribution would not demonstrate the discontinuous behavior it shows here. Figure 4.2 Eddy viscosity contours for a circular arc bump in a two-dimensional channel for two different turbulence models. 4-6 Figure 4.3 Distance above the wall of the Fmax peak and the Fmax value for a twodimensional channel with a circular arc bump. 4.6 Grid Sensitivity for a Flat Plate with Adiabatic Wall The initial wall spacing of the computational grid and the grid-stretching ratio (the ratio of the change in grid spacing normal to the wall) can affect the accuracy of the Baldwin-Lomax model. Figure 4.4 shows the sensitivity of the skin friction to initial wall spacing for a flat plate. The grid-stretching ratio was 1.2 for all these cases. The plots include the theoretical skin friction curves of White and of Spalding. 4-7 Figure 4.4 Flat plate skin friction predictions for the Baldwin-Lomax turbulence model for varying initial wall grid point spacings. The calculated boundary layer is seen to become fully turbulent around a length Reynolds number (Rex) of 1x106. The results for y+=0.2 and y+=1 are virtually identical indicating a grid independent solution. The y+=5 solution shows some small divergence from the y+=1 solution at the lower length Reynolds numbers while the y+=10 solution shows large differences from the other solutions. Predicted velocity profiles for the flat plate boundary layer for various initial wall grid point spacings are shown in Fig. 4.5. The velocity profile shows little effect of the initial spacing for all but the y+=10 profile. All of the profiles but the y+=10 profile are in good agreement with the theoretical profile from Spalding. Note that the theoretical profile does not include the law-of-the wake (see Chapter 3), and hence the predicted profiles diverge from the theoretical profile in the wake region of the boundary layer. The predicted eddy viscosity for various initial wall spacings is shown in Fig. 4.6. 4-8 Figure 4.5 Flat plate boundary layer profiles predicted by the Baldwin-Lomax turbulence model for varying initial wall grid point spacings. Figure 4.6 Eddy viscosity distribution predicted by the Baldwin-Lomax turbulence model for varying grid initial wall spacings. Here again it is seen that the y+=0.2 and the y+=1.0 results are almost identical. The y+=5 and y+=10 results show the solutions are no longer grid independent at larger wall spacings. 4-9 The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig. 4.7. All of these solutions used an initial wall spacing of y+=1. Figure 4.7 The effect of grid stretching ratio on the skin friction for a flat plate boundary layer using the Baldwin-Lomax turbulence model. There seems to be very little effect of grid stretching for this case indicating that the initial wall spacing is the more critical parameter for skin friction predictions for flat plates with the Baldwin-Lomax turbulence model. This is also the case for the velocity profile, as seen in Fig. 4.8. The eddy viscosity does change as the stretching ratio increases as shown in Fig. 4.9. It is interesting to note that a wide range of eddy viscosity distributions have little effect on skin friction and the velocity profile for a flat plate boundary layer. This insensitivity to the absolute eddy viscosity level is probably a major reason why the eddy viscosity concept has worked so well in practice. 4-10 Figure 4.8 The effect of grid stretching ratio on the velocity profile for a flat plate boundary layer using the Baldwin-Lomax turbulence model. Figure 4.9 The effect of grid stretching ratio on the eddy viscosity distribution for a flat plate boundary layer using the Baldwin-Lomax turbulence model. 4.7 Grid Sensitivity for Axisymmetric Bump A second example of the grid sensitivity of the Baldwin-Lomax turbulence model that includes a pressure gradient is the NASA Ames transonic axisymmetric 4-11 bump experiment (Ref. 7). The geometry is shown in Fig. 4.10. The model consisted of a sharp-lipped hollow cylinder with a 15.2 cm outer surface diameter. The bump was a circular arc 20.3 cm long and 1.9 cm high that begins 60.3 cm downstream of the cylinder leading edge. The upstream intersection of the bump and cylinder was faired with a circular arc. The test was run at a Mach number of 0.875 and a chord Reynolds number of 2.67x106. M=0.875 Shock Recirculation D=15.2 cm c=20.3 cm h=1.9 cm Figure 4.10 Geometry for the transonic axisymmetric bump. The effect of initial grid spacing on the pressure coefficient distribution along the bump is shown in Fig. 4.11. The stretching ratio was 1.2 for these cases. The pressure coefficient seems to be relatively insensitive to the initial grid spacing, with the y+=10 and y+=20 curves being slightly displaced from the other curves. The velocity distribution at the aft junction of the bump and the cylinder (x/c=1) is shown in Fig. 4.12. The y+=20 solution predict a larger velocity in the reverse flow region than the other solutions. Grid stretching effects on the pressure coefficient distribution along the bump is shown in Fig. 4.13. The initial grid spacing was 1.2 for these cases. The pressure distribution coefficient changes slightly as the grid-stretching ratio is increased to 1.5. The solution in the separated region differs greatly for a grid-stretching ration of 2.0. The effect on the velocity distribution at x/c=1 is shown in Fig. 4.14. As with increasing initial grid spacing, increasing the grid spacing increases the size and the magnitude of the separated flow region. 4-12 Figure 4.11 The effect of initial grid wall spacing on the pressure coefficient for the axisymmetric bump using the Baldwin-Lomax turbulence model. Figure 4.12 The effect of initial grid wall spacing on the velocity distribution at x/c=1 for the axisymmetric bump using the Baldwin-Lomax turbulence model. 4-13 Figure 4.13 The effect of grid stretching on the pressure coefficient for the axisymmetric bump using the Baldwin-Lomax turbulence model. Figure 4.14 The effect of grid stretching on the velocity profile at x/c=1 for the axisymmetric bump using the Baldwin-Lomax turbulence model. The results for the Ames axisymmetric bump indicate that the grid-stretching ratio is a critical parameter when pressure gradients are present in the flow. The stretching ratio should probably be kept between 1.2 and 1.3 to assure that the grid can capture the pressure gradient effects. This is true in both structured and unstructured grids. 4-14 4.8 Grid Sensitivity for a Flat Plate with Heat Transfer Calculating heat transfer accurately can be more difficult than predicting skin friction. This can be seen in the subsonic flat plate example when the wall temperature is specified to be 1.5 times the free-stream temperature. The sensitivity of the skin friction and heat transfer result with varying initial grid wall spacing is shown in Fig. 4.15 and Fig. 4.16. The grid stretching ratio was fixed at 1.2 for these results. Both the skin friction and heat transfer seem to be relatively insensitive to the wall spacing for wall spacings less than y+=5 when no pressure gradient is present. Figure 4.15 The effect of wall spacing on the skin friction on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. 4-15 Figure 4.16 The effect of wall spacing on the heat transfer (Stanton number) on a flat plate using the Baldwin-Lomax turbulence model. The velocity and temperature profiles for a length Reynolds number (Rex) of 1.0x107 are shown in Fig. 4.17 and 4.18 respectively. Both the velocity and temperature profiles are relatively insensitive to wall spacing for this model when no pressure gradient is present. Fig. 4.17 The effect of wall spacing on the velocity profile on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. 4-16 Fig. 4.18 The effect of wall spacing on the temperature profile on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. Grid stretching effects on skin friction and heat transfer predictions are shown in Fig. 4.19 and 4.20. The effect of grid stretching on the velocity and temperature profiles is shown in Fig. 3.22 and 3.23 respectively. The initial wall spacing was held at y+=0.1 for these calculations. The results reach a grid independent result for a stretching ratio of less than 1.3. Figure 4.19 The effect of grid stretching on the skin friction on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. 4-17 Figure 4.20 The effect of grid stretching on heat transfer (Stanton number) on a flat plate using the Baldwin-Lomax turbulence model. Figure 4.21 The effect of grid stretching on the velocity profile on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. 4-18 Figure 4.22 The effect of grid stretching on the temperature profile on a flat plate with heat transfer using the Baldwin-Lomax turbulence model. 4.9 Grid Sensitivity for a Nozzle with Heat Transfer Flow through a supersonic nozzle with a constant temperature wall can serve as a test case for evaluating the performance of the turbulence model in the presence of strong pressure gradients. Back, Massier, and Gier8 measured the wall pressure distribution and heat transfer for a converging-diverging nozzle with a throat diameter of 0.0458 meters and an exit diameter of 0.1227 meters. Highpressure air was heated by the internal combustion of methanol and flowed along a cooled constant area duct with a length of 0.4572 meters and a diameter of 0.355 meters before entering the nozzle. The nozzle geometry and boundary conditions are shown in Fig. 4.23. The gas could be treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The nozzle exit Mach number was approximately 2.5. The molecular viscosity and thermal conductivity were assumed to vary according to Sutherland’s law. Pt=150psia Tt=1418oR Twall=0.5Tt Figure 4.23 Nozzle geometry. Supersonic Outlet 4-19 The initial wall spacing was varied with a grid stretching ratio of 1.2 in the boundary layer. The grid spacing was held constant in the core of the nozzle. Comparisons of the predicted and measured pressure along the nozzle are shown in Fig. 4.24 for varying initial wall spacings. The throat is located at x=0.091 meters. The pressure distribution is seen to be insensitive to the initial wall spacing. The heat transfer at the wall is shown in Fig. 4.25. The results are in poor agreement with the data and are quite sensitive to the initial wall spacing, especially for wall spacing greater than y+=1. Figure 4.24 The effect of wall spacing on the pressure distribution for a supersonic nozzle with heat transfer using the Baldwin-Lomax turbulence model. Figure 4.25 The effect of wall spacing on the wall heat transfer distribution for a supersonic nozzle using the Baldwin-Lomax turbulence model. 4-20 R/Rthroat The poor performance of the Baldwin-Lomax model for this case is due to the difficulty in choosing the proper Fmax peak as was described in Section 4.5. The F function at the nozzle throat is shown in Fig. 4.26. Notice that multiple peaks are present in the function. It is extremely difficult to predict which, if any, of these three peaks is the proper value for this case. The eddy viscosity distribution at the nozzle throat is shown in Fig. 4.27. The Spallart-Allmaras eddy viscosity distribution is included in Fig. 4.26 for comparison. The Baldwin-Lomax eddy viscosity is much lower than the Spalart-Allmaras distribution and is seen to cut off prematurely. This low value of eddy viscosity results in a low prediction of the wall heat transfer. Figure 4.26 The F function at the nozzle throat. 4-21 Fig. 4.27 The eddy viscosity distribution at the nozzle throat. 4.10 Summary Algebraic turbulence models were popular in the 1970’s and 1980’s because of their simplicity and robustness. As Navier-Stokes CFD applications became more complex in the 1990’s these models began to lose popularity because of accuracy limitations for flows that contain multiple shear layers or boundary layer separation. The difficulties encountered when multiple peaks occur in the F function have been demonstrated. These models have also lost favor in unstructured grid applications since they require a velocity profile over multiple grid points aligned with the flow for successful application. The eddy viscosity predicted by an algebraic turbulence model is only a function of the local velocity profile used to generate the F function. Thus the eddy viscosity relates directly to the local instantaneous vorticity field of that profile and cannot model the transport of turbulence by the flow. This makes these models a poor choice for unsteady flows and for flows where turbulent transport is important. 4-22 Baldwin-Lomax Application Hints 1. The Baldwin-Lomax model requires that the F function be well defined. This normally requires that at least three points be located within the sublayer (y+<10). The first point off the wall should be located about y+<5 for pressure distributions, y+<2 to obtain reasonable skin friction values, and y+<0.5 for heat transfer. The grid stretching normal to the wall should not exceed 1.3. Improved heat transfer results can be obtained by using a constant spacing for the first three cells off the wall. 2. In order to reduce the probability of finding a second peak well off the wall, it is usually good to limit the number of points over which the F function is calculated. 3. Care should be taken not to divide viscous regions such as boundary layers when dividing the computational domain for blocked or chimera applications since the entire velocity profile is required to properly define the Fmax and Udiff quantities. Chapter 4 References: 1. Prandtl, L., “Bericht uber Untersuchugen zur Ausgebildeten Turbulenz,” Z. Angew Math., Meth., Vol. 5, pp. 136-139. 2. Cebeci, T, and Smith, A. M. O., Analysis of Turbulent Boundary Layers, Academic Press, New York, 1974. 3. Baldwin, B. S. and Lomax, H., “Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows,” AIAA-78-0257, Jan. 1978. 4. Klabenoff, P. S., “Characteristics of Turbulence in a Boundary Layer with Zero Pressure Gradient,” NACA Rep. 1247, 1955. 5. Degani, D. and Schiff, L. B., “Computation of Supersonic Viscous Flows Around Pointed Bodies at Large Incidence,” AIAA-83-0034, Jan. 1983. 6. Buning, P. G., Jesperson, D. C., Pulliam, T. H., Klopfer, G. H., Chan, W. C., Slotnik, J. P., Krist, S. E., and Renze, K. J., “OVERFLOW User’s Manual,” Version 1.8l, July 1999. 7. Johnson, D. A., “Predictions of Transonic Separated Flow with an EddyViscosity/Reynolds Shear Stress Closure Model,” AIAA-85-1683, July 1985. 8. Back, L. H., Massier, P. F., and Gier, H. L., “Convective Heat Transfer in a Convergent-Divergent Nozzle,” International Journal of Heat Mass Transfer, Vol. 7, 1964, pp. 549-568. 5-1 5.0 One-Equation Turbulence Models 5.1 Theory The general form for all transport turbulence models is similar to the form of a species equation in reacting flows ∂Z ∂Z 1 ∂ ⎡ ∂Z ⎤ +Ui = ⎢(ν + ν t ) ⎥ + P( Z ) − D( Z ) ∂t ∂xi σ ∂x j ⎣ ∂xi ⎦ (5.1) where Z is a general turbulence variable. The terms on the left hand side of Eq. 5.1 represent the convective transport of Z. The first term on the right hand side of Eq. 5.1 is the diffusion of Z assuming a gradient-diffusion model. The functions P(Z) and D(Z) on the right hand side of Eq. 5.1 represent the production and destruction of Z. The production function for turbulence models are usually functions of the eddy viscosity and the fluid strain or vorticity. A functional relationship between Z and the eddy viscosity is required to complete the model. Several one-equation transport models for turbulence have appeared over the years. The most popular one-equation models have been the Baldwin-Barth1 model and the Spalart-Allmaras model2. The Spalart-Allmaras model will be discussed in some detail here. 5.2 Spalart-Allmaras Model The Spalart-Allmaras turbulence model was derived using empirical relationships, dimensional analysis, and Galilean invariance. The goal was to produce a turbulent transport model that was fast, numerically stable, and reasonably accurate for both shear layers and boundary layers. The model uses a turbulence variable ν~ that has the dimensions of viscosity. The model can be written as [ ] ∂ν~ ∂ν~ 1 2 +Ui = ∇ ⋅ ((ν + ν~ )∇ν~ ) + C b 2 (∇ν~ ) + P(ν~ ) − D(ν~ ) ∂t ∂xi σ (5.2) where the production term is given by ν~ ⎛ ⎞ P(ν~ ) = C b1ν~⎜ Ω + 2 2 f v 2 ⎟ κ d ⎝ ⎠ (5.3) and the dissipation is given by ⎛ ν~ ⎞ D(ν~ ) = C w1 f w ⎜ ⎟ ⎝d⎠ 2 (5.4) 5-2 Here Ω is the magnitude of the vorticity, d is the distance to the nearest wall, and fv2 and fw are given by f v2 = 1 − χ 1 + χf v1 (5.5) ⎡ 1+ C6 ⎤ f w = g ⎢ 6 w36 ⎥ ⎣ g + C w3 ⎦ 1 6 (5.6) and the remaining functions are given by f v1 = χ= χ3 (5.7) χ 3 + C v31 ν~ ν (5.8) ( g = r + C w2 r 6 − r r= (Ωκ ) (5.9) ν~ 2 d 2 + ν~f v 2 ) (5.10) The remaining constants are given in Table 5.1. Cb1 0.1355 Cb2 0.622 σ κ 2/3 0.41 Cw1 Cb1 / κ + (1 + Cb 2 ) / σ 2 Cw2 0.3 Cν1 7.1 Table 5.1 Coefficients for the Spalart-Allmaras model. The eddy viscosity is defined as ν t = ν~f v1 (5.11) The original model also includes a boundary layer transition model. The trip function can be used for simple geometries, but is difficult to apply for complex geometries and is usually omitted. The turbulent transition location can also be fixed by setting the production term to zero in the region of the flow upstream of the desired transition location. The wall boundary condition is ν~ = 0 . At free stream boundaries, ν~ is set to a small number, generally ν/10. The turbulence variable ν~ is usually limited to be greater than the free stream value in the field so that it will not go negative. The model will self-initialize from this small value of ν~ in the freestream. The 5-3 numerical implementation for implicit algorithms is discussed in detail in Ref. 2. Generally the convective terms are treated with upwind differencing. The diffusion term is treated with central differencing. The full source term jacobians ∂P(ν~ ) ∂D(ν~ ) ) should be used in implicit formulations. When the r function ( ~ , ∂ν ∂ν~ becomes large, the fw function becomes a constant, so r can be clipped at a value of 10 to avoid numerical overflow problems. The Spalart-Allmaras model suffers from a deficiency common to transport type models in that it fails to predict the reduction in shear layer growth rate with increasing jet Mach number. The model will predict the incompressible shear layer growth rate regardless of the jet Mach number. Paciorri and Sabetta3 suggest a compressibility correction for this model. 5.3 Rotation/Streamline Curvature Corrections Conventional linear eddy viscosity models have difficulty predicting flows with large system rotation or streamline curvature. The Spalart-Allmaras model falls into this category. The vorticity reaches a local maximum in the core of a vortex causing the eddy viscosity to increase rapidly. This leads to excessive dissipation of the vortex core. To remedy this problem the production term (Eq. 5.3) be modified by multiplying by a correction factor4. ν~ ⎛ ⎞ P(ν~ ) = Cb1ν~⎜ Ω + 2 2 f v 2 ⎟ Fr1 (r *) κ d ⎝ ⎠ (5.12) where Fr1(r*) is a correction for rotational flow given Fr1 (r*) = [1 + Cvor min (0, r * −1)] (5.13) and r* = S +ε (5.14) Ω +ε where |S| is the magnitude of the strain. The strain is defined as ⎛ ∂u ∂u j S ij = 0.5⎜ i + ⎜ ∂x j ∂xi ⎝ ⎞ ⎟ ⎟ ⎠ (5.15) The magnitude of the strain is given by S 2 = 2S ij S ij (5.16) 5-4 The ε in Eq. (8) is a small number used as a threshold value so that the turbulence model returns to its baseline form in areas of low vorticity and strain such as near free stream boundaries. The ε value can be defined as ε = 0.5 U ref (5.17) Lref Where Uref and Lref are the reference velocity and length scale for the calculation respectively. This correction form implies that curvature effects are a function of the ratio of the local strain magnitude and vorticity magnitude. The strain goes to zero in the core of a vortex while the vorticity reaches a local maximum. Thus Eq. (5.13) converts the production term of the SA model into a dissipation term in the vortex core region if Cvor is greater than one. Ref. 4 recommended Cvor=2. The modified form of the production term has almost no effect in boundary layers and shear layers since the magnitude of the strain and the magnitude of the vorticity are almost equal for these flows and the correction term goes to zero. This correction is generally referred to as the approximate Spalart-Allmaras rotational correction (ASARC). Shur, et al.5 introduced a correction to the Spalart-Allmaras turbulence model for rotating and curved flows. The correction is applied to the production term in a manner similar to the ASARC correction. The correction has the form [ ] 2r * 1 − c r 3 tan −1 (cr 2 ~ r ) − c r1 Fr1 (r*, ~ r ) ) = (1 + c r1 ) 1+ r * (5.18) where 2Ω ik S jk ~ r = D4 ⎛ DS ij c ⎞ ⎜ ⎟ ⎜ Dt + ε imn S jn + ε jmn S in Ω m ⎟ ⎝ ⎠ ( ) (5.19) The remaining term D is define as ( D 2 = 0.5 S 2 + Ω 2 ) (5.20) The recommended values for the constants are cr1=1.0, cr2=12.0 and cr3=1.0. The Lagrangian derivative of the strain tensor in Eq. (5.19) is included to obtain the proper value for the curvature correction. This term makes the correction expensive and cumbersome to implement. Some implementations of this correction ignore the contribution of the time derivative of the Lagrangian derivative and the reference frame rotation to simplify the coding. This is referred as the SARC correction. 5-5 An example of the standard SA, SARC, and ASARC models applied to flow in a u-duct6 is shown in Fig. 5.1. The distance along the wall is denoted by s. The standard SA model under-predicts the pressure coefficient and the skin friction at the exit of the u-bend. Both the SARC and ASARC models improve the comparison with data. Figure 5.1. Pressure coefficient and skin friction on the outer wall of a u-duct.. The standard SA model also overdamps the flow in the core of a vortex. This tends to cause vortices to be damped out prematurely in vortex dominated flows such as wing or rotor tip vortices. The production term for the SA model (Eq. 5.3) is based on the vorticity magnitude. The vorticity magnitude reaches a local maximum in a vortex core, hence the model tends to produce much more eddy viscosity than is required to properly simulate the flow in the vortex core. The SARC and ASARC correction terms are negligible in simple shear layers and boundary layers and do not affect the model since the magnitude of the strain is approximately equal to the magnitude of the vorticity. The correction term transforms the production term into a dissipation term in the core of a vortex where the magnitude of the strain is zero. This additional dissipation tends to reduce the eddy viscosity in the core of the vortex and significantly improves the simulation of vortex dominated flows. The wing tip vortex on a NACA 0012 wing4 tip at an angle-of-attack of 10o is shown in Fig. 5.2 for both the original and modified SA production terms. The tangential velocity predicted by the modified SA models is in much better agreement with the data. Fig. 5.2b also includes results for the hybrid RANS/LES Detached Eddy Simulation (DES) models that are discussed in Chapter 9. 5-6 Tangential Velocity Figure 5.2a Axial velocity for a NACA 0012 wing tip vortex. Vortex Radius Fig. 5.2b Vortex tangential velocity. 5.4 Grid Sensitivity for a Flat Plate with Adiabatic Walls The initial wall spacing of the computational grid and the grid-stretching ratio can affect the accuracy of the Spalart-Allmaras model. Figure 5.3 shows the sensitivity of the skin friction to initial wall spacing for a flat plate. The gridstretching ratio was 1.2 for all these cases. The plots include the theoretical skin friction curves of White and of Spalding. 5-7 The boundary layer is seen to become fully turbulent around a length Reynolds number (Rex) of 1x106. The results for y+=0.2 and y+=1 are virtually identical indicating a grid independent solution. The y+=5 solution show some small divergence from the y+=1 solution at the lower length Reynolds numbers while the y+=10 solution shows large differences from the other solutions. Predicted velocity profiles for the flat plate boundary layer for various initial wall grid point spacings are shown in Fig. 5.4. The velocity profile shows little effect of the initial spacing for all but the y+=10 profile. All of the profiles but the y+=10 profile are in good agreement with the theoretical profile from Spalding. Note that the theoretical profile does not include the law-of-the wake (see Chapter 3), and hence the predicted profiles diverge from the theoretical profile in the wake region of the boundary layer. The predicted eddy viscosity for various initial wall spacings is shown in Fig. 5.5. Figure 5.3 Flat plate skin friction predictions for the Spalart-Allmaras turbulence model for varying initial wall grid point spacings. 5-8 Figure 5.4 Flat plate boundary layer profiles predicted by the Spalart-Allmaras turbulence model for varying initial wall grid point spacings. Figure 5.5 Eddy viscosity distribution predicted by the Spalart-Allmaras turbulence model for varying grid initial wall spacings. Here again it is seen that the y+=0.2 and the y+=1.0 results are almost identical. The y+=5 and y+=10 results show the solutions are no longer grid independent at larger wall spacings. The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig. 5.6. All of these solutions used an initial wall spacing of y+=1. 5-9 There seems to be very little effect of grid stretching for these cases indicating that the initial wall spacing is the more critical parameter for skin friction predictions for flat plates with the Spalart-Allmaras turbulence model. This is also the case for the velocity profile, as seen in Fig. 5.7. The eddy viscosity does change as the stretching ratio increases as shown in Fig. 5.8. It is interesting to note that a wide range of eddy viscosity distributions have little effect on skin friction and the velocity profile for a flat plate boundary layer. Figure 5.6 The effect of grid stretching ratio on the skin friction for a flat plate boundary layer using the Spalart-Allmaras turbulence model. Figure 5.7 The effect of grid stretching ratio on the velocity profile for a flat plate boundary layer using the Spalart-Allmaras turbulence model. 5-10 Figure 5.8 The effect of grid stretching ratio on the eddy viscosity distribution for a flat plate boundary layer using the Spalart-Allmaras turbulence model. 5.5 Grid Sensitivity for and Axisymmetric Bump A second example of the grid sensitivity of the Spalart-Allmaras turbulence model that includes a pressure gradient is the NASA Ames transonic axisymmetric bump experiment described in the previous chapter. The effect of initial grid spacing on the pressure coefficient distribution along the bump is shown in Fig. 5.9. The stretching ratio was 1.2 for these cases. The pressure coefficient seems to be relatively insensitive to the initial grid spacing for y+<10. The y+=20 result shows a significant difference from the other results in regions where a pressure gradient is present. The velocity distribution at the aft junction of the bump and the cylinder (x/c=1) is shown in Fig. 5.10. The y+=20 solution predict a larger velocity in the reverse flow region than the other solutions. Grid stretching effects on the pressure coefficient distribution along the bump is shown in Fig. 5.11. The initial grid spacing was 1.2 for these cases. The pressure distribution coefficient changes slightly as the grid-stretching ratio is increased to 1.5. The solution in the separated region differs greatly for a gridstretching ration of 2.0. The effect on the velocity distribution at x/c=1 is shown in Fig. 5.12. As with increasing initial grid spacing, increasing the grid spacing increases the size and the magnitude of the separated flow region. 5-11 Figure 5.9 The effect of initial grid wall spacing on the pressure coefficient for the axisymmetric bump using the Spalart-Allmaras turbulence model. Figure 5.10 The effect of initial grid wall spacing on the velocity distribution at x/c=1 for the axisymmetric bump using the Spalart-Allmaras turbulence model. 5-12 Figure 5.11 The effect of grid stretching on the pressure coefficient for the axisymmetric bump using the Spalart-Allmaras turbulence model. Figure 5.12 The effect of grid stretching on the velocity profile at x/c=1 for the axisymmetric bump using the Spalart-Allmaras turbulence model. As was seen with the Baldwin-Lomax model, grid stretching is a critical parameter for this turbulence model when adverse pressure gradients are present. Care should be taken to keep the grid-stretching ratio between 1.2 and 1.3. 5-13 5.6 Grid Sensitivity for a Flat Plate with Heat Transfer Calculating heat transfer accurately can be more difficult than predicting skin friction. This can be seen in the subsonic flat plate example when the wall temperature is specified to be 1.5 times the free-stream temperature. The sensitivity of the skin friction and heat transfer result with varying initial grid wall spacing is shown in Fig. 5.13 and Fig. 5.14. The grid stretching ratio was fixed at 1.2 for these results. Both the skin friction and heat transfer seem to be relatively insensitive to the wall spacing for initial wall spacings less than y+=1. Figure 5.13 The effect of wall spacing on the skin friction on a flat plate with heat transfer using the Spalart-Allmaras turbulence model. The effect of wall spacing on the velocity and temperature profiles for a length Reynolds number (Rex) of 1.0x107 is shown in Fig. 5.15 and 5.16 respectively. The profiles are relatively insensitive to wall spacing for initial wall spacings less than y+=1.0. 5-14 Figure 5.14 The effect of wall spacing on the heat transfer (Stanton number) on a flat plate using the Spalart-Allmaras turbulence model. Fig. 5.15 The effect of wall spacing on the velocity profile for a flat plate with heat transfer using the Spalart-Allmaras turbulence model. 5-15 Fig. 5.16 The effect of wall spacing on the temperature profile for a flat plate with heat transfer using the Spalart-Allmaras turbulence model. Grid stretching effects on skin friction and heat transfer predictions are shown in Fig. 5.17 and 5.18. The effect on velocity and temperature profiles for a length Reynolds number (Rex) of 1.0x107 are shown in Fig. 5.19 and 5.20. The initial wall spacing was held at y+=0.1 for these calculations. The results reach a grid independent result for stretching ratios of less than 1.3. Figure 5.17 The effect of grid stretching on the skin friction on a flat plate with heat transfer using the Spalart-Allmaras turbulence model. 5-16 Figure 5.18 The effect of grid stretching on heat transfer (Stanton number) on a flat plate using the Spalart-Allmaras turbulence model. Figure 5.19 The effect of grid stretching on the velocity profile on a flat plate with heat transfer using the Spalart-Allmaras turbulence model. 5-17 Figure 5.20 The effect of grid stretching on the temperature profile on a flat plate with heat transfer using the Spalart-Allmaras turbulence model. 5.7 Grid Sensitivity for a Nozzle with Heat Transfer Flow through a supersonic nozzle with a constant temperature wall can serve as a test case for evaluating the performance of the turbulence model in the presence of strong pressure gradients. Details of the geometry and boundary conditions for the converging-diverging supersonic nozzle are given in Chapter 4. High-pressure air was heated by the internal combustion of methanol and flowed along a cooled constant area duct before entering the nozzle. The gas could be treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The nozzle exit Mach number was 2.5. The molecular viscosity and thermal conductivity were assumed to vary according to Sutherland’s law. The grid initial wall spacing was varied and the grid stretching ratio was held at 1.2 in the boundary layer. A uniform grid was used in the nozzle core. Predicted wall pressure distribution results for varying initial wall spacings are shown in Fig. 5.21. The pressure distribution is seen to be insensitive to the initial wall spacing. Predicted wall heat transfer is shown in Fig. 5.22. The results are somewhat sensitive to the initial wall spacing for values of y+ less than one. The predicted results diverge rapidly from the data for wall spacing greater than y+=1. 5-18 Figure 5.21 The effect of wall spacing on the pressure distribution for a supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence model. Figure 5.22 The effect of wall spacing on the wall heat transfer for a supersonic nozzle using the Spalart-Allmaras turbulence model. The grid stretching ratio was varied in the boundary layer while the initial grid spacing held at y+=0.5. A uniform grid was used in the nozzle core. Predicted wall pressure distribution results for varying initial wall spacings are shown in Fig. 5.23. The pressure distribution is seen to be insensitive to the grid stretching 5-19 ratio. Predicted wall heat transfer is shown in Fig. 5.24. The results are insensitive to the grid stretching ratios less than 1.3. The predicted results diverge slightly from the data for grid stretching ratios greater than 1.3. Figure 5.23 The effect of grid stretching ratio on the pressure distribution for a supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence model. Figure 5.24 The effect of grid stretching ratio on the wall heat transfer for a supersonic nozzle using the Spalart-Allmaras turbulence model. 5-20 5.8 Summary The Spalart-Allmaras turbulence model is a very stable and generally reasonably accurate model for a wide range of turbulent flows. The model is relatively easy to implement in both structured and unstructured Navier-Stokes codes. One drawback is the model requires the calculation of the distance to the nearest wall for all field points. This can be an expensive computation, especially for unstructured grid codes. The model has been used extensively for threedimensional geometries and is well documented. The model has been used with some success for some unsteady flows. The example given here are for subsonic and supersonic flows. An excellent validation source for hypersonic flows can be found in Ref. 7. Spalart-Allmaras Application Hints 1. The first point off the wall should be located about y+=1 to obtain reasonable skin friction values and about y+=0.5 for heat transfer. The grid stretching normal to the wall should not exceed 1.3. 2. The eddy viscosity should be limited so that it will not run away in some complex applications. Generally a limit of νt/ν=200,000 is acceptable. 3. Care should be taken not to divide viscous regions such as boundary layers when dividing the computational domain for blocked or overset applications since the model requires the distance from the nearest wall. 4. This model tends to smear out three-dimensional vortical flows. Rotation and curvature corrections (SARC and ASARC) can significantly improve the results. 5. The model can overdamp some unsteady flows. 6. The model contains no corrections for compressibility and will overpredict the growth rate of high speed shear layers. Chapter 5 References: 1. Baldwin B. S. and Barth, T. J.,”A One-Equation Turbulence Transport Model for High Reynolds Number Wall-Bounded Flows,” AIAA-91-0610, Jan. 1991. 2. Spalart, P. R. and Allmaras, S. R., “A One-Equation Turbulence Model for Aerodynamic Flows,” AIAA-92-0439, Jan. 1992. 3. Paciorri, R. and Sabetta, F., “Compressibility Correction for the SpalartAllmaras Model in Free-Shear Flows,” Journal of Spacecraft and Rockets, Vol. 40, No. 3, May-June 2003. 4. Dacles-Mariani, J., Zilliac, G., Chow, J., and Bradshaw, P., “Numerical/Experimental Study of a Wingtip Vortex in the Near Field,” AIAA Journal, Vol. 33, No. 9, September 1995, pp. 1561-1568. 5-21 5. Shur, M. L., Strelets, M. K., Travin, A. K., and Spalart, P. R., “Turbulence Modeling in Rotating and Curved Channels: Assessing the Spalart-Shur Correction,” AIAA Journal, Vol. 38, No. 5, May 2000, pp.784-792. 6. Nichols, R., “Algorithm and Turbulence Model Requirements for Simulating Vortical Flows,” AIAA-2008-0337, Jan. 2008. 7. Roy, J. and Blottner, F., “Methodology for Turbulence Model Validation: Application to Hypersonic Flows,” Journal of Spacecraft and Rockets, Vol. 40, No. 3, May-June 2003. 6-1 6.0 Two-Equation Turbulence Models Two-equation turbulence models have been used for almost fifty years. There are a tremendous number of these models in the literature. Most of these models solve a transport equation for turbulent kinetic energy (k) and a second transport equation that allows a turbulent length scale to be defined. The most common forms of the second transport equation solve for turbulent dissipation (ε) or turbulent specific dissipation (ω). Some of these two-equation models are valid down to the wall (low Reynolds number models) and some are only valid outside the inner region of the boundary layer (high Reynolds number models). All of these models have their own strengths and weaknesses. Several of the currently popular models are described in this chapter. 6.1 Theory A transport equation for turbulent kinetic energy can be derived from NavierStokes equations and has the form ∂k ∂k ∂ +Ui =− ∂t ∂xi ∂x j ⎡ ⎛ u' u' ' ⎞⎤ ∂U i ∂u i' ∂u i' ⎢u i' ⎜ j j + p ⎟⎥ − u i' u 'j −ν ρ ⎟⎠⎥ ∂x j ∂x j ∂x j ⎢ ⎜⎝ 2 ⎣ ⎦ (6.1) The left hand side of Eq. 6.1 contains the temporal change and convective terms for the turbulent kinetic energy. The first term on the right hand side of Eq. 6.1 is the diffusive transport of the turbulent stresses (Dk). The second term on the right hand side of Eq. 6.1 is the production of turbulent kinetic energy (Pk). The last term on the right hand side of Eq. 6.1 represents the dissipation of turbulent kinetic energy (ε). All of the terms on the right hand side must be modeled. Applying the Boussinesq approximation (Eq. 2.22) the production term becomes ⎛ ∂U ∂U j Pk = ν t ⎜ i + ⎜ ∂x ∂xi ⎝ j ⎞ ∂U i 2 ∂U i ∂U i 2 ∂U i 2 ∂U i ⎟ − k = ν t S ij S ij − k − k = ν t S ij ⎟ ∂x 3 ∂xi ∂x j 3 ∂xi ⎠ j 3 ∂xi (6.2) The last term in Eq. 6.2 is zero for incompressible flow since the divergence of ∂U the velocity ( i ) is zero. This term is also neglected in many compressible flow ∂xi applications since it can cause numerical difficulties near strong shocks and expansions. Sij is the strain rate defined as 1 ⎛⎜ ∂ui ∂u j S ij = + 2 ⎜⎝ ∂x j ∂xi ⎞ ⎟ ⎟ ⎠ 2 Kato and Launder1 suggested redefining the production term as (6.3) 6-2 Pk = υ t S ij Ω ij (6.4) Where Ωij is the vorticity defined as 1 ⎛⎜ ∂ui ∂u j ⎞⎟ − 2 ⎜⎝ ∂x j ∂xi ⎟⎠ Ω ij = 2 (6.5) The purpose of this modification is to reduce the tendency of many twoequation models to overpredict the turbulent production in regions with large normal strain, i.e. regions with strong acceleration or deceleration. This correction is often used in stagnation regions such as at the leading edge of a wing or blade. The diffusion term can be modeled using the gradient-diffusion approximation as Dk = ⎡⎛ νt ⎢⎜⎜ν + σk ⎣⎝ ∂ ∂xi ⎞ ∂k ⎤ ⎟⎟ ⎥ ⎠ ∂xi ⎦ (6.6) where σk can be thought of as a turbulent Schmidt number. Relationships for the eddy viscosity (νt) and the turbulent dissipation (ε) are required for closure. The square root of the turbulent kinetic energy can be used to represent a velocity scale for the large-scale turbulent motion. Using this in an eddy viscosity relationship yields νt ≈ kL (6.7) where L is a turbulent length scale. A relationship for the turbulent length scale of the large scale eddies is given in Eq. 1.13. Using this definition of the turbulent length scale in Eq. 6.7 produces ν t = Cμ k2 ε (6.8) where Cμ is a constant. The system of equations would be closed with a relationship for the turbulent dissipation (ε). The two features that distinguish different two equation models are the treatment of the wall and the form of the turbulent dissipation equation. It is possible to derive a transport equation for turbulent dissipation from the Navier-Stokes equations, but the equation contains fluctuating correlation terms that cannot be 6-3 easily modeled. Turbulent dissipation cannot be measured directly experimentally, but must be inferred from measurements of other correlations. Because of the difficulty in deriving a turbulent dissipation equation based on physics, turbulence modelers have traditionally resorted to using a transport equation with empirically determined production and dissipation terms. The two most popular equations used to close the turbulence model are based on the dissipation ε2 and a turbulence variable ω defined by Wilcox3 as ω= ε (6.9) Cμ k The eddy viscosity for the Wilcox formulation is νt = k (6.10) ω 6.2 Traditional k-ε Models The standard form of the ε transport equation is ∂ ∂ε ∂ε = + Ui ∂xi ∂xi ∂t ⎡⎛ ν t ⎞ ∂ε ⎤ ε ε2 ⎟ +R ⎥ + Cε 1 f ε 1 Pk − Cε 2 f ε 2 ⎢⎜⎜ν + σ ε ⎟⎠ ∂xi ⎦ k k ⎣⎝ (6.11) The production and dissipation terms of the dissipation equation are formed from the production and dissipation terms of the turbulent kinetic energy equation scaled by ε /k and multiplied by empirically determined constants and wall damping functions (Cε1fε1 and Cε2 fε2). An additional damping term must be included for the eddy viscosity in the k-ε equations for applications near walls so that k and ε will have the proper behavior in the near wall region. The constants Cε1, Cε2, and Cμ are empirically determined from comparisons with experimental data. In grid-generated turbulence the diffusion term and the production term in Eq. 6.1 and Eq. 6.11 are negligible and Cε2 can be determined from the measured rate of decay of the turbulent kinetic energy k behind the grid. The value of Cε2 was found to lie between 1.8 and 2.0. For shear-layers in local equilibrium (production=dissipation), Eqs. 6.1 and 6.8 can be combined to produce ⎛ u ′v ′ ⎞ ⎟ C μ = ⎜⎜ ⎟ ⎝ k ⎠ 2 (6.12) Measurements in equilibrium flows suggest that Cμ~0.09. In the near-wall region of the logarithmic region of a boundary layer turbulent production and dissipation 6-4 are approximately equal and the convection of ε is negligible. assumptions Eq. 6.11 reduces to Cε 1 = Cε 2 − Under these κ2 (6.13) σ ε Cμ Note that all of these empirical constants are defined for near equilibrium conditions, and thus are not truly valid when production and dissipation are not of the same size. This is often a source of numerical difficulties when solving these models and leads to the use of ad hoc relationships to control the behavior of the turbulence model for cases when production is much greater than dissipation. The general form of the eddy viscosity relationship for low Reynolds number models is ν t = Cμ f μ k2 (6.14) ε Many k-ε turbulence transport models exist. The coefficients of some of the more popular are shown in Table 6.1. The wall damping functions are shown in Table 6.2. Model Standard JonesLaunder2 Chien4 Speziale5 Cε1 1.44 1.45 Cε2 1.92 1.92 σk σε 1.0 1.0 1.35 1.47 1.80 1.83 1.0 1.0 κ 1.3 1.3 Cμ 0.09 0.09 0.43 0.43 1.3 1.56 0.09 0.09 0.42 0.41 Table 6.1 Coefficients for low Reynolds number k-ε models. Model Standard JonesLaunder2 Chien4 fε1 1 2ν ⎛ ∂ k ⎜ 1+ ε ⎜⎝ ∂y 2ν k 1+ 2 yε fε2 1 ⎞ ⎟ ⎟ ⎠ 2 2νν ⎛ ∂ 2 u ⎞ 1 − .3 exp(− R ) + 2 t ⎜⎜ 2 ⎟⎟ ε ⎝ ∂y ⎠ 2 t ⎡ ⎛ R 2 ⎞⎤ Cε 2 ⎢1 − .222 exp⎜⎜ − t ⎟⎟⎥ ⎝ 36 ⎠⎦ ⎣ + ⎛ y+ ⎞ 2ν k ⎜⎜ − ⎟⎟ exp y 2ε ⎝ 2 ⎠ fμ 1 2 ⎤ ⎡ 3.4 exp ⎢ − 2⎥ ⎣ (1 + 0.02 Rt ) ⎦ ⎛ y+ ⎞ ⎟⎟ 1 − exp⎜⎜ − ⎝ 87 ⎠ 6-5 Speziale5 1 ⎡ ⎛ Rt2 ⎞⎤ ⎜⎜ − ⎟⎟⎥ 1 − . 222 exp ⎢ 36 ⎝ ⎠⎦ ⎣ ⎡ ⎛ y + ⎞⎤ ⎟⎟⎥ * ⎢1 − exp⎜⎜ − ⎝ 4.9 ⎠⎦ ⎣ + ⎛ ⎞ ⎜1 + 3.45 ⎟ tanh ⎛⎜ y ⎞⎟ ⎜ 70 ⎟ ⎜ Rt ⎟⎠ ⎝ ⎠ ⎝ 2 Table 6.2 Wall damping functions for low Reynolds number k-ε models. The last column of Table 6.1 also contains the von Karman constant (κ) calculated from Eq. 6.10. All of the models shown here produce κ values slightly above the usually accepted value of 0.4-0.41 with the exception of the Speziale model. The wall boundary condition for the turbulent kinetic energy is k wall = 0 (6.15) Since the turbulent kinetic energy is zero at the wall, then the eddy viscosity is also zero at the wall. A wall boundary condition for the dissipation (ε) can be found from the turbulent kinetic energy equation ε wall = ∂ ∂xi ⎛ ∂k ⎞ ⎜⎜ν ⎟⎟ ⎝ ∂xi ⎠ (6.16) This boundary condition can produce numerical stability problems, so Speziale4 suggested the boundary condition be replaced with the form ε wall ⎛∂ k = 2ν ⎜⎜ ⎝ ∂xi ⎞ ⎟ ⎟ ⎠ 2 (6.17) The turbulent kinetic energy is usually assigned to the free stream turbulence level at free stream boundaries if it is known. If the free stream turbulence level is not known, then the turbulent kinetic energy is set to a small number at the free stream boundary. The eddy viscosity is also set to a small number at free stream boundaries. νt < 0.1 ν∞ ∞ (6.18) 6-6 Once the free stream turbulent kinetic energy and eddy viscosity are known, the free stream value of turbulent dissipation may be calculated from the definition of eddy viscosity. 6.3 k-ω Model Using Wilcox’s formulation, a transport equation for ω can be written as ∂ ∂ω ∂ω = + Ui ∂x i ∂x i ∂t ⎡⎛ ν t ⎞ ∂ω ⎤ ⎡ κ2 ⎜ ⎟ ⎢ + − ν + β ⎥ ⎢⎜ σ ω ⎟⎠ ∂x i ⎦ ⎢ σ ϖ Cμ ⎣⎝ ⎣ ⎤ω ⎥ Pk − βω 2 ⎥⎦ k (6.19) The ω transport equation is essentially the ε transport equation (Eq. 6.11) with the addition of a cross diffusion term of the form ⎛ ν Cross Diffusion = ⎜⎜ν + t σε ⎝ ⎞ ∂ε ∂k ⎟⎟ ⎠ ∂xi ∂xi (6.20) The k-ω models are derived for wall bounded flows and require no additional wall damping terms when used in boundary layer flows. The coefficients for Eq. 6.19 are given in Table 6.3. β σω 0.075 2 Cμ 0.09 κ 0.41 Table 6.3 Coefficients for the k-ω model. The wall boundary condition for the turbulent kinetic energy is also given by Eq. 6.15. The definition of ω (Eq. 6.9) suggests that ω should go to infinity at the wall since the turbulent kinetic energy (k) is going to zero. Wilcox, based on asymptotic arguments, suggested that ω be given the value ω wall = 60ν C μ (Δy 2 ) (6.20) at the wall. Here Δy is the normal distance to the first point off the wall. The following free stream boundary conditions are recommended ν t∞ U = 0.001, ω∞ = 10 ∞ , k ∞ = ν t∞ω∞ ν∞ L (6.21) The k-ω model has been shown to be sensitive to the free stream value of ω chosen. 6-7 6.4 SST Model In practice, the k-ε models are generally more accurate in shear type flows and are well behaved in the far field. The k-ω models are more accurate and much more numerically stable in the near wall region. Recognizing that each model has its strength and weakness and that the forms of the equations are similar, Menter6 suggested a blended model that is a k-ω model near the wall and transitions to a k-ε model away from the wall. The model has the form ∂k ∂k ∂ ⎡ ∂k ⎤ +Ui = ⎢(ν + σ kν t ) ⎥ + Pk − C μ ωk ∂t ∂xi ∂xi ⎣ ∂xi ⎦ (6.22) ∂ω ∂ω ∂ ⎡ ∂ω ⎤ ω 2 +Ui = ⎢(ν + σ ων t ) ⎥ + γ Pk − βω + ∂t ∂xi ∂xi ⎣ ∂xi ⎦ k (1 − F1 ) 2σ ω ∂k ∂ω ω ∂xi ∂xi (6.23) The coefficients are blended forms of the two baseline models. The blending function F1 is defined as ( F1 = tanh arg14 ) (6.24) where ⎡ ⎛ k 500ν arg 1 = min ⎢ max⎜ , ⎜ C ωy y 2ω μ ⎢⎣ ⎝ ⎞ 4 ρσ ω 2 k ⎤ ⎟, ⎟ CD y 2 ⎥ ⎥⎦ kω ⎠ (6.25) Here y is the normal distance to the wall and CDkω is the positive portion of the cross-diffusion term ⎛ ⎞ 1 ∂k ∂ω CDkω = max⎜⎜ 2 ρσ ω 2 ,10 − 20 ⎟⎟ ω ∂x i ∂x i ⎝ ⎠ (6.26) The coefficients σk, σω, γ, and β in Eq. 6.23 are computed with the general form φ = F1φ1 + (1 − F1 )φ 2 (6.27) where the φ1 corresponds to coefficients from the k-ω model and φ2 corresponds to coefficients of the k-ε model. The eddy viscosity is calculated from νt = a1 k max (a1ω , ΩF2 ) (6.28) 6-8 where F2 is given by ( F2 = tanh arg 22 ) (6.29) and ⎛ 2 k 500ν arg 2 = max⎜ , ⎜ C ω y 2ω ⎝ μ ⎞ ⎟ ⎟ ⎠ (6.30) Here Ω is the magnitude of the vorticity vector. The first term in the denominator of Eq. 6.28 comes from the traditional definition of eddy viscosity for a k-ω model (Eq. 6.10). The second term in the denominator represents an attempt to improve the performance of the model for adverse pressure gradients and comes from a one equation turbulence modeled developed by Bradshaw7 based on a relationship for shear stress in a boundary layer similar to Eq. 6.12 − u ' v' = a1k (6.31) The second term in the denominator of Eq. 6.28 reduces the eddy viscosity in the logarithmic and wake region of the boundary layer when adverse pressure gradients are present. The use of the shear stress relationship of Eq. 6.31 in the definition of the eddy viscosity (Eq. 6.28) has given rise to calling this model the shear stress transport (SST) model. The coefficients used with the SST model are given Tables 6.4 and 6.5. Set 1 k-ω Set 2 k-ε β σk σω 0.075 0.85 0.5 β1 − Cμ C μ (Cε 2 − 1) 1 γ σ ω1κ 2 Cμ Cε 1 − 1 0.857 Table 6.4 SST model functions. Cε1 1.44 Cε2 1.92 Cμ 0.09 κ 0.41 a1 0.31 Table 6.5 SST model coefficients. The boundary conditions for SST model are the same as the k-ω model. The SST model is relatively insensitive to the free stream value of ω. 6-9 6.5 RNG Model One of the major criticisms of the k-ε model is that it is not derived from the Navier-Stokes equations in any systematic fashion. Yakhot and Orszag8 applied Renormalization Group (RNG) methods to derive the k-ε equations. In this approach, an expansion is made about an equilibrium state with known Gaussian statistics using the correspondence principle that the effects of the mean strains can be represented by a random force. Bands of high wave numbers (small scales) are systematically removed and space is rescaled in a manner analogous to that used in phase transitions. The successive removal of larger scales ultimately leads to a k-ε model of the form νt ∂k ∂k ∂ ⎡⎛ + Ui = ⎢⎜⎜ν + ∂t ∂xi ∂xi ⎣⎝ σk ⎞ ∂k ⎤ ⎟⎟ ⎥ + Pk − ε ⎠ ∂xi ⎦ (6.32) ν t ⎞ ∂ε ⎤ ∂ε ∂ε ∂ ⎡⎛ ε ε2 ⎟⎟ + Ui = +R ⎢⎜⎜ν + ⎥ + Cε 1 Pk − Cε 2 ∂t ∂xi ∂xi ⎣⎝ σ ε ⎠ ∂xi ⎦ k k (6.33) where R is defined as R=− ζ (1 − ζ / 4.38) Pk ε (1 + 0.012ζ 3 ) k ζ = (2 S ij Sij )1 / 2 (6.34) k (6.35) ε The parameter ζ is the ratio of the turbulent to mean strain time scale. The remaining constants are given in Table 6.6. Cε1 1.42 Cε2 1.68 Cμ 0.085 σk σe 0.719 0.719 Table 6.6 RNG model coefficients. The RNG model is similar to the traditional k-ε model with the exception of the additional term R and the lower value of Cε2. These two differences decrease both the rate of production of turbulent kinetic energy and the rate of dissipation of ε, leading to lower values of the eddy viscosity. In free shear flows, ζ (and thus R) is zero. The model presented is the high Reynolds number form. The additional terms required for the low Reynolds form can be found in Ref. 9. 6-10 6.6 Numerical Implementation The flux terms on the right hand side of Eqs. 6.1, 6.11, and 6.19 are usually treated with upwind differences. The diffusion terms are generally treated with central differences. The turbulent kinetic energy production term Pk is limited to be less than 20ε in the turbulent kinetic energy source term to keep the model from “running away” in complex flows. The maximum value of eddy viscosity is also capped (normally < 2x105ν) to limit non-physical behavior of the model. The turbulent quantities (k, ε, and ω) are limited to be positive in the field. The exact source term Jacobian term can be used for the ε or ω equations for implicit computations. Use of the exact source term for the Jacobian of the turbulent kinetic energy equation usually results in numerical instability. It is recommended that the turbulent kinetic energy source term Jacobian be replaced with TKE source term Jacobian = − 2 Pk + ε k (6.36) The two equations need not be solved fully coupled, and can be lagged from the mean flow equations. For time accurate implicit solutions the turbulence equations should be solved within a Newton iteration loop with the mean flow equations to assure that the solution is locally converged at each time step. 6.7 Compressibility Correction for Shear Layers The standard k-ε and k-ω turbulence models will over predict the mixing for highspeed compressible shear layers. The models tend to predict the incompressible growth rate for the shear layer rather than decreasing the growth rate as compressible effects increase. A common method of correcting the k-ε model for this effect is given in Ref. 10. The compressible turbulent kinetic energy equation is written as μt ∂ρk ∂ρU i k ∂ ⎡⎛ + = ⎢⎜⎜ μ + ∂t ∂xi ∂xi ⎣⎝ σk ⎞ ∂k ⎤ ⎟⎟ ⎥ + ρPk − ρ (ε + ε c ) + p" d " ⎠ ∂xi ⎦ (6.37) where ε c = α 1 M t2 ε (6.38) Mt is the turbulent Mach number defined as Mt = 2k γRT (6.39) 6-11 The last term in Eq. 6.37 represents the pressure dilatation. modeled as p" d " = −α 2 ρPk M t2 + α 3 ρεM t2 This term is (6.40) Based on direct numerical simulation, the recommended values for the constants in Eqs. 6.38 and 6.40 are α 1 = 1.0,α 2 = 0.4, α 3 = 0.2 (6.41) The compressibility correction must also be included in the SST model since the SST model reduces to the standard k-ε model away from the wall. A method for incorporating the compressibility correction on the SST model is given in Ref. 11. Further corrections for high temperature jet flows are described in Ref. 12. The need for a compressibility correction can be seen in predictions of the supersonic axisymmetric jet shear layer of Eggers13. The jet exit Mach number was 2.22 and the jet exit static pressure was matched to the quiescent outer air. Fig. 6.1 shows the predicted and measured axial velocity on the jet centerline. Fig. 6.2 contains velocity comparisons at various downstream locations in the jet. Figure 6.1. Centerline axial velocity for the Ref. 12 jet. 6-12 Figure 6.2a. Velocity profile at X/R=25 for the Ref. 12 jet. Figure 6.2b. Velocity profile at X/R=50 for the Ref. 12 jet. 6-13 Figure 6.2c. Velocity profile at X/R=100 for the Ref. 12 jet. Without the compressibility correction all of the turbulence models predict a faster growth rate than is seen in the experimental data. The SST model with the compressibility correction is seen to do a good job for accounting for the compressibility effect. 6.8 Initializing Turbulence Values for a Given Profile Shape There is really no unique way to construct turbulence variables for a given velocity profile shape for an inflow boundary condition for transport type turbulence models. Here are two simple approaches. Often a combination of these approaches is required in real world applications. One approach is to specify a turbulent intensity and a turbulent length scale for a given profile. The turbulent intensity is specified as a percentage (a) of the local time-averaged fluid velocity u ′ = aU (6.42) The turbulent intensity is normally less than twenty per cent of the mean flow velocity. A turbulent length scale (L) must then be specified. Some typical choices are shown in Table 6.7. Flow L Boundary Layer 0.41δ Plane Jet Round Jet 0.53δ 0.44δ Plane Wake 0.95δ Pipe 0.25D 6-14 Table 6.7 Typical values for the turbulent length scale L. In Table 6.7, δ is defined as the distance between points where the velocity differs from the free stream velocity by one percent of the maximum velocity difference across the layer. For symmetrical flows and boundary layers, δ is the distance from the symmetry plane or wall to the one percent point at the outer edge. D is the internal pipe diameter. The turbulent kinetic energy can be found from k= 3 2 3 2 u ′ = (aU ) 2 2 (6.43) The turbulent length scale is related to the eddy viscosity by the KolmogorovPrandtl expression μt = Cμ ρ k L (6.44) The traditional definition of eddy viscosity for two-equation turbulence models is given by μt = Cμ ρ k2 ε (6.45) Thus the turbulent dissipation is defined as ε= k 3/ 2 L (6.46) Although this method is relatively simple to implement, it is not very accurate for boundary layers and looses meaning for complex velocity profile shapes. The second method takes advantage of the turbulence equilibrium assumption (turbulent production = turbulent dissipation). This requires an estimate of the eddy viscosity along the profile. The eddy viscosity profile can be obtained from a simple algebraic turbulence model such as Baldwin-Lomax or a one-equation turbulence model such as Spalart-Allmaras. The turbulent dissipation is then given by ε = Pk (6.47) Finally, the turbulent kinetic energy is obtained from the definition of eddy viscosity in Eq. (6.45) 6-15 k= μtε Cμ ρ (6.48) This method works well for general applications and can be useful for starting k-ε models. These models often have difficulties generating eddy viscosity when a solution is started from a uniform free stream condition. 6.9 Rotation and curvature correction A rotational and curvature correction for the SST model can also be formulated in a manner similar to the SARC and ASARC corrections used in the SA oneequation model. The correction has a much smaller effect because the SST uses strain rather than vorticity in the definition of the production term. The production term used in both the k and ω equations (Eq. 6.22 and 6.23) is redefined as Pk = μ t S ij S ij Fr1 (r*, ~ r) (6.49) The corrections defined in Eq. (5.13) and Eq. (5.18) are used with Eq. (6.49) to provide a general three-dimensional rotation and curvature correction for the SST model in this study. The Eq. (5.13) and Eq. (5.18) corrections are limited 0.05. 6-24 Figure 6.16 The effect of wall spacing on the velocity profile on a flat plate with heat transfer using the SST turbulence model. Figure 6.17 The effect of wall spacing on the temperature profile on a flat plate with heat transfer using the SST turbulence model. Grid stretching effects on skin friction and heat transfer predictions are shown in Fig. 6.18 and 6.19. The velocity and temperature profiles for a length Reynolds number (Rex) of 1x107 are shown in Fig. 6.20 and 6.21. The initial wall spacing was held at y+=0.1 for these calculations. The results reach a grid independent result for stretching ratios less than 1.3. 6-25 Figure 6.18 The effect of grid stretching on the skin friction on a flat plate with heat transfer using the SST turbulence model. Figure 6.19 The effect of grid stretching on heat transfer (Stanton number) on a flat plate using the SST turbulence model. 6-26 Figure 6.20 The effect of grid stretching on the velocity profile on a flat plate with heat transfer using the SST turbulence model. Figure 6.21 The effect of grid stretching on the temperature profile on a flat plate with heat transfer using the SST turbulence model. 6.13 Grid Sensitivity for a Nozzle with Heat Transfer Flow through a supersonic nozzle with a constant temperature wall can serve as a test case for evaluating the performance of the turbulence model in the presence of strong pressure gradients. Details of the geometry and boundary 6-27 conditions for the converging-diverging supersonic nozzle are given in Chapter 4. High-pressure air was heated by the internal combustion of methanol and flowed along a cooled constant area duct before entering the nozzle. The gas could be treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The nozzle exit Mach number was 2.5. The molecular viscosity and thermal conductivity were assumed to vary according to Sutherland’s law. The grid initial wall spacing was varied and the grid stretching ratio was held at 1.2 in the boundary layer. A uniform grid was used in the nozzle core. Predicted wall pressure distribution results for varying initial wall spacings are shown in Fig. 6.22. The pressure distribution is seen to be insensitive to the initial wall spacing. Predicted wall heat transfer is shown in Fig. 6.23. The results are somewhat sensitive to the initial wall spacing for values of y+ less than one. The predicted results diverge rapidly from the data for wall spacing greater than y+=1. Figure 6.22 The effect of wall spacing on the pressure distribution for a supersonic nozzle with heat transfer using the SST turbulence model. 6-28 Figure 6.23 The effect of wall spacing on the wall heat transfer for a supersonic nozzle using the SST turbulence model. The grid stretching ratio was varied in the boundary layer while the initial grid spacing held at y+=0.5. A uniform grid was used in the nozzle core. Predicted wall pressure distribution results for varying initial wall spacings are shown in Fig. 6.24. The pressure distribution is seen to be insensitive to the grid stretching ratio. Predicted wall heat transfer is shown in Fig. 6.25. The results are insensitive to the grid stretching ratios less than 1.3. The predicted results diverge slightly from the data for grid stretching ratios greater than 1.3. 6-29 Figure 6.24 The effect of grid stretching ratio on the pressure distribution for a supersonic nozzle with heat transfer using the SST turbulence model. Figure 6.25 The effect of grid stretching ratio on the wall heat transfer for a supersonic nozzle using the SST turbulence model. 6.14 Summary Two-equation models have been the industry standard for many years. These models have been used for a wide variety of internal and external flows. Many 6-30 variants are available that have been tuned for specific classes of flows. These models are more difficult to implement than are the one-equation models and can often be less numerically stable. For these reasons many code developers have avoided two-equation models. The examples given here are mainly for subsonic and supersonic applications. An excellent summary of hypersonic validation of these models can be found in Ref. 15. Two-Equation Model Application Hints 1. 2. 3. 4. 5. 6. The first point off the wall should be located about y+=1 to obtain reasonable skin friction values and about y+=0.5 for heat transfer. The grid stretching normal to the wall should not exceed 1.3. A constant spacing should be used for the first three cells off the wall for heat transfer calculations. The eddy viscosity should be limited so that it will not run away in some complex applications. Generally a limit of νt/ν=200,000 is acceptable. Care should be taken not to divide viscous regions such as boundary layers when dividing the computational domain for blocked or overset applications since the model requires the distance from the nearest wall. k-ε models are often difficult to start. It is usually easier to come off the ground with another turbulence model such as Baldwin-Lomax or Spalart-Allmaras and switch to the k-ε model after νt/ν>1.0 in the computational domain. This can be done by first calculating the turbulent kinetic energy production Pk with Eq. 6.2. Assuming the turbulence is in local equilibrium, then the dissipation can be found from ε =Pk. The turbulent kinetic energy k can then be found from the definition of eddy viscosity in Eq. 6.8. This model can overdamp some unsteady flows. Compressibility corrections should be included for high-speed shear layer flows. Chapter 6 References: 1. Kato, M. and Launder, B. E., "The Modeling of Turbulent Flow Around Stationary and Vibrating Square Cylinders", Proc. 9th Symposium on Turbulent Shear Flows, Kyoto, August 1993, pp. 10.4.1-10.4.6. 2. Jones, W. P., and Launder, B. E., “The Calculation of Low-Reynolds Number Phenomena with a Two-Equation Model of Turbulence,” Int. J. of Heat and Mass Transfer, Vol. 16, 1973, pp. 1119-1130. 6-31 3. Wilcox, D. C., “Reassessment of the Scale-Determining Equation for advanced Turbulence Models, AIAA Journal, Vol. 26, No. 11, 1988, pp. 1299-1310. 4. Chien, J. Y., “Predictions of Channel Boundary-Layer Flows with a LowReynolds-Number Turbulence Model,” AIAA Journal, Vol. 20, Jan. 1982, pp. 33-38. 5. Speziale, C., Abid, R., and Anderson, E., “A Critical Evaluation of TwoEquation Models for Near Wall Turbulence,” ICASE Report No. 90-46, June 1990. 6. Menter, F. R. and Rumsey, C. L., “Assessment of Two-Equation Turbulence Models for Transonic Flows,” AIAA-94-2343, June 1994. 7. Bradshaw, P., Ferriss, D. H., and Atwell, N. P., “Calculation of Boundary Layer Development Using the Turbulent Energy Equation,” Journal of Fluid Mechanics, Vol. 28, 1967, pp. 593-616. 8. Yakhot, V., Orszag, S., Thangam, S., Gatski, T., and Speziale, C. G., “Development of Turbulence Models for Shear Flows by a Double Expansion Technique,” Physics of Fluids A, Vol. 4, No. 7, 1992, pp. 15101520. 9. Barber, T., Choi, D., Nedungadi, A., Orszag, S., Konstantinov, A., and Staroselsky, I., “Advanced Turbulence Models for High-Mach Number Wall Bounded Flows,” AIAA-98-0536, Jan. 1998. 10. Speziale, C., Abid, R., and Anderson, E., “A Critical Evaluation of TwoEquation Models for Near Wall Turbulence,” ICASE Report No. 90-46, June 1990. 11. Suzen, Y. B. and Hoffmann, K. A., “Investigation of Supersonic Jet Exhaust Flow by One- and Two-Equation Turbulence Models,” AIAA-980322, Jan. 1998. 12. Abdol-Hamid, K. Pao, S. Massey, S., and Elmiligui, A., “Temperature Corrected Turbulence Model for High Temperature Jet Flow,” AIAA-20034070, Jun. 2003. 13. Eggers, J. “Velocity Profiles and Eddy Viscosity Distributions Downstream of a Mach 2.22 Nozzle Exhausting to Quiescent Air,” NASA-TN D-3601, September 1966. 14. Nichols, R., “Algorithm and Turbulence Model Requirements for Simulating Vortical Flows,” AIAA-2008-0337, Jan. 2008. 15. Roy, J. and Blottner, F., “Methodology for Turbulence Model Validation: Application to Hypersonic Flows,” Journal of Spacecraft and Rockets, Vol. 40, No. 3, May-June 2003. 7-1 7.0 Reynolds and Algebraic Stress Models 7.1 Reynolds Stress Models One means of obtaining closure for the Reynolds or Favre averaged NavierStokes equations is to investigate higher moments of the equations themselves. Transport equations for each of the Reynolds stresses can be derived from the momentum equations. These equations are called the second moment equations. The second moment equations can be obtained by multiplying the momentum equations by Uj and then Reynolds averaging the resulting equations. The incompressible form of the second moment equation is ∂ ui' u 'j ∂t ∂ ∂x k + Uk ∂ ui' u 'j ∂x k ⎛ ⎞ ∂U j ⎟ ⎜ ' ' ∂U i ' ' = −⎜ u j u k + ui u k ⎟− ∂x k ∂x k ⎟ ⎜ ⎝ ⎠ ( ) ⎡ ' ' ' ⎤ p' ' ∂ ' ' ( ν u u u u u uiδ jk + u 'jδ ik )⎥ + − + ⎢ i j k i j ρ ∂xk ⎢⎣ ⎥⎦ ' p ' ⎛⎜ ∂ui' ∂u j + ρ ⎜⎝ ∂x j ∂xi (7.1) ' ⎞ ∂u 'j ⎟ − 2ν ∂ui ⎟ ∂x k ∂x k ⎠ The first term on the right hand side of Eq. 7.1 is the turbulent production term. The second term on the right hand side is the diffusion of the Reynolds stresses. The third term on the right side is called the pressure strain. The pressure strain term redistributes energy among the turbulent stress components. The last term on the right hand side is the turbulent dissipation. The balance of these source terms is critical to the proper simulation of a turbulent flow. The second moment equations produce a host of new unknowns and higher order correlations that must be modeled in order to close the system. These higher order correlations are difficult, and in many cases impossible, to determine from experimental measurements. Direct Numerical Simulation (DNS) has been used to determine these quantities in many cases. Unfortunately DNS is restricted to low Reynolds numbers and simple geometries at the present time, so the data base for determining these higher order correlations, whether experimentally or numerically generated, is limited. Eq. 7.1 may be rewritten as ∂ ui' u 'j ∂t + Uk ∂ ui' u 'j ∂x k = Pij + ∂ ∂ ( ( Dijkv ) + Dijkt ) + Π ij − ε ij ∂x k ∂xk The production term is defined using the exact form (7.2) 7-2 ⎛ ⎞ ∂U j ⎟ ⎜ ' ' ∂U i ' ' Pij = −⎜ u j uk + ui u k ⎟ ∂x k ∂xk ⎟ ⎜ ⎝ ⎠ (7.3) Molecular diffusion also follows the exact definition Dijkv = ν ( ) ∂ ' ' ui u j ∂x k (7.4) There are several models suggested for the turbulent diffusion. Harlow1 suggested t = Cs Dijk ' ' k ⎛⎜ ' ' ∂ ui u j uk ul ∂x l ε⎜ ⎝ ⎞ ⎟ ⎟ ⎠ Daly and (7.5) Hanjalic and Launder2 used the form Dijkt = C s ' ' ∂ ui' u 'j ⎞⎟ ∂u ' u ' k ⎛⎜ ' ' ∂ u j uk ui u l + u 'j ul' k i + uk' ul' ∂xl ε⎜ ∂xl ∂xl ⎟ ⎠ ⎝ (7.6) Mellor and Herring3 suggests ' ' ' ' k 2 ⎛⎜ ∂u j uk ∂ uk' ui' ∂ ui u j D = Cs + + ε ⎜ ∂xi ∂x j ∂x k ⎝ t ijk ⎞ ⎟ ⎟ ⎠ (7.7) Many models for the pressure-strain term have been developed since the pioneering work of Launder, Reece, and Rodi4. A general form is given by 1 ⎞ ⎛ Π ij = −C1bij + C1′ε ⎜ bik bkj − bmn bmnδ ij ⎟ + C2 kS ij 3 ⎠ ⎝ 2 ⎞ ⎛ +C3k ⎜ bik S jk + b jk S ik − bmn S mnδ ij ⎟ 3 ⎠ ⎝ +C 4 k (bikW jk + b jkWik ) (7.8) where the stress anisotropy tensor is given by bij = ui′u ′j 1 − δ ij 2k 3 (7.9) 7-3 The strain rate tensor is defined as S ij = 1 ⎛⎜ ∂ui ∂u j + 2 ⎜⎝ ∂x j ∂xi ⎞ ⎟ ⎟ ⎠ (7.10) and the rotation tensor is given by Wij = 1 ⎛⎜ ∂ui ∂u j − 2 ⎜⎝ ∂x j ∂xi ⎞ ⎟ ⎟ ⎠ (7.11) If the turbulence is assumed to be locally isotropic (i.e. the same amount of turbulent energy is dissipated by each of the turbulent energy components u i' 2 ) then the dissipation term in Eq. 7.1 can be reduced to 2 3 ε ij = εδ ij (7.12) where ε is the total rate of energy dissipation. A transport equation for the dissipation must also be solved to close the system of equations. The dissipation equation can be written as ∂ε ∂ε ε ε2 ∂ + Uk = Cε 1 Pε − Cε 2 + ∂t k k ∂x k ∂x k ⎛ ∂ε ⎞ ∂ ⎟⎟ + ⎜⎜ν (Dεtj ) x x ∂ ∂ k ⎠ i ⎝ (7.13) The first term on the right hand side of Eq. 7.13 is the production of dissipation. The second term represents the dissipation of the dissipation. The last two terms are the molecular and turbulent diffusion of dissipation respectively. The production of dissipation is given by ⎛ ∂U i ⎞⎟ Pε = −⎜ ui' u 'j ⎜ ∂x j ⎟⎠ ⎝ (7.14) The dissipation equation turbulent diffusion term is modeled by ⎛k ∂ε ⎞⎟ Dεtj = Cε ⎜ uk' ul' ⎜ε ∂x j ⎟⎠ ⎝ (7.15) Typical values for the constants in the above equations for the Launder-ReeceRodi4 (LRR), Gibson-Launder5 (GL) and Speziale-Sarker-Gatski6 (SSG) RSM models are 7-4 Model LRR GL SSG C1 1.5 3.6 3.4+1.8P/ε C1′ 0.0 0.0 4.2 C2 0.4 0.8 0.8-1.3(bijbij)1/2 C3 0.6 1.2 1.25 C4 0.0 1.2 0.4 Cs 0.25 0.11 0.11 Cε 1 1.44 1.44 1.44 Cε 2 1.92 1.92 1.83 Cε 0.15 0.11 0.11 Table 7.1 Coefficients for three Reynolds stress models. A system of seven transport equations (six Reynolds stress components and one turbulent dissipation equation) must be solved at each time step to provide the Reynolds stresses for the Navier-Stokes equations. The above correlations must be further modified when applied in the presence of a wall, and hence many of these models utilize wall functions (see Chapter 10) to avoid the complications and the additional grid points introduced by the presence of the wall. The Reynolds stresses produced by these models must meet certain constraints for the models to be applicable to a large range of problems without significant tuning. Two of the more important constraints are tensor invariance and realizability. Tensor invariance requires the replaced terms to have the same tensor form as the original terms. This will insure that the modeled terms transform properly in different coordinate systems. Schumann7 introduced the realizability constraint. This constraint requires the equation for the turbulent stresses to have the property that all the component energies of the turbulent kinetic energy (the diagonal terms of the Reynolds stress tensor) remain nonnegative and all off-diagonal components of the Reynolds stress tensor satisfy Schwartz’s inequality. This can be written as ui' ui' ≥ 0 (7.16) 2 ui' ui' u 'j u 'j − ui' u 'j ≥ 0 (7.17) ( ) 2 u1' u1' ⎛⎜ u2' u2' u3' u3' − u2' u3' ⎞⎟ − u1' u2' u1' u2' u3' u3' − u2' u3' u1' u3' + ⎠ ⎝ ( ) (7.18) u u u u u u −u u u u ≥0 ' ' 1 3 ' ' 1 2 ' ' 2 3 ' ' 2 2 ' ' 1 3 The modeling of the right hand side source terms in Eq. 7.2 can often cause numerical stability problems with the second moment closure models. The difficulty in deriving a turbulent dissipation equation was discussed in Chap. 6 and is a source of weakness in this approach of modeling turbulence as well a weakness in the application of two-equation turbulence models. Because of the numerical stability problems, the extra computational time requirements, and the limited improvement in solution accuracy over lower order models (Ref. 8), second moment closure is not commonly used for many complex CFD problems 7-5 today. Some cases where second moment closures have shown improvement over lower order models can be found in Ref. 9. 7.2 Algebraic Stress Models Under certain assumptions the Reynolds stress transport equations (Eq. 7.1, 7.2) can be reduced to a system of algebraic equations that require knowledge of the turbulent kinetic energy (k) and the turbulent dissipation (ε). This class of turbulence models is called an Algebraic Stress Model (ASM). ASM models reduce the closure problem to solving two transport equations and a system of algebraic equations. This is significantly faster than solving the full Reynolds stress transport equations set. The algebraic equations include most of the models and assumptions that were used to solve the full equation set. Rodi10 suggested that the transport of u i' u 'j is proportional to the transport of the turbulent kinetic energy k, the proportionality factor being the ratio u i' u 'j /k (which is not a constant). This yields Dui' u'j Dt − Diff ( ui' u'j ) = ' ' ui' u'j ⎛ Dk ⎞ ui u j (Pr oduction − ε ) − = Diff ( k ) ⎜ ⎟ k k ⎝ Dt ⎠ (7.19) Equation 7.19 is valid when the temporal and spatial change in u i' u 'j /k is small compared with the change of u i' u 'j itself. The algebraic stress model of Rodi is given by ⎡ ⎤ ⎢ (1 − C3 ) ⎛⎜ − u ' u ' ∂U i − u ' u ' ∂U i + 2 δ u ' u ' ∂U i ⎞⎟ ⎥ j l ij i j ⎢ ⎜ i l ∂x ∂xl 3 ∂x j ⎟⎠ ⎥ ε 2 l ⎝ ' ' ⎢ ⎥ ui u j = k δ ij + ⎢3 ⎥ ' ' ∂U i − ui u j ⎢ ⎥ ∂x j ⎢ ⎥ −1 C1 + ε ⎣⎢ ⎦⎥ (7.20) A more recent algebraic stress model derived by Abid, Morrison, Gatski, and Speziale11 (AMGS) is given by k k⎛ 1 1 ⎡⎛ ⎞⎤ ⎞ − ρ ui′u ′j = 2 μt* ⎢⎜ S ij − S kk δ ij ⎟ + α 4 (SikWkj + S jkWki ) − α 5 ⎜ Sik S kj − S kl S kl δ ij ⎟⎥ ε ε⎝ 3 3 ⎠⎦ (7.21) ⎠ ⎣⎝ 2 − ρkδ ij 3 where 7-6 μ t* = α 1 η= ξ= k ε k ε ( ) ( ) 3 1 + η 2 + 0.2 η 6 + ξ 6 k2 ρ 3 + η 2 + 6η 2ξ 2 + 6ξ 2 + η 6 + ξ 6 ε (7.22) (S S ) 1/ 2 ij (7.23) ij (W W ) 1/ 2 ij (7.24) ij ⎛4 ⎞g α1 = ⎜ − C 2 ⎟ ⎝3 ⎠2 g α 4 = (2 − C 4 ) 2 g2 4 α 2 = (2 − C3 ) 2 α 5 = (2 − C3 )g α 3 = (2 − C 4 ) 2 g2 4 (7.25) 1 g= (C1 / 2 ) + C5 − 1 The remaining coefficients for the model of Rodi and AMGS are given in Table 7.2. Model Rodi AMGS C1 1.5 6.8 C2 0.0 0.36 C3 0.6 1.25 C4 0.0 0.40 C5 0.0 1.88 Table 7.2 Coefficients for the ASM models of Rodi and AMGS Rodi’s model in Eq. 7.20 does not include the necessary terms to be valid near walls and requires wall functions (described in Chapter 10). The AMGS model is valid down to the wall. The solution obtained from an ASM model is not independent of the choice of model used to provide k and ε. The pressure coefficient on the Ames axisymmetric bump for a free stream Mach number of 0.875 is shown for the AMGS ASM model using a k-ω and k-ε two equation models for closure in Fig. 7.1. The velocity distribution at the trailing edge of the bump is shown in Fig. 7.2. Note the solutions are quite different for the same ASM model. 7-7 Figure 7.1 Pressure coefficient distribution for the axisymmetric bump at M=0.875. Figure 7.2 Velocity Distribution at the trailing edge of the bump for the axisymmetric bump at M=0.875. 7.3 Grid Sensitivity for a Flat Plate with Adiabatic Walls The initial wall spacing of the computational grid and the grid-stretching ratio can affect the accuracy any turbulence model. Figure 7.2 shows the sensitivity of the 7-8 skin friction to initial wall spacing for a flat plate for the SSG RSM model. Flat plate results for the AMGS ASM model are shown in Fig. 7.3. The grid-stretching ratio was 1.2 for all these cases. The plots include the theoretical skin friction curve of White. The boundary layer is seen to become fully turbulent around a length Reynolds number (Rex) of 1x106 for both models. These models are relatively insensitive to wall spacing below a y+=10. Figure 7.3 Flat plate skin friction predictions for the SSG RSM turbulence model for varying initial wall grid point spacings. 7-9 Figure 7.4 Flat plate skin friction predictions for the AMGS ASM turbulence model for varying initial wall grid point spacings. Predicted velocity profiles for the flat plate boundary layer for various initial wall grid point spacings are shown in Fig. 7.5 and 7.6. The velocity profile shows little effect of the initial spacing for all but the y+=10 profile for the ASM model. Figure 7.5 Flat plate boundary layer profiles predicted by the SSG RSM turbulence model for varying initial wall grid point spacings. 7-10 Figure 7.6 Flat plate boundary layer profiles predicted by the AMGS ASM turbulence model for varying initial wall grid point spacings. The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig. 7.7 and 7.8 All of these solutions used an initial wall spacing of y+=1. Figure 7.7 The effect of grid stretching ratio on the skin friction for a flat plate boundary layer using the SSG RSM turbulence model. 7-11 Figure 7.8 The effect of grid stretching ratio on the skin friction for a flat plate boundary layer using the AMGS ASM turbulence model. The RSM model shows little sensitivity to grid stretching below a ration of 2. The model has difficulty converging as can be seen by the spikes in the solution. The ASM solution shows sensitivity beginning at a ration of 1.5. The boundary layer profile predicted by both of these models show the same sensitivities to skin friction as the skin friction. This can be seen in Fig.7.9 and Fig. 7.10. Figure 7.9 The effect of grid stretching ratio on the velocity profile for a flat plate boundary layer using the SSG RSM turbulence model. 7-12 Figure 7.10 The effect of grid stretching ratio on the velocity profile for a flat plate boundary layer using the AMGS ASM turbulence model. 7.4 Grid Sensitivity for an Axisymmetric Bump A second example of the grid sensitivity of these turbulence models that includes a pressure gradient is the NASA Ames transonic axisymmetric bump experiment described in the Baldwin-Lomax Chapter 4. The effect of initial grid spacing on the pressure coefficient distribution along the bump for the RSM and ASM models is shown in Fig. 7.11 and Fig. 7.12. The stretching ratio was 1.2 for these cases. The pressure coefficient seems to be relatively insensitive to the initial grid spacing for the ASM model, with the y+=10 and y+=20 curves being slightly displaced from the other curves. The velocity distribution at the aft junction of the bump and the cylinder (x/c=1) are shown in Fig. 7.13 and 7.14. The trends are similar to those seen with the pressure coefficient . Grid stretching effects on the pressure coefficient distribution along the bump is shown in Fig. 7.15 and Fig. 7.16. The initial grid spacing was 1.2 for these cases. Both models show sensitivity to grid stretching ration in the separated flow region. The velocity profiles shown in Fig. 7.17 and 7.18 show a sensitivity to stretching ratio. The sensitivity of the Reynolds shear stress predicted by the SSG RSM model to both initial wall spacing and grid stretching ratio is shown in Fig. 7.19 and Fig. 7.20. The shear stress prediction is shown to be extremely sensitive to stretching ratio. 7-13 Figure 7.11 The effect of initial grid wall spacing on the pressure coefficient for the axisymmetric bump using the SSG RSM turbulence model. Figure 7.12 The effect of initial grid wall spacing on the pressure coefficient for the axisymmetric bump using AMGS ASM turbulence model. 7-14 Figure 7.13 The effect of initial grid wall spacing on the velocity distribution at x/c=1 for the axisymmetric bump using the SSG RSM turbulence model. Figure 7.14 The effect of initial grid wall spacing on the velocity distribution at x/c=1 for the axisymmetric bump using the AMGS ASM turbulence model. 7-15 Figure 7.15 The effect of grid stretching on the pressure coefficient for the axisymmetric bump using the SSG RMS turbulence model. Figure 7.16 The effect of grid stretching on the pressure coefficient for the axisymmetric bump using the AMGS ASM turbulence model. 7-16 Figure 7.17 The effect of grid stretching ratio on the velocity distribution at x/c=1 for the axisymmetric bump using the SSG RSM turbulence model. Figure 7.18 The effect of grid stretching ratio on the velocity distribution at x/c=1 for the axisymmetric bump using the AMGS ASM turbulence model. 7-17 Figure 7.19 The effect of initial wall spacing on the shear stress distribution at x/c=1 for the axisymmetric bump using the SSG RSM turbulence model. Figure 7.20 The effect of grid stretching ratio on the shear stress distribution at x/c=1 for the axisymmetric bump using the SSG RSM turbulence model. 7-18 RSM and ASM Model Application Hints 1. 2. 3. RSM models generally suffer from numerical stability issues. This is due to the complexity of the modeled terms. This is the primary reason that these models are not regularly used in large scale applications. ASM models are sensitive to the transport equations used to obtain k and ε. These models are generally less numerical stable than traditional two-equation models. The choice of the pressure strain term is critical to the performance of both RSM and ASM models. More complex pressure strain models are generally more accurate but less numerically stable. Chapter 7 References: 1. Daly, B. and Harlow, F., “Transport Equations of Turbulence,” Physics of Fluids, 1970, pp. 2634-2649. 2. Hanjalic, K. and Launder, B., “Contribution Towards a Reynolds-Stress Closure for Low-Reynolds-Number Turbulence,” Journal of Fluid Mechanics, Vol. 74, 1976, pp. 593-610. 3. Mellor, G. and Herring, H., “A Survey of the Mean Turbulent Field Closure Models,” AIAA Journal, Vol. 11, 1973, pp. 590-599. 4. Launder, B., Reece, G., and Rodi, W., “Progress in the development of a Reynolds-Stress Turbulence Closure,” Journal of Fluid Mechanics, 1975, Vol. 68, Part 3, pp. 537-566. 5. Gibson, M. and Launder, B., “On the Calculation of Horizontal, Turbulent, Free Shear Flows under Gravitational Influence,” Journal of Heat Transfer, Vol. 98, 1976, pp. 81-87. 6. Speziale, C., Sarker, S., and Gibson, T., “Modeling the Pressure-Strain Correlation of Turbulence: An Invariant Dynamical Systems Approach,” Journal of Fluid Mechanics, Vol. 227, 1991, pp. 245-272. 7. Schumann, U., “Realizability of Reynolds Stress Turbulence Models,” Physics of Fluids, Vol. 20, pp. 721-725, 1977. 8. Leschziner, M., Batten, P. and Craft, T., “Reynolds-Stress Modeling of Transonic Afterbody Flows,” The Aeronautical Journal, June 2001, pp. 297306. 9. Nallasamy, M., “Turbulence Models and Their Applications to the Prediction of Internal Flows: A Review,” Computers and Fluids, Vol. 15, pp. 151-194, 1987. 10. Rodi, W., “Turbulence Models and Their Application in Hydraulics – A State of the Art Review,” International Association for Hydraulic Research, 2nd Ed., Feb. 1984. 11. Abid, R., Morrison, J., Gatski, T., and Speziale, C., “Prediction of Aerodynamic Flows with a New Explicit Algebraic Stress Model,” AIAA Journal, Vol. 34, No. 12, Dec. 1996. 8-1 8.0 Large Eddy Simulation Theoretically the Navier-Stokes equations can be used to simulate turbulent flows. The computational grid used in such a simulation would have to be fine enough to allow the smallest turbulent length scales to be realized and the computational time step would have to be small enough to simulate the highest frequencies of the turbulent spectrum. In practice such computations are prohibitively expensive for high Reynolds numbers. The Reynolds averaged Navier-Stokes (RANS) equations were derived assuming that all of the unsteadiness due to the turbulent nature of the flow could be modeled with empirically derived correlations. This reduces the time and length scales that must be simulated, but also limits the applicability of the simulation for unsteady flows. Large eddy simulations (LES) were developed to extend the simulation of unsteady flows beyond DNS. The desired result of an LES computation is to obtain a DNS equivalent solution for the large-scale turbulence on a much coarser grid than is required for DNS. An LES simulation requires: 1. A grid fine enough to discretize the small nearly isotropic scales of the turbulence A low dissipation numerical scheme A filter function to determine the division of the turbulent spectrum into grid realized and subgrid regions A subgrid turbulence model 2. 3. 4. A true LES simulation is more than a high Reynolds number computation run without a turbulence model. Although the resulting solution from such a simulation may resemble turbulent flow, the resulting solution will most likely not represent an equivalent DNS solution. 8.1 The Filtering Operation The filtering operation is critical to LES. Consider a filtering operation with a uniform characteristic filter width Δ (which implies isotropic grid elements). Leonard1 defined the following filter in physical space ∞ φ ( x ) = ∫ G ( x − ξ )φ (ξ )dξ −∞ (8.1) Note that in this chapter the over bar represents a filtered quantity, not a time averaged quantity as in the previous chapters. The filtering operation is a spatial operation as opposed to the Reynolds averaging operation discussed in Chapter 2 that is a temporal operation. The original function φ is then decomposed into a filtered field (or grid resolved) and a “subgrid” term φ = φ +φ′ (8.2) 8-2 The function G defined in Eq. 8.1 is the “filter function”. The filter function may be any function defined on an infinite domain that satisfies the following requirements (Ref. 1): 1. G ( −ξ ) = G (ξ ) ∞ 2. ∫ G(ξ )dξ = 1 −∞ 3. G (ξ ) → 0 as ξ → 0 such that all moments ∞ ∫ G(ξ )ξ n dξ (n ≥ 0) exist −∞ ⎛ Δ Δ⎞ 4. G (ξ ) is “small” outside ⎜ − , ⎟ ⎝ 2 2⎠ One important and useful feature of this choice of filter is ∂φ ∂φ = ∂x ∂x (8.3) Three common filters are ⎛ πξ ⎞ sin⎜ ⎟ ⎝ Δ ⎠ Spectral cutoff filter: G (ξ ) = πξ (8.4) Δ Gaussian filter: G (ξ ) = ⎛ Cξ 2 ⎞ 1 C exp⎜⎜ 2 ⎟⎟ Δ π ⎝ Δ ⎠ ⎧1 Δ ⎪⎪ Δ if ξ < 2 Box filter: G (ξ ) = ⎨ ⎪ ⎩⎪ 0 otherwise (8.5) (8.6) The spectral cutoff filter is normally applied in spectral space as ϕˆ (k ) = Gˆ (k )ϕˆ ( k ) (8.8) 8-3 2π ⎧ ⎪⎪1 if k < Δ Gˆ (ξ ) = ⎨ ⎪ ⎪⎩ 0 otherwise (8.9) The constant C in the Gaussian filter (Eq. 8.5) is somewhat arbitrary, and values from 2 to 6 have been used in practice. Note that the box filter (Eq. 8.6) is valid on both finite and infinite domains. For this reason box filters are often used to relate DNS to LES for physical space numerical schemes. There are advantages to using spectral and pseudo-spectral solution methodologies for LES. The first advantage is that the approximate field, which is discrete in spectral space, is a finite sum of continuous functions in physical space. Thus the approximate field, the spectral filter, and their derivatives are continuous functions in physical space. The spectral filters defined in Eqs. 8.8 and 8.9 have the following additional useful properties φ =φ (8.10) φ′ = 0 (8.11) This says that filtering a grid resolved quantity yields the original grid resolved quantity and that the spectral filter of a subgrid quantity is equal to zero. Another advantage of the spectral approach is that derivatives yield an exact value for the approximate field. In other words, if one is approximating a function using N Fourier modes, then as long as one has at least 2(N+1) points in physical space, it does not matter if one has 20 or 20 million points: the value of the derivative at any given point is the same. Therefore one can take a flow field generated with a spectral DNS code, filter and coarsen the solution, and produce a flow field which satisfies the governing equations for LES. There are drawbacks to spectral filters. First, they have limited applicability in realistic flow situations because of boundary constraints. Also, spectral filters are non-positive (i.e. the filter function is negative at some points in space). The latter condition results in a subgrid stress tensor that does not satisfy the Reynolds stress realizability conditions outlined in Chapter 2. Both Gaussian filters (Eq. 8.5) and box filters (Eq. 8.6) are positive functions, and thus the Reynolds stress realizability conditions outlined in Chapter 2 will be satisfied on the resulting subgrid stress tensor. Unfortunately Eq. 8.10 and Eq. 8.11 are not valid for these filters. Furthermore, if a spatial numerical solver is to be used the filter must also be discretized. Therefore it is necessary to examine the discretized system that is solved in the actual code because discrete systems 8-4 sometimes have different properties than the original continuous systems from which they were derived. 8.2 Derivation of the LES Equations The Favre averaged or filtered Navier-Stokes equations were derived in Chapter 2. The Favre filter can be defined for any variable as ~ φ = ρφ ρ (8.12) Note that the over bar signifies a spatially filtered quantity, and the tilde represents a Favre filtered, or grid resolved, quantity. Thus the Favre filter may be thought of as a density weighted filter in space. Applying this filtering operation to the Navier-Stokes equations, and assuming that the filtering commutes with the derivative operation, the LES equations are ∂ρ ∂ρ u~i + =0 ∂t ∂xi (8.13) ~ ⎞ ⎛ ~ ∂ρ u~i ∂ρ u~i u~ j ∂p ∂τ ij ∂ (τ ij − τ~ij ) − ∂ ρ ⎜⎜ ui u j − u~i u~ j ⎟⎟ + =− + + ∂t ∂x j ∂xi ∂x j ∂x j ∂x j ⎠ ⎝ (8.14) ~ ∂ρ E ∂ ∂ ⎛⎜ ~ ~~ ⎞⎟ ∂ ~ ρ E + p ui = − ρ Eui − Eui − pui − pu~i + ⎟ ∂xi ∂t ∂x i ∂xi ⎜ ⎝ ⎠ ~ ~ ~ ~ ∂u jτ ij ∂ ∂ ⎛ ~ ∂T ⎞ ∂ ⎛ ∂T ~ ∂T ⎞ ⎟ ⎜⎜ κ ⎟⎟ + ⎜⎜ κ u jτ ij − u~ jτ~ij + −κ + + ∂xi ⎟⎠ ∂xi ∂xi ∂xi ⎝ ∂xi ⎠ ∂xi ⎝ ∂xi (8.15) ( ( ) ( ) ) The resolved viscous stress tensor in Eqn. 8.14 takes the form ⎡ 1 ⎛ ∂u~i ∂u~ j ⎞ 1 ∂u~k ⎤ ~ ~ ⎟− τ ij = 2 μ ⎢ ⎜⎜ δ + ⎟ 3 ∂x ij ⎥ 2 x x ∂ ∂ i ⎠ k ⎢⎣ ⎝ j ⎥⎦ (8.16) The viscosity and thermal conductivity is assumed to be calculated from the ~ ~ Favre averaged temperature ( T ). The total energy ( E ) is defined as 1 ~ E=~ e + u~i u~ j + k sgs 2 (8.17) The subgrid kinetic energy, ksgs, is defined as the effect of the subgrid scales on the kinetic energy of the resolved field 8-5 k sgs = 1 ⎛⎜ ~ ~ ~ ⎞⎟ ui u i − u i ui ⎟ 2⎜ ⎠ ⎝ (8.18) Finally, the LES equation of state is ~ p = ρ RT (8.19) As with Reynolds averaging, the filtering operation produces terms that must be modeled in order to close the equation set. The LES momentum equation, Eqn. ∂ 8.14, contains two such terms. The first, (τ ij − τ~ij ) , represents the difference ∂x j in the viscous terms between Favre and “straight” filtering. The second term, ∂ ⎛⎜ ~ ~ ~ ⎞⎟ ρ ui u j − ui u j , arises from the convection terms and has been observed to ⎟ ∂x j ⎜ ⎠ ⎝ behave much like a viscous stress. Thus it is frequently redefined in terms of the “subgrid stress tensor” τ ijsgs defined as ⎛ ~ ⎞ ⎜ ⎝ ⎟ ⎠ τ ijsgs = ρ ⎜ ui u j − u~i u~ j ⎟ (8.20) The LES energy equation, Eq. 8.15, contains four subgrid terms. The first term, ∂ ⎛⎜ ~ ~~ ⎞⎟ ρ Eui − Eui , arises from the convection term. The second term is a ⎟ ∂xi ⎜ ⎝ ⎠ ∂ pui − pu~i , which is also related to the pressure-velocity correlation term, ∂xi ∂ u jτ ij − u~ jτ~ij , is a viscous subgrid term convection term. The third term, ∂xi similar to the term in the LES momentum equation. This term represents the transfer of energy due to subgrid viscous forces. The fourth term, ~ ∂ ⎛ ∂T ~ ∂T ⎞ ⎜⎜ κ ⎟⎟ , is a heat flux subgrid term. The subgrid kinetic energy, ksgs, −κ ∂xi ⎝ ∂xi ∂xi ⎠ used in the definition of the Favre filtered total energy in Eq. 8.7 must also be defined for closure. ( ( ) ) The subgrid stress tensor, τ ijsgs is usually treated in a similar manner as the Boussinesq approximation for the RANS equations (Chapter 2) 8-6 ⎛~ ⎝ 1~ 3 ⎞ ⎠ τ ijsgs = −2 ρ ν t ⎜ S ij − S kk δ ij ⎟ + 2 sgs ρ k δ ij 3 (8.21) ~ In the above equation, S ij represents the resolved rate of strain tensor ~ 1 ⎛ ∂u~ ∂u j ~ S ij = ⎜ i + 2 ⎜⎝ ∂x j ∂xi ⎞ ⎟ ⎟ ⎠ (8.22) Note that a subgrid eddy viscosity, νt, has been introduced in Eq. 8.21. The subgrid eddy viscosity accounts for the turbulence that cannot be resolved on the computational grid. Introducing the subgrid stress tensor into Eq. 8.14 and assuming that the straight filtered stress tensor is equal to the Favre filtered stress tensor (τ ij = τ~ij ) the LES momentum equation becomes sgs ~ ∂ρ u~i ∂ρ u~i u~ j ∂p ∂τ ij ∂τ ij + =− + − ∂t ∂x j ∂x i ∂x j ∂x j (8.23) The LES energy equation (Eq. 8.15) can be further simplified as follows. Using the definition of total enthalpy H =E+ p (8.24) ρ two of the subgrid terms can be combined to form ∂ − ∂xi ~ ⎛ ⎞ ~ ~~ ⎞⎟ ⎜ ρ Eu − ρ E~u~ ⎟ − ∂ pu − pu~ = ∂ ⎛⎜ ρ Hu − ρ H ui i i i i i ⎟ ∂xi ⎜ ⎟ ∂xi ⎜ ⎝ ⎠ ⎝ ⎠ ( ) (8.25) This term can then be modeled with an eddy diffusion model as ∂ ∂xi ~ ⎛ ~ ∂ ⎛ ∂H ⎞ ~~ ⎞⎟ sgs ⎜ ρ Hu − ρ H ⎜ ⎟ −c ρ k Δ ui ≈ i ⎜ ⎟ ∂xi ⎜⎝ e ∂xi ⎟⎠ ⎝ ⎠ (8.26) where Δ is the local grid spacing and ce is an empirical coefficient. The thermal ~ ∂ ⎛ ∂T ~ ∂T ⎞ ⎜⎜ κ ⎟ , is neglected based on the assumption that conductivity term, −κ ∂xi ⎝ ∂xi ∂xi ⎟⎠ the thermal conductivity is locally constant and that the Favre filtering is equivalent to the straight filtering. The viscous transport term is modeled as 8-7 ( ) ∂ ∂ ⎛ ∂k sgs ⎞ ⎜μ ⎟ u jτ ij − u~ jτ~ij ≈ ∂xi ∂xi ⎜⎝ ∂xi ⎟⎠ (8.27) This term is generally small and is often neglected in practice. Thus the modeled LES energy equation may now be written as ~ ∂ρ E ∂ ∂ ~ + ρ E + p ui = ∂t ∂xi ∂xi ~ ∂ ⎛ ~ ∂k sgs ⎞ ∂ ⎛ ~ ∂T ⎞ ⎜κ ⎟ ⎜μ ⎟+ ∂xi ⎜⎝ ∂xi ⎟⎠ ∂xi ⎜⎝ ∂xi ⎟⎠ ( ) ~ ~~ ⎛ ∂H ⎞ ∂u jτ ij ⎜⎜ ce ρ k sgs Δ ⎟⎟ + + ∂xi ⎠ ∂xi ⎝ (8.28) Inspection of the modeled LES equations (Eqs. 8.13, 8.23, and 8.28) shows that these equations are quite similar to the RANS equations incorporating the Boussinesq approximation derived in Chapter 2. It should be noted that the time averaging process used to derive the RANS equations and the spatial filtering operations used to derive the LES equations are quite different, and hence the terms in the two equation sets are not the same. Three things are required to close the modeled LES equations: 1. The subgrid turbulent kinetic energy ksgs 2. The subgrid eddy viscosity νt 3. The empirical coefficient in the energy equation, ce 8.3 Smagorinsky Model As with the RANS equations, there are numerous models for the LES eddy viscosity. One of the earliest models for the LES eddy viscosity was proposed by Smagorinski2 and is given by ν t = (C s Δ g )2 (2 Sij S ij ) 2 ~ ~ 1 (8.29) where Cs is the Smagorinsky coefficient and Δg is the local grid spacing. This model ignores the subgrid kinetic energy ksgs. The Smagorinsky model provides a simple closure for the LES equations and has been used effectively for a number of applications. But, like its algebraic counterpart in the RANS regime, it has been shown to be lacking in simulating complex turbulent flows. The Smagorinsky model is not valid near walls, and a wall damping term (Ref. 3) is often added to the model. The wall-damped form is given by 8-8 ( ( ν t = Cs Δ g 1 − e − y + / 25 )) (2S~ S~ ) 1 2 ij 2 (8.30) ij The coefficient Cs must be defined for this model. In practice no universal value of Cs exist. 8.4 Dynamic Smagorinsky Model Several investigators have attempted to dynamically calculate Cs based on equating the highest wave number resolved turbulent stresses with the subgrid stress. One dynamic LES closure model outlined in Ref. 4 can be described as follows. The subgrid stress tensor is written as τ ijsgs = −C (x, t )Δ2 S (u~ ) S ij (u~ ) (8.31) where C(x,t) is the Smagorinsky coefficient to be determined dynamically and Δ is the grid spacing. For this purpose a second spatial filter, called the test-filter, of width larger than the grid filter is introduced. We choose the test-filter scale Δ = 2Δ . This filter generates a second set of resolvable-scale fields (denoted by ( ) ). Then, the dynamic SGS model of the τ ij at the grid filter and of the Tij at the test-filter are written as: δ ij ⎧ 2 ⎪τ ij − τ kk = −2C Δ S S ij 3 ⎪ = −2Cβ ij ⎪ ⎨ ⎪T − δ ij T = −2C Δ 2 S S ij ⎪ ij 3 kk ⎪ = −2Cα ij ⎩ (8.32) A least squares method is developed to predict C(x,t) as: C ( x, t ) = − ( ) l ij α ij − β ij 1 2 α −β α −β mn mn mn mn ( )( ) (8.33) with l aij = l ij − (δ ij / 3)l kk = −2Cα ij + 2C β ij . 8.5 k-Equation Model Another approach to closing the LES equations is to introduce a transport equation for the subgrid kinetic energy, ksgs. The exact subgrid kinetic energy equation may be written 8-9 ⎛ ∂p ~ ∂p ⎞ ⎛ ∂τ ij ~ ∂τ ij ⎞ ∂ρ k sgs ∂ ⎟− ⎟⎟ + ⎜ ui ( + − ui − ui ρ k sgs u~i ) = −⎜⎜ ui ⎜ ⎟ ∂t ∂xi ∂ ∂ ∂ ∂ x x x x i i ⎠ i i ⎠ ⎝ ⎝ ∂u~ jτ ijsgs ∂u~ j ~~ ∂ sgs ρkui − ρ k ui − + τ ij ∂xi ∂xi ∂xi ( ) (8.34) This equation contains several subgrid terms. The equation is modeled in a similar manner as the RANS turbulent kinetic energy equation. The modeled equation is 3 ∂u~ j ( ∂ρ k sgs ∂ ∂ ⎛ ∂k sgs ⎞ k sgs ) 2 sgs ~ sgs (ρk ui ) = ∂x ⎜⎜ ρν t ∂x ⎟⎟ − τ ij ∂x − ρcε Δ + ∂t ∂xi i ⎝ i ⎠ i (8.35) where the terms on the right-hand-side are the diffusion, production, and dissipation of the subgrid kinetic energy respectively. An algebraic expression for the dissipation is generally used instead of a transport equation as in twoequation RANS turbulence models. The turbulent viscosity is given by ν t = cν k sgs Δ (8.36) The LES energy equation dissipation coefficient cε and the eddy viscosity coefficient cν are determined dynamically by equating the highest wave number resolved turbulent properties with the subgrid properties. The dynamic process for determining these coefficients is described in detail in Ref. 5. 8.6 Inflow Turbulence Boundary Condition Flows that include a developing boundary layer at the inflow plane of a simulation require a boundary condition that includes the specification of the inflow velocity fluctuations. Soteriou and Ghoniem6 show the development of an incompressible mixing layer with and without perturbation at the inflow plane. Compared to the perturbed case, the development of the mixing layer without perturbation is delayed significantly. This is because disturbances within the boundary layer are amplified and lead to the Kelvin-Helmholtz instability that eventually causes the shear layer to roll up. If these disturbances are not present in the boundary layer, then the shear layer will develop more slowly until disturbances are generated through minute numerical errors. Several approaches have been taken to add the perturbations to the inflow. The simplest approach is to add “white noise” to the velocity at the inflow plane as was done by Comte et al7. This approach is somewhat unphysical in that the perturbations have no correlation in space or time. A second approach is to add perturbations at discreet frequencies as was done by Soteriou and Ghoniem6. This can be effective for studying forced flows, but it must be remembered that turbulence is broadband and does not confine itself to discrete frequencies. A 8-10 third approach to this problem is to attempt to reconstruct the unsteady flow at the inflow using correlations for the perturbations. One example of this can be found in Ref. 5, in which a box of frozen turbulence is generated and saved. This box is scaled and applied as perturbations at the inflow of the simulation. 8.6 Other LES References LES methods are still evolving as researches investigate methods for simulating the near-wall flow physics and methods to improve the performance of the subgrid models used in LES simulations. Fureby, et al8, provide a good summary of the LES methods currently being employed and investigated. Ref. 9 includes a wide range of experimental test cases for validation of LES simulations. 8.7 Spatial Mixing-Layer Example The mixing layer experiment of Samimy and Elliot8 can be used to demonstrate the behavior of LES models. The data is also included in Ref. 10. The experimental setup had an upper supersonic stream and a lower subsonic stream mixing at a matched static pressure. Laser Doppler Velocimeter (LDV) measurements were made at several downstream locations between the trailing edge of the splitter plate and a station 210 mm downstream of the splitter plate. The flow parameters are given in Table 8.1. T0, K P01, kPa M1 M2 Mc U1, m/sec U2/U1 ρ2/ρ1 δ1, mm 291.0 314.0 1.80 0.51 0.52 479.5 0.355 0.638 8.0 Table 8.1 Flow parameters for the spatial mixing-layer case The solutions were run for 30,000 time steps to wash out initial transients and then statistics were taken over the next 60,000 time steps. The time step used was 5.0x10-8 seconds. The flow solver used for these calculations was 4th order in space and 2nd order in time. The computational grid dimensions were 181 x 121 x 61. The inflow plane perturbations were specified using the “box of turbulence” approach from Ref. 5. Further details of the experiment and the computations can be found in Ref. 10. All of the turbulence models investigated predict roughly the same mixing-layer thickness for this case as shown in Fig 8.1. The Hybrid k-ε model in Fig. 8.1 is the multi-scale hybrid RANS/LES turbulence model described in Chapter 9. The DES model is the Spalart-Allmaras DES hybrid model that is also described in Chapter 9. The traditional RANS models (k-ε and Spalart-Allmaras) have the worst agreement with the experimental data. 8-11 Figure 8.1 Mixing-layer thickness The streamwise velocity 210 mm downstream of the splitter plate is shown in Fig. 8.2. Again all of the models are in reasonable agreement with the experimental data. Figure 8.2 Nondimensional streamwise velocity at x=210 mm 8-12 The streamwise turbulence intensity, Reynolds stress, streamwise velocity fluctuation skewness, and streamwise velocity fluctuation flatness at x=210 mm are shown in Fig. 8.3-8.6. The performance of the traditional RANS models is seen to deteriorate as the order of the turbulent velocity fluctuation correlation increases. The LES and hybrid RANS/LES models are all in reasonable agreement with the experimental data. Figure 8.3 Streamwise turbulence intensity at x=210 mm Figure 8.4 Reynolds stress profiles at x=210 mm 8-13 Figure 8.5 Stream wise velocity fluctuation skewness at x=210 mm Figure 8.6 Streamwise velocity fluctuation flatness at x=210 mm. Chapter 8 References: 1. Leonard, A., “Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows,” Advances in Geophysics, Vol. 18A, pp. 237-348, 1974. 2. Smagorinsky, J. S., “General Circulation Experiments with the Primitive Equations I: The Basic Experiment,” Monthly Weather Review, Vol. 91, pp. 99-164, 1963. 3. Piomelli, U., Zang, T. A., Speziale, C. G., and Hussaini, M. Y., Physics of Fluids, Vol. 2, pp. 257, 1990. 4. Wasistho, B., Haselbacher, A., Najjar, F.M., Tafti, D., Balachandar, S. and Moser, R.D., “Direct and Large Eddy Simulations of Compressible Wall- 8-14 Injection Flows in Laminar, Transitional and Turbulent Regimes,” AIAA 2002-4344, 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, July 7-10, 2002, Indianapolis, IN. 5. Nelson, C. C., “Simulations of Spatially Resolved Compressible Turbulence Using a Local Dynamic Subgrid Model,” Ph.D. Thesis, Georgia Institute of Technology, 1997. 6. Soteriou, M. and Ghoniem, A., “Effects of the Free-Stream Density Ratio on Free and Forced Spatially Developing Shear Layers,” Physics of Fluids, Vol. 7, No. 8, 1995, pp. 2036-2051. 7. Comte, P., Lesieur, M., Laroche, H., and Normand, X., “Numerical Simulations of Turbulent Plane Shear Layers,” in Turbulent Shear Flows 6, Springer-Velag, 1989, pp. 360-380. 8. Fureby, C., Alin,N., Wikstrom, N., Menon, S., Svanstedt, N. and Persson, L., “Large-Eddy Simulation of High-Reynolds-Number Wall-Bounded Flows,” AIAA Journal, Vol. 42, No. 3, March 2004, pp 457-468. 9. “A Selection of Test Cases for the Validation of Large-Eddy Simulations of Turbulent Flows,” AGARD-AR-345, April 1998. 10. Samimy, M. and Elliot, G., “Effects of Compressibility on the Characteristics of Free Shear Layers,” AIAA Journal, Vol. 28, No. 3, 1990, pp.439-445. 11. Nelson, C. and Nichols, R., “Evaluation of Hybrid RANS/LES Turbulence Models Using an LES Code,” AIAA-2003-3552, June 2003. 9-1 9.0 Hybrid RANS/LES Models 9.1 Theory Several investigators have noted a limitation with RANS turbulence models when applied to unsteady flows. The turbulence models produce too much eddy viscosity and over-damp the unsteady motion of the fluid. The problem is inherent in the construction of the turbulence models because of the assumption that all scales of the unsteady motion of the fluid are to be captured and modeled by the turbulence model. One approach to overcoming the shortcomings of RANS models applied to unsteady flows is to spatially filter the RANS turbulence models such that the eddy viscosity does not include the energy of grid resolved turbulent scales. The spatially filtered RANS turbulence model may be thought of as a subgrid model for very large turbulent eddies. This class of turbulence models has been called hybrid RANS/LES models because they incorporate aspects of both forms of turbulence modeling. It is desirable that the spatial filter function chosen not degrade the performance of the turbulence model when the largest turbulent scales present are below the resolution of the grid as is often the case in current aircraft CFD applications. There are two possibilities for developing hybrid models from existing RANS models. A first approach would be to modify the production or dissipation source terms for existing turbulence model differential equations to include additional terms to adjust the local eddy viscosity so that it does not include the grid realized contribution. This means that the turbulence quantities normally transported by the RANS turbulence model will be created or destroyed based on the local grid resolution. A second approach would be to solve the existing RANS turbulence model in the normal manner and then filter the results to determine the level of eddy viscosity that will be used in the solution of the Navier-Stokes equations. This approach is somewhat simpler to develop since it requires no tuning of the differential transport equations and can be easily extended to different RANS turbulence models. Both approaches have been applied to unsteady flows. One example of the first form of a hybrid RANS/LES turbulence model based on the Spalart-Allmaras one-equation turbulence model can be found in Ref. 1. The standard SpalartAllmaras turbulence model contains a destruction term for eddy viscosity that is inversely proportional to the square of the distance from the wall (d). Ref. 1 suggests replacing the wall distance (d) in the destruction term with ~ d = min (d , CDES Lg ) (9.1) 9-2 where CDES is a constant of O(1) and Lg is a grid length scale defined by Lg = max(Δx, Δy, Δz ) (9.2) where Δx, Δy, and Δz are the local grid lengths. For two-dimensional or axisymmetric problems, the out of plane direction is simply ignored in defining the local grid length. The modified destruction term has the effect of decreasing the eddy viscosity in regions of tight grid spacing. This modification causes the Spalart-Allmaras RANS turbulence model to behave like a Smagorinsky LES turbulence model (see Chapter 8) when the grid spacing (Lg) is less than the distance from the wall, which is generally the case outside of the boundary layer. Note that the transition from RANS to LES does not include any turbulent length scale dependence, but is solely a function of the local grid spacing. Ref. 1 introduces the term Detached-Eddy Simulation (DES) to describe this model. DES has been applied to a number of unsteady flow problems including flow over a sphere2, flow over a delta wing3 and flow over an aircraft4. A similar modification to an SST two-equation model is given by Strelets5. In this hybrid model the dissipation term (ε) in the turbulent kinetic energy equation (k) is replaced by 3 ω k 2 ε= * = = β k min (Lt , C DES L g ) ε Lg ⎞ ⎛ ⎟⎟ min⎜⎜1.0, C DES L t ⎠ ⎝ (9.3) where Lt = k β *ω = k 3 2 ε C DES = (1 − F1 )C DESKE + F1C DESKW (9.4) (9.5) here ω is the specific dissipation, β *=0.09, CDESKE=0.61, and CDESKW=0.78. The turbulent dissipation ε defined in Eq. 9.3 is effectively increased when the grid size length scale Lg (Eq. 9.2) is less than the turbulent length scale Lt (Eq. 9.4). This causes the eddy viscosity to be reduced in these regions. Unlike the oneequation SA-DES model, this model does include a dependency on the local turbulent length scale. The turbulent length scale used in this model is a function of the filtered (subgrid) turbulent quantities. This model behaves like a kequation LES subgrid model when the turbulent length scale is greater than the grid length scale, and the dissipation equation is decoupled from the kinetic energy equation in this region. 9-3 Nichols and Nelson6 give an example of the second approach for developing a hybrid RANS/LES turbulence model which they have designated a multi-scale model. The method was implemented in conjunction with the SST two-equation turbulence model and is termed a multi-scale model. The turbulent length scale used in this effort is defined by ( 3/ 2 Lt = max 6.0 ν tRANS / Ω , k RANS / ε RANS ) (9.6) where νtRANS is the unfiltered (large scale) eddy viscosity and Ω is the local mean flow vorticity. This length scale is a mixture of the traditional turbulent scale 3/ 2 /εRANS) and the definition definition for two-equation turbulence models ( k RANS usually associated with algebraic turbulence models ((νtRANS/Ω)1/2). The turbulent length scale definition could be easily adapted to other types of turbulence models. The subgrid turbulent kinetic energy is defined as k LES = k RANS f d (9.7) The damping function is defined as f d = (1 + tanh (2π (Λ − 0.5))) / 2 (9.8) where Λ= 1.0 ⎛L 1.0 + ⎜ t ⎜ Lg ⎝ ⎞ ⎟ ⎟ ⎠ 4/3 (9.9) where Lg is defined in Eq. 9.2. The eddy viscosity is then calculated from ν t = ν tRANS f d + (1 − f d )ν tLES (9.10) where the LES based subgrid eddy viscosity is given by ν tLES = min(0.0854 Lg k LES ,ν tRANS ) (9.11) The multi-scale hybrid model behaves like a traditional SST model on the RANS end of the spectrum and transitions to a nonlinear k-equation model on the LES end of the spectrum. The transition function in Eq. 9.10 was chosen to allow a smooth transition from the standard RANS turbulence model to the LES subgrid model. This hybrid RANS/LES approach can easily be extended to other RANS turbulence models with little to no alteration providing that a length scale can be derived for that model and a value for the turbulent kinetic energy can be approximated by that model. The multi-scale hybrid model, like the SST-DES 9-4 model, transitions from RANS to LES as a function of the ratio of the local turbulent length scale predicted by the RANS model and the local grid spacing rather than being a function of the grid spacing alone as is the case for the Spalart-Allmaras DES model. Nelson and Nichols evaluated the Spalart-Allmaras DES and the multi-scale hybrid turbulence models as subgrid turbulence models for LES applications in Ref. 7. The hybrid models were applied to a high-speed shear layer using a true LES flow solver. The results indicate that the hybrid models perform as well as more complicated LES subgrid models for this application. Nichols8 investigated these three hybrid RANS/LES models applied to a circular cylinder and a generic weapons bay. These results are shown in the following sections. 9.2 Circular Cylinder Unsteady three-dimensional calculations were performed for the vortex shedding from a circular cylinder for M=0.2 and Red=8x106. Three computational grids were used in the simulation: Fine – 401x201x201 Mid – 201x101x101 Coarse – 101x51x51 The mid and coarse level grids were constructed by removing every other point from the fine and mid level grids respectively. The grids have a span of ten cylinder diameters. Periodic boundary conditions were applied at the side planes of the grid. The initial wall spacing was 2x10-4 diameters for the fine grid, which corresponds to a y+ of about 20. Wall function boundary conditions were used in the calculations. Calculations with three different time steps were run on the mid grid with the multi-scale hybrid turbulence model to determine an acceptable time step for this simulation. The results are shown in Fig. 9.1. The primary shedding frequency for this example is just less than 60 Hz (Strouhal number=0.28). These solutions correspond to 50, 200, and 800 time steps per shedding cycle. The results from the largest time scale are quite different from the smaller time step results. The spectral peak predicted with the large time step occurs at a much lower frequency than the peak predicted using the smaller time steps. Based on these results all subsequent simulations were performed using the 9.0E-4 second time step with 200 time steps per shedding cycle. 9-5 Figure 9.1. Normal force power spectral density for 3 different time steps using the SST-MS model on the medium grid. The data in Ref. 9 was obtained using both air (γ=1.4) and Freon (γ=1.13) as a test medium. Calculations were performed on the medium grid with the SST-MS hybrid model to assess the affect of test medium. The normal force spectral results are shown in Fig. 9.2. This indicates that the test medium is not an important factor at the conditions of this study. All subsequent solutions were performed using air as the test medium. Figure 9.2. Normal force power spectral density for air and Freon using the SSTMS hybrid model. 9-6 The calculations were run 10,000 iterations and the final 4096 time steps were statistically analyzed. The centerline plane of the three grids is shown in Fig. 9.3. Fine Grid Medium Grid Coarse Grid Figure 9.3. Cylinder grid centerline. A comparison of the normal force power spectral density for the three hybrid models, the Spalart-Allmaras model, and the SST model for the medium grid are shown in Fig. 9.4. The hybrid models show similar trends in that they all have a rather broad spectral peak at a Strouhal number of about 0.28. The two RANS models have a much narrower spectral peak at a similar Strouhal number. The RANS models have much lower energy away from the peak. This is an indication that the RANS models are providing too much damping of the unsteady solutions. Figure 9.4. Normal force power spectral density on the medium grid. 9-7 Instantaneous x-vorticity isosurfaces computed with the SST-DES hybrid turbulence model are shown in Fig. 9.5 for all three grids. The turbulent structure in the wake of the cylinder is clearly evident in the fine grid solution. The turbulent structure is reduced in the mid grid solution. There is very little structure present in the coarse grid solution. Fine Grid Medium Grid Coarse Grid Figure 9.5. Instantaneous vorticity isosurfaces colored by Mach number for the SST-DES hybrid turbulence model. Instantaneous contours of Mach number and eddy viscosity for the SA-DES, SST-DES, and SST-MS hybrid turbulence models computed on the three grid levels are shown in Figs. 9.6-9.11. a. Fine b. Mid c. Coarse Figure 9.6. Instantaneous Mach number contours for the SA-DES hybrid turbulence model. 9-8 a. Fine b. Mid c. Coarse Figure 9.7. Instantaneous eddy viscosity contours for the SA-DES hybrid turbulence model. a. Fine b. Mid c. Coarse Figure 9.8. Instantaneous Mach number contours for the SST-DES hybrid turbulence model. a. Fine b. Mid c. Coarse Figure 9.9. Instantaneous eddy viscosity contours for the SST-DES hybrid turbulence model. 9-9 a. Fine b. Mid c. Coarse Figure 9.10. Instantaneous Mach number contours for the SST-MS hybrid turbulence model. a. Fine b. Mid c. Coarse Figure 9.11. Instantaneous eddy viscosity contours for the SST-MS hybrid turbulence model. All three hybrid turbulence models have similar trends on the three grid levels. Significant turbulent structure can be seen in the fine grid solutions with both large-scale and small-scale turbulent structure present. The mid level solutions have large-scale structure present, but the smaller turbulent scales are absent. The mid level grids also have structure that appears to be more periodic than does the fine grid solution. The coarse level solutions show almost no turbulent structure and produce an almost steady state wake away from the cylinder. The level of eddy viscosity in the wake differs for the three turbulent models. The eddy viscosity for the SA-DES model shown in Fig. 9.7 indicates that this model effectively shuts off the eddy viscosity outside the boundary layer. The eddy viscosity for the SST-DES and the SST-MS models is reduced as the grid is refined as would be expected for these models. The SST-DES model tends to predict higher levels of eddy viscosity in the wake region than does the SST-MS model. The SST-DES model also predicts higher eddy viscosities along the 9-10 edges of the wake than either the SA-DES or SST-MS models. Both the SSTDES and SST-MS models are tending toward a RANS type solution in the far wake of the coarse grid solution. The ratio of the turbulent length scale to the local grid length scale for the SSTMS model is shown in Fig. 9.12 for the three grid levels. a. Fine b. Mid c. Coarse Figure 9.12. Instantaneous ratio of turbulent length scale to the grid length scale for the SST-MS hybrid turbulence model. The length scale ratio seems to scale with the grid refinement for the mid to fine grid results, which indicates that the turbulent length scales predicted on the mid and fine grids are similar. The mid level grid results in a length scale ratio of greater than two for most of the wake region. This is the desired performance trend with grid refinement for hybrid turbulence models. The coarse grid results do not seem to continue the scaling trend. The power spectral densities (PSD) of the axial normal force coefficient are shown in Figs. 9.13-9.15 for the three hybrid turbulence models. The mid and fine level grid results are in general agreement for the SA-DES hybrid model. The location and intensity of the primary spectral peak is still varying with the SST-DES and SST-MS hybrid models. This may be due to the fact that the SST-DES and SST-MS hybrid models use a ratio of the local turbulent length scale to the local grid length scale to determine the subgrid eddy viscosity, while the SA-DES model only uses the local grid length scale. The coarse grid solutions are seen to be significantly different and the energy is contained in narrower peaks for all of the hybrid models. 9-11 Figure 9.13. Power spectral density of the normal force coefficient for the SADES hybrid model. Figure 9.14. Power spectral density of the normal force coefficient for the SSTDES hybrid model. 9-12 Figure 9.15. Power spectral density of the normal force coefficient for the SSTMS hybrid model. The average integrated drag coefficient (Cd), standard deviation of the lift coefficient (σ(Cl)), and the peak lift coefficient Strouhal number (St) for the 3 hybrid models computed on the three grid levels are shown in Table 9.1. The mid and fine level grid results are in general agreement with each other and the experimental data9,10,11. The results are also in reasonable agreement with the Ref. 12 results obtained at a Reynolds number of 3.0x106. Model SA SST SA-DES SST-DES SST-MS SA SST SA-DES SST-DES SST-MS SA-DES SST-DES SST-MS DATA Ref. 9 DATA Ref. 10 DATA Ref. 11 SA-DES Ref. 12 Grid Coarse Coarse Coarse Coarse Coarse Mid Mid Mid Mid Mid Fine Fine Fine Average Cd 0.266 0.192 0.298 0.216 0.196 0.511 0.517 0.626 0.589 0.556 0.518 0.590 0.585 0.505-0.540 0.79 0.41-0.51 σ(Cl) 0.163 0.103 0.253 0.171 0.115 0.428 0.520 0.186 0.274 0.135 0.182 0.212 0.238 0.0575-0.0770 0.06-0.13 St 0.232 0.244 0.232 0.250 0.250 0.281 0.281 0.256 0.293 0.281 0.281 0.256 0.268 0.306-0.308 0.27 0.29 0.33-0.35 Table 9.1 Force coefficient and Strouhal number predictions on the 3D circular cylinder. 9-13 9.3 WICS Bay Unsteady Navier-Stokes calculations were also performed for the WICS (Weapons Internal Carriage and Separation) L/D=4.5 weapons bay13 for M=0.95 and Re=2.5x106/ft. The wind tunnel data were obtained in the AEDC four-foot transonic wind tunnel. The weapons bay was 18 in. long, 4 in. wide, and 4 in. deep. The computational geometry was a flat plate that extended 15 in. upstream of the bay to match the experiment and 25 in. downstream of the bay. The sides of the computational grid extended 50 in. on either side of the bay centerline. The full bay geometry was modeled using wall functions. The wall spacing was chosen as 0.0075 in., which corresponds to a y+ of 50 on the upstream plate. This wall spacing was used for all the grids in the grid refinement study. Grid resolution effects were investigated by modifying the bay grid and the grids in the vicinity of the bay grid. The grids upstream and to the sides of the bay were not modified. Details of the grid systems and grid spacings used in this study are shown in Table 9.2. The centerline plane of the three bay grids is shown in Figure 9.16. Grid Fine Medium Coarse Total Points 1.8x106 1.1x106 7.9x105 Bay Grid Dimensions 121x61x61 71x41x41 61x31x31 Bay Grid Δxmax 0.3 in. 0.6 in. 0.75 in. Bay Grid Δymax 0.1 in. 0.2 in. 0.3 in. Bay Grid Δzmax 0.1 in. 0.2 in. 0.3 in. Table 9.2. Parameters for the grid refinement study. A time step study was performed using the SST-MS hybrid model and the medium grid set. Time steps of 8.0x10-6 seconds to 8.0x10-5 seconds were evaluated. A time step of 8.0x10-6 seconds corresponds to 260 time steps per cycle for the primary spectral mode of the cavity at 480 hz. The calculations were run 2000 steps to remove the initial transients and all unsteady results were processed over the last 8192 time steps. The results are compared to data (Ref. 12) for the K16 and K18 transducer locations in the bay in Fig. 9.17. K16 is located on the bay ceiling centerline 0.275 inches from the back wall. K18 is located on the bay back wall centerline 0.725 inches from the bay opening. The sound pressure levels predicted by the hybrid turbulence models includes only the contribution of the grid resolved turbulent scales. Subgrid turbulent energy, which includes most of the energy in the boundary layer, is not included. Thus the spectrum predicted using the hybrid models should have lower energy at higher frequencies (i.e. should have a more rapid roll off) than the experimental data. 9-14 Figure 9.16. Centerline plane of the bay grids. Figure 9.17a. Sound pressure level spectra at the K16 location for varying time step. 9-15 Figure 9.17b. Sound pressure level spectra at the K18 location for varying time step. The agreement with the measured spectra is seen to improve at the higher frequencies as the time step is reduced. The solution with the largest time step predicts the primary spectral peaks at too low a frequency. The solution with the smallest time step provides reasonable agreement for the first four spectral peaks. The solution using 1.6x10-5 second time step has reasonable agreement with the first three spectral peaks. All subsequent solutions were performed using the 1.6x10-5 second time step. A comparison of the sound pressure level for the three hybrid models, the Spalart-Allmaras model, and the SST model for the medium grid are shown in Fig. 9.18. The hybrid models show similar trends throughout the spectral range. The two RANS models predict the first two spectral peaks, but are well below the data away from the peaks. This is an indication that the RANS models are providing too much damping of the unsteady solutions. 9-16 Figure 9.18a. Sound pressure level spectra at the K16 location for the medium grid. Figure 18b. Sound pressure level spectra at the K18 location for the medium grid. Instantaneous Mach number and eddy viscosity contours on the bay centerline are shown in Figs. 9.19-9.21. 9-17 a. Fine b. Medium c. Coarse Figure 9.19a. Instantaneous Mach number contours on the bay centerline using the SA-DES hybrid model. a. Fine b. Medium c. Coarse Figure 9.19b. Instantaneous eddy viscosity contours on the bay centerline using the SA-DES hybrid model. a. Fine b. Medium c. Coarse Figure 9.20a. Instantaneous Mach number contours on the bay centerline using the SST-DES hybrid model. 9-18 a. Fine b. Medium c. Coarse Figure 9.20b. Instantaneous eddy viscosity contours on the bay centerline using the SST-DES hybrid model. a. Fine b. Medium c. Coarse Figure 9.21a. Instantaneous Mach number contours on the bay centerline using the SST-MS hybrid model. a. Fine b. Medium c. Coarse Figure 9.21b. Instantaneous eddy viscosity contours on the bay centerline using the SST-MS hybrid model. The turbulent structure shown in the Mach number contours is seen to dramatically increase with grid resolution for all three hybrid models. The shear layer is seen to grow and then burst as it traverses the bay. The eddy viscosity decreases with increasing grid resolution as is expected for all three hybrid models. The SST-DES hybrid model has the largest levels of eddy viscosity, and the maximum eddy viscosity is located in the shear layer. The ratio of the 9-19 turbulence length scale to the grid length scale (Lt/Lg) for the SST-MS model is shown in Fig. 9.22. The contours scale with the grid refinement as is expected for this model. The time averaged pressure coefficient on the bay ceiling is shown in Fig. 9.23 for the three hybrid models. The sound pressure level on the bay ceiling is shown in Fig. 9.24. The solutions are similar for the three models for both of these properties. There is some sensitivity to grid density for the average pressure coefficient on the downstream end of the ceiling. The medium and fine grids produce similar results for the sound pressure. The coarse grid sound pressure is only slightly different from the finer grids. This indicates that average pressure coefficient and sound pressure level are not overly sensitive to grid density. Fine Medium Coarse Figure 9.22. Instantaneous ratio of turbulent length scale to the grid length scale contours on the bay centerline using the SST-MS hybrid model. 9-20 Figure 9.23a. Average pressure coefficient on the WICS ceiling using the SADES model. Figure 9.23b. Average pressure coefficient on the WICS ceiling using the SSTDES model. 9-21 Figure 9.23c. Average pressure coefficient on the WICS ceiling using the SSTMS model. Figure 9.24a. Sound pressure level on the WICS ceiling using the SA-DES model. 9-22 Figure 9.24b. Sound pressure level on the WICS ceiling using the SST-DES model. Figure 9.24c. Sound pressure level on the WICS ceiling using the SST-MS model. The sound pressure level spectra at the K16 and K18 locations for the three hybrid models are shown in Figs. 9.25-9.27. The medium and fine grid spectra are in reasonable agreement with the first three spectral peaks for all three hybrid 9-23 models. The coarse grid spectra are not quite as good, but are not totally unacceptable. The fine grid solutions underpredict the primary spectral peak, while the coarse and medium grids generally overpredict the primary spectral peak. Spectral results seem to be more sensitive to grid density than the time averaged quantities. Figure 9.25a. Sound pressure level spectrum at the K16 location using the SADES model. Figure 9.25b. Sound pressure level spectrum at the K18 location using the SADES model. 9-24 Figure 9.26a. Sound pressure level spectrum at the K16 location using the SSTDES model. Figure 9.26b. Sound pressure level spectrum at the K18 location using the SSTDES model. 9-25 Figure 9.27a. Sound pressure level spectrum at the K16 location using the SSTMS model. Figure 9.27b. Sound pressure level spectrum at the K18 location using the SSTMS model. 9.4 Delayed Detached Eddy Simulation (DDES) Hybrid RANS/LES models were developed to provide a RANS solution in the boundary layer along a body and to transition to LES away from the body. The 9-26 functions used to transition from RANS to LES in these hybrid models can transition too early and lead to solutions that are neither LES nor RANS. This usually occurs in small pockets of separated flow such as a shock induced separation on an airfoil. If the grid is fine enough and the turbulent intensities are high enough in this region the hybrid model will begin to reduce the eddy viscosity level. The reduced level may not be low enough to allow the flow to become unsteady. Recently a correction has been proposed for the DES models to force the models to remain turbulent in the boundary layer. The new formulation is called Delayed Detached Eddy Simulation (DDES). The DDES14 version of the SA model modifies Eq. (9.1) to force the modified wall distance parameter to equal the distance from the wall further into the boundary layer. This delays the transition from RANS to LES. The new modified distance function is given by ~ d = d − f d max(0, d − C DES Δ ) (9.12) where ( f d ≡ 1 − tanh [8rd ]3 ) (9.13) and rd ≡ ν t +ν Sκ 2 d 2 (9.14) fd is designed to be 1 in the LES region where rd << 1, and 0 elsewhere. rd can be thought of as a ratio of the turbulent length scale defined by ν t / S and the wall distance d. This length scale ratio will work for simple turbulent flows, but may have difficulty in more complicated flows since the eddy viscosity is a transported quantity while the strain is a locally derived quantity. It is not difficult to imagine flows where a large value of eddy viscosity is transported into a region of relatively small strain. This could cause the DDES model to transition from LES to RANS mode away from the body. A DDES version of the SST model15 can be derived by replacing Lg in Eq. (9.3) with a wall corrected grid length scale Lgcor defined as L gcor = Lt − f D (Lt − L g ) where (9.15) 9-27 3 ⎡⎛ L ⎞ ⎤ f D = 1 − tanh ⎢⎜1.5 t ⎟ ⎥ d ⎠ ⎥ ⎢⎣⎝ ⎦ (9.16) Here d is the distance to the wall. The turbulent dissipation in the turbulent kinetic energy equation is effectively increased when the grid size length scale Lg is less than the turbulent length scale Lt. The wall correction increases the size of Lgcor near the wall and causes the turbulence model to remain in RANS mode in the near wall region without regard to the local grid size in the boundary layer. Note that the turbulent length scale used here is based on the ratio of the transported turbulent quantities k and ω only. The pressure distribution for the transonic flow over a NACA 001216 wing test can be used to show the benefit of the DDES correction. The free-stream Mach number was 0.799, the angle-of-attack was 2.26 degrees, and the Reynolds number based on chord was 9x106. A 381x81 two dimensional grid was used for this study with 301 points on the airfoil. The pressure coefficient on the airfoil is shown in Fig. 9.28. All of the solutions were run in an unsteady mode, but all of the models produced steady state solutions for this case. The SA and SST RANS model produce similar results for this case. The SA DES model moves the shock on the upper surface upstream and produces a weak shock on the lower surface. This is due to a decrease in eddy viscosity in the boundary layer. The SA DDES model produces a solution that is closer to the SA model result. The SST DES model moves the shock well upstream on the upper surface and produces a shock on the lower surface. The SST DDES model produces a solution similar to the SA DDES model. 9-28 Figure 9.28a Pressure coefficient for a NACA 0012 airfoil for SA models. Figure 9.28b Pressure coefficient for a NACA 0012 airfoil for SST models. ~ The SA DES and SA DDES model modified wall distance ( d ) is shown in Fig. 9.29 for the separated flow region behind the shock at x/c=0.75 on the wing upper surface. The boundary layer edge is approximately at z/c=0.045. Note ~ that d = d for the SA model. The SA DDES model is seen to delay the transition to LES to a much larger distance from the wall than does the SA DES model. The effective turbulent length scale used for transitioning the SST DES and SST DDES models at the same location is shown in Fig. 9.30. The SST DES is seen to transition very near the wall. The SST DDES modified turbulent length scale moves the RANS to LES transition much further away from the wall. 9-29 Figure 9.29 Wall distance parameter for the SA, SA DES, and SA DDES turbulence models. L Figure 9.30 Turbulent length scale for the SST, SST DES and SST DDES turbulence models. 9-30 9.5 Summary These simulations indicate a fundamental difficulty in verification and validation for unsteady flows. As the grid is refined, smaller scale turbulent structures are resolved in the solution. This process will continue until the grid is refined below the Komologrov scale and all of the turbulent scales are grid resolved. The limiting grid refinement case cannot be approached with current computing hardware for most real world high Reynolds problems of interest. Comparison of statistical quantities derived from the unsteady solutions can be made to assure that the important features of the unsteady flow have been captured with a given grid, but the nature of the statistical analysis makes it difficult to assess whether a true “grid convergence” has been achieved. The time step study indicated that about 200 time steps per primary shedding cycle were required for temporal accuracy with the SST-MS hybrid model. Using a larger time step causes the primary spectral peak to occur at a lower frequency than predicted with the smaller time step and than is seen in the data. The simulations on the mid and fine grid systems are in good general agreement indicating that a level of grid convergence can be achieved for the large turbulent scales. The turbulent length scale to grid length scale ratio was greater than two in the wake region of the cylinder and in the bay. This may serve as a rule of thumb for grid resolution for hybrid model applications. Although the primary shedding frequency for the cylinder is approximately right for the coarse grid, the integrated forces and the spectral shape are quite different from the mid and fine grid solutions. This indicates that these turbulence models can produce solutions that appear to capture the relevant physics, but not capture the physical details of the flow. The coarse grid solutions were totally non-physical. The hybrid RANS/LES turbulence models are relatively new and will need to be exercised for a wide variety of problems to determine their accuracy before they become an accepted tool for fluid modelers. They seem to offer much for unsteady flow applications, but issues such as grid sensitivity need to be further addressed. Hopefully more effort will go into these models in the near future so that they can be matured for use in everyday applications. 9-31 Hybrid RANS/LES Application Hints 1. These models are very new and have not been used extensively. They should be used with care and only when validation calculations can be performed. 2. These models are for unsteady applications, and should not be used with local time stepping or with other non-time accurate algorithms. 3. Turbulent flows are three dimensional, and hence these models should be used only in 3D (see Ref. 17). 4. Because of the unsteady nature of these models, they may require a large number of time-steps to obtain a statistically stationary solution for analysis. 5. These models may be sensitive to the computational mesh because the filter function is inversely proportional to the grid spacing. A rule-of-thumb is that the ratio of the turbulent scale to grid length scale should be greater than two in the region of interest when using hybrid models. 6. As with all unsteady applications, care should be taken to be sure the time step is small enough to temporally capture the unsteady phenomena of interest. 7. The models can produce solutions that are neither RANS or LES when applied to flows with small pockets of flow separation such as shock induced separation on an airfoil. Chapter 9 References: 1. 2. 3. 4. 5. 6. 7. Spalart, P., Jou, W-H. Strelets, M., and Allmaras, S. "Comments on the Feasibility of LES for Wings and on a Hybrid RANS/LES Approach," First AFOSR Conference on DNS/LES, August 1997, C. Liu and Z. Liu editors, Greyden Press, Columbus, Ohio. Constatinescu, G. and Squires, K., “LES and DES Investigations of the Turbulent Flow over a Sphere,” AIAA-2000-0540, Jan. 2000. Mitchell, A., Morton, S., and Forsythe, J., “Analysis of Delta Wing Vortical Substructures using Detached-Eddy Simulation,” AIAA-2002-2968, Jun. 2002. Forsythe, J., Squires, K., Wurtzler, K., and Spalart, P, “Detached-Eddy Simulation of Fighter Aircraft at High Alpha,” AIAA-2002-0591. Strelets, M., "Detached Eddy Simulation of Massively Separated Flows," AIAA-2001-0879, Jan. 2001. Nichols, R. and Nelson, C., “Applications of Hybrid RANS/LES Turbulence Models,” AIAA-2003-0083, Jan. 2003. Nelson, C. and Nichols, R., “Evaluation of Hybrid RANS/LES Turbulence Models using an LES Code,” with C. C. Nelson, AIAA-2003-3552, Jun. 2003. 9-32 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. Nichols, R., “Comparison of Hybrid RANS/LES Turbulence Models for a Circular Cylinder and a Cavity,” AIAA Journal, Vol. 44, No. 6, June 2006, pp. 1207-1219. Jones, G. W., Cincotta, J. J., and Walker, R. W., ”Aerodynamic Forces on a Stationary and Oscillating Circular Cylinder at High Reynolds Number,” NASA-TR-R-300, October 1968. Roshko, A., “Experiments on the Flow Past a Circular Cylinder at Very High Reynolds Number,” Journal of Fluid Mechanics, Vol. 10, Part 3, May 1961, pp. 345-356. Schlichting, H., Boundary Layer Theory, 7th Edition, McGraw-Hill, New York, 1979. Travin, A., Shur, M., Strelets, M. and Spalart, P., “Detached-Eddy Simulations Past a Circular Cylinder,” Flow, Turbulence and Combustion, Vol. 63, pp 293-313, 1999. Dix, R. E. and Bauer, R. C., “Experimental and Theoretical Study of Cavity Acoustics,” AEDC-TR-99-4, May 2000. Spalart, P. R., Deck, S., Shur, M. L., Squires, K. D., Strelets, M. Kh., and Travin, A., “A New Version of Detached-Eddy Simulation, Resistant to Ambiguous Grid Densities,” Theor. Comput. Fluid Dyn., Vol. 20, 2006, pp. 181-195. Nichols, R. H., “A Comparison of Hybrid RANS/LES Turbulence Models for a Generic Weapons Bay With and Without a Spoiler,” AIAA-2008-6629, Aug. 2008. Maksymiuk, C. M. and Pulliam, T. H., “Viscous Transonic Airfoil Workshop Results Using ARC2D,” AIAA Paper 87-0415, Jan. 1987. Shur, M., Spalart, P., Squires, K., Strelets, M., and Travin, A., “Three Dimensionality in Reynolds-Averaged Navier-Stokes Solutions around TwoDimensional Geometries,” AIAA Journal, Vol. 43, No. 6, June 2005, pp. 1230-1242. 10-1 10.0 Wall Function Boundary Conditions Turbulent transport models that are applicable all the way to the wall are called low Reynolds number turbulence models since they are valid at low turbulent Reynolds numbers (ρk2/(με)). Low Reynolds number models require tight grid spacing near the wall to resolve the large gradients in velocity, temperature, and turbulent quantities near the wall. This can lead to large grid requirements for 3D Navier-Stokes applications. The stability of a numerical algorithm is limited by the smallest cell size in a grid, so the small cells near the wall can severely restrict the maximum allowable time step for a problem. The near-wall turbulence damping terms are generally expressed in exponential functions. This can lead to difficulties in converging the numerical solution in the near wall region. Turbulence models that do not include wall correction terms in the differential equations are called high Reynolds number turbulence models. High Reynolds number turbulence models rely on empirically derived algebraic models of the near-wall region of the boundary layer to provide boundary condition information to the mean flow Navier-Stokes equations at the first point off the wall. These empirically derived relationships are called wall functions. Since the high gradient region near the wall is modeled with these empirical relationships, the first point off the all may be placed much farther away than with low Reynolds number models. This reduces the number of points required to discretize a flow field and increases the maximum allowable time step. Wall functions have been used by turbulence modelers for a number of years. Wall functions were a matter of necessity during the years before computers began to mature, and are still a convenient and efficient means of modeling many turbulent flows. Most wall function implementations are based on empirical relations for incompressible adiabatic flows. Some implementations include a decoupled relation for the temperature distribution in the lower portion of the boundary layer that allows the wall functions to be applied to flows with heat transfer at the wall. The wall function boundary condition formulation developed here is based on a coupled velocity and temperature boundary layer profile that is applicable to incompressible flows, compressible flows, and flows with heat transfer. The wall functions described here are compatible with most transport type turbulence models. 10.1 Theory There are six fundamental assumptions used in the development of wall function boundary conditions for compressible flows: 1. Analytical expressions are available for the velocity and temperature profiles in the lower part of the boundary layer. 10-2 2. Analytical expressions are available for the turbulent transport variables at the first point off the wall. The pressure is constant in the lower part of the boundary layer. ∂u The shear stress ( τ = τ w = (μ + μ t ) ) is constant in the lower part of ∂y the boundary layer. ∂T The heat transfer ( q = qw = (k + k t ) ) is constant in the lower part of ∂y the boundary layer. There are no chemical reactions (i.e. the chemistry is frozen) in the lower part of the boundary layer. 3. 4. 5. 6. Several empirical relationships for the velocity in a boundary layer are available for use in wall function boundary conditions. A typical boundary layer velocity profile is shown in Fig. 3.2. The simplest profile form is the adiabatic incompressible universal law-of-the-wall u+ = 1 κ log( y + ) + B (10.1) where u+ and y+ are defined in Eqs. 3.10 and 3.11. The universal law-of-thewall is valid for the y+ range of about 30 4) in Ref. 9. In all cases, the model performed satisfactorily. This model was later implemented in GASP4 and again was found to be quite accurate. 3) Wilcox k-ω turbulence model11,12: 11-4 The low Re k-ω model developed by Wilcox was used to predict the transition region. The prediction of the transition region was obtained by tripping the boundary layer at a given point by decreasing the value of the dissipation so as to destabilize the boundary layer and cause transition. The application of this model in Ref. 4 showed that it was not very easy to trip the boundary layer at the desired location due to the sensitivity to the initial conditions. This model predicted a short transition length and over-predicted the peak heat transfer for some cases. 4) Schmidt and Patankar production term modifications13: Schmidt and Patankar have developed modifications to the production term in the turbulent kinetic energy (TKE) equation of the Lam and Bremhorst k-ε model14. These modifications limited the production of the kinetic energy. For the use of this model, a trial and error method was needed to make transition occur at the desired position by varying the inlet conditions. The results did not compare well to the experimental cases in Ref. 4, and the method was found to be very sensitive to the grid spacing near the wall. Due to the defects in this model a few modifications were suggested in Ref. 4. Since it was found that the model was difficult to trigger turbulence transition, a spot with high TKE was introduced into the boundary layer. This spot then grew and caused transition to take place. In order to improve the prediction of the length of the transition region, an exponential function was used for the maximum allowable production of TKE. These modifications improved the results for some cases but gave worse results for other cases. 5) Algebraic transition model15,16: The algebraic transition region model was developed at ONERA/CERT and is described in Arnal15,16. The form of the model of Singer et al.17,18 was implemented in Ref. 4. This model predicts transition by multiplying the eddy viscosity by a transition function before adding it to the fluid viscosity. This function was found to be related to the momentum thickness growth. As a result, in test cases with severe adverse pressure gradients, where the momentum thickness decreases, the model did not produce transition. Theoretically this model should be compatible with any turbulence model. However, it was found that this model did not perform well with two equation models. In Ref. 4 the model was used with the Baldwin-Lomax model. Corrections to the calibration of the transition function for high speed flows were also suggested in Ref. 4. The new model predicted the cases tested better than the original model. 6) Linear combination transition model10: This model was developed by Dey and Narasimha19 and is based on the concept that the transition flow is a combination of the laminar and turbulent flow fields. The contribution from laminar and turbulent values is proportioned based on the intermittency factor developed by Dhawan and Narasimha10 mentioned above. This model requires that a complete laminar flow simulation be run first. This simulation is followed by a turbulent one, with the turbulent boundary layer 11-5 starting at the point of transition. The model then uses these two solutions to generate the transitional solution. For example, the mean velocity (U) and skin friction coefficient (Cf) are given by U = (1 − Γ )U L + ΓU T (11.3) C f = (1 − Γ )C fL + ΓC fT (11.4) In the above equations the subscripts L and T stand for values in the laminar and turbulent boundary layers, respectively. The peak heat transfer was not predicted in the test cases simulated in Ref. 4. The main difficulty in getting accurate results with this model was that one of the modeling constants needs to be modified from case to case to obtain good results. Many researchers such as Abid20 have used the intermittency function from the linear combination model as an algebraic transition region function to proportion the amount of the eddy viscosity added to the fluid viscosity. The results using this method were found to be very similar to the linear combination model mentioned above, but there are some noticeable differences4. The transition length was always under-predicted. For the cases with no pressure gradient and adverse pressure gradients, the heat transfer predicted at the end of transition and through the turbulent region was significantly high. In addition to the above models, efforts have been conducted to modify existing turbulence models for turbulence transition. In Ref. 21, the performance of the Spalart Allmaras (S-A)22 and the Baldwin-Barth (B-B)23 one-equation models and three two-equation models for simulating hypersonic transition were evaluated. The two equation models assessed in Ref. 21 included a low Re k-ε model with the modifications of Nagano and Hishida24, the hybrid k-ω model of Menter (SST)25, and the Wilcox k-ω model26. The Sandia Advanced Code for Compressible Aerothermodynamics Research and Analysis (SACCARA) was used to evaluate these models in Ref. 21 using two flow cases. The first case was the flow over a flat plate at Mach 8 with flow conditions corresponding to an altitude of 15 km, where a perfect gas was assumed. The second flow case considered was the flow over a re-entry flight vehicle at Mach 20 and an altitude of 24.4 km, where real gas effects need to be taken into account. The method employed in Ref. 21 to specify transition from laminar to turbulent flow was as follows. The turbulence transport equations were solved over the entire domain, with a transition plane specified by the user. Upstream of this plane, the effective viscosity was simply the laminar value whereas at downstream the effective viscosity was the sum of the laminar and turbulent viscosities. An advantage of this approach was that the turbulence transport equations were solved over the whole domain, thus promoting turbulent behavior downstream of the transition plane. On the other hand, if the turbulence source terms were simply turned on after the transition plane, the turbulence model might not transition to turbulent flow until farther downstream, depending on the free stream turbulence level. However, a disadvantage of the approach was that a discontinuity in the total 11-6 viscosity (laminar plus turbulent) could occur at the transition plane. All models except the S-A model and the low Re k-ε model predicted the transition at the correct location for the flat plate case at Mach 8. For these two models, the free stream turbulence values needed to be increased. All the models tested in Ref. 21 provided the correct skin friction levels for this case. For the reentry flight vehicle, the wall heat flux predicted by the S-A model, the SST model, and the Wilcox k-ω model were found to be in reasonable agreement with the experimental data. The B-B and the Low Re k-ε model greatly over-predicted the heat flux in the turbulent region. 11.3 Transition Models with Onset Prediction Capability These models not only simulate the characteristic of the transition region, but also predict the onset of transition. 1) k-ζ turbulence/transition model: The k-ζ turbulence model27 was used to study the effect of high disturbance environments (HIDE) on the transitional simulations carried out in conventional hypersonic facilities. Since HIDE cannot be described by linear stability theory, a minimum heat flux criterion was used to determine onset of transition. This is done by assuming initial transition onset points and employing linear interpolation for interior points. After running a few iterations, the minimum heat flux criterion is employed to find the locations where the wall heat flux is a minimum. The solution is independent of the initial guess as long as the initial transition points are ahead of the actual locations. This approach is similar to the WHH model mentioned above. In this case, the eddy viscosity was modified using the formula μ e = (1 − Γ )μ nt + Γμ t (11.5) μ nt is the eddy viscosity due to the non-turbulent fluctuations, and can be calculated as μ nt = 0.09 ρkτ nt (11.6) where τnt is the non-turbulent time scale. In addition, the dissipation time scale in the turbulent kinetic energy equation was also chosen as the combination of time scales of turbulent and non-turbulent fluctuations. The time scale for calculating the eddy viscosity and the dissipation time scale were derived for three different transition mechanisms: cross flow instabilities, second mode instabilities, and HIDE. These three mechanisms were selected because they were believed to be responsible for transition over 3D bodies in conventional hypersonic wind tunnels28. The simulations were performed on an elliptic cone at the Mach number of 7.93 and were compared with experimental results. It was concluded from the results that HIDE had a higher impact on the transition mechanism than the other two mechanisms. The main disadvantage of this model is that it does 11-7 not solve the non-turbulent fluctuations using transport equations, which limits the flexibility of this method. Also, this model requires an initial guess for the transition location. 2) Papp and Dash model29: Papp and Dash proposed a concept analogous to that of the WHH model was used. The SSGZ-J k-ε model developed by So et al.30 was implemented with compressibility corrections for hypersonic flows. An additional transport equation was solved for the non-turbulent fluctuations. The non-turbulent fluctuations included the first- and second-mode mechanisms. The location of the onset of transition31 was said to be the minimum distance along the surface for 1 ν lt Re t ≡ ≥ 1 , where Cμ is the turbulent viscosity coefficient (0.09), ν is the fluid Cμ ν kinematic viscosity, and νlt is the eddy viscosity due to the non-turbulent fluctuations which can be calculated as ν lt = C μ klτ nt (11.6) where τnt is the viscosity time scale obtained for different transition mechanisms, and kl is the laminar turbulent kinetic energy, which is obtained from a transport equation in the Papp and Dash model29. This model was incorporated into a Reynolds-Averaged Navier-Stokes (RANS) flow solver by multiplying the turbulent eddy viscosity by the intermittency before adding to the fluid viscosity. In all cases simulated, the transition onset was properly obtained. However, in some cases the peak in heat transfer was not reproduced correctly. This has been attributed to the algebraic nature of the intermittency function used. This is the biggest disadvantage of the model. 3) Suzen and Huang model [32]: This model uses a transport equation for the intermittency factor. This equation not only reproduces the intermittency distribution of Dhawan and Narasimha10, but also gives a realistic variation of the intermittency in the cross-stream direction. The intermittency transport equation includes source terms from two different models: the Steelant and Dick model33 and the Cho and Chung model34. The model is incorporated into the Navier-Stokes solvers by simply multiplying the eddy viscosity obtained from the turbulence calculations with the intermittency factor. The Menter’s shear stress transport (SST25) model was used to calculate the turbulent quantities. The onset of transition was determined by comparing the local Reynolds number with a transition onset Reynolds number (Reθ) calculated using the correlation of Huang and Xiong35, where Reθ is a function of the free stream turbulent intensity and the acceleration parameter. This model was tested for zero- and variable-pressure gradient flows with different free stream turbulence intensities. The numerical result showed good agreement with the experimental data of Savill36,37. This model is not a single 11-8 point model since it uses the free stream turbulence intensity value to calculate the transition onset Reynolds number, which requires global parameters. 4) Walters and Leylek model38: This model is based on the concept that bypass transition is caused by very high amplitude stream wise fluctuations. These fluctuations are very different from turbulent fluctuations. Mayle and Schulz39 proposed a second kinetic energy equation to describe these fluctuations. This kinetic energy was called laminar kinetic energy kL. In the near-wall region, the turbulent kinetic energy (TKE, kT) was split into small-scale energy and large scale energy. The small-scale energy (kT,s) contributes directly to the turbulence production, and the large-scale energy (kT,l) contributes to the production of laminar kinetic energy. These two energies can be calculated from the kT based on the turbulent length scale. The eddy viscosities based on both scales are calculated from the respective-scale kinetic energies. For the onset of transition, a parameter is calculated from kT, the kinematic viscosity and the wall distance. When this parameter exceeds a certain threshold, transition is assumed to start. The onset of transition is associated with the reduction of kL and the consequent increase of kT (indicating the breakdown of laminar fluctuations into turbulence). This model was incorporated into a RANS flow solver for the calculation of the total eddy viscosity and eddy thermal diffusivity to account for contributions from the small-scale as well as large-scale turbulent kinetic energies. For all test cases simulated, the model responded correctly to increases in the free stream intensity. It yielded reasonable results for cases with high pressure gradients and streamline curvatures. Advantages of this method are that it is very simple to implement it into the existing CFD codes since it is based on a RANS framework. This is a single point transition model meaning that it requires only local information, which makes this method easily applicable to unstructured and parallel computations. The low-Re k-ε models are typically not calibrated for transition prediction, but provide the transition location as a by-product of their viscous sublayer formulation. Since this transition model is developed based on the low-Re k-ε model, the embedded viscous sublayer formulation coupled with the added transition prediction capability cannot be calibrated independently. Hence, a change in the transition formulation would affect the solution in the fully turbulent region. In addition, it is generally observed that these models are not flexible enough to sufficiently cover the wide range of transition mechanisms observed in reality40,41. 5) Local Correlation Based Transition Model (γ − Reθ model or LCTM) 40,41: This model is based on the vorticity Reynolds number. The vorticity Reynolds number is an extremely important parameter and is a local property and can be easily calculated in CFD codes. The maximum value of the vorticity Reynolds number in a boundary layer profile is directly proportional to the momentum thickness Reynolds number. The vorticity Reynolds number is used in triggering transition instead of directly using the momentum thickness Reynolds number. This model solves a transport equation for intermittency and also a transport equation for the Reynolds number based on the transition onset momentum 11-9 thickness. The first transport equation includes two terms that control production. These are Flength, a parameter which controls the length of transition zone, and Reθc which is the momentum thickness Reynolds number at the point where the intermittency starts to increase in the boundary layer. These two variables are calculated from empirical functions of the transition momentum thickness Reynolds number ( Reθ t ). A second transport equation is required to solve for Reθ t and to include the non-local influence of the turbulence intensity, which varies with the free stream turbulence kinetic energy and the free stream velocity. In the case of flows with boundary layer separation, this transition model is modified so that the intermittency is allowed to exceed the unity when the boundary layer separates. This event results in larger production of kinetic energy leading to correct prediction of reattachment41. This model is applied by modifying the production and destruction terms of the original SST model using the intermittency. The model was validated with many complicated 2D and 3D configurations. In all cases, good agreement with the experimental data was obtained. This model offers two main advantages: 1) it is based on local variables; 2) it is very flexible and can be used for any mechanism as long as the empirical correlation can be formulated. However, the empirical correlations used with this model are proprietary. 6) The model of Lodefier et al.42: This model is also based on the concept of pre-transitional fluctuations similar to the Walters and Leylek model. However, this model uses the concept of intermittency to describe the transition region. The intermittency equation used was proposed by Steelant and Dick43. The production term of this intermittency equation was modified by multiplying it with a new factor which is used to locate the start of transition. This factor is zero before the start of transition and rapidly goes to unity after the onset point. Similar to LCTM, the vorticity Reynolds number is used in triggering transition instead of directly using the momentum thickness Reynolds number ( Reθ t ). Unlike LCTM, the equation used to calculate the critical value of Reθ t for transition is calculated from the local free stream turbulence intensity and not from a transport equation. The empirical correlation used for Reθ t does not include a pressure gradient term. The model is incorporated into the SST model both by multiplying the eddy viscosity with the intermittency and by modifying the production terms of the k and ω equations. These modifications are used to ensure that the turbulence quantities have small non-zero values at the start of transition as in the concept of pre-transitional fluctuations. The main disadvantage of this model is that it uses the free stream intensity to determine the onset of transition, which makes the model non-local unlike LCTM. 7) Model of Lian and Shyy44: This model was developed for simulation of flow around the wing of a micro air vehicle (MAV). The approach used in this model was to couple an 11-10 incompressible RANS solver with the en method. The k-ω model of Wilcox [11] was selected for modeling turbulence in the RANS solver. This coupling is accomplished as follows. The computation is started with the solution of the RANS equations; however, the eddy viscosity is not added to the effective viscosity. The boundary layer parameters required for the solution of the en method are extracted from the Navier-Stokes solutions to evaluate the amplification factor. Once the threshold value of the n-factor is reached, the flow is allowed to become turbulent by multiplying the eddy viscosity with the intermittency factor and adding it to the effective viscosity. The intermittency in this case is calculated from an empirical formula. The en method employed in this case is based on the assumptions that the initial disturbance is small and that the boundary layer is thin. 8) Model of Arthur and Atkin3: This method is based on linear stability theory (en method) applied within a RANS framework. The overall process is as follows. The viscous flow over the configuration of interest is first calculated with an initial guess of the transition onset location. A series of pressure distributions is extracted from the RANS solution at different “line- of-sight” positions across the span. These pressure distributions are fed into a boundary layer code to predict the boundary layer parameters with great accuracy and fidelity. The stability analysis, together with some “n” factor criterion is conducted to yield the transition location. This information is then passed onto the RANS solver for further solution. This process is continued until the transition location and the pressure distribution are converged. For flows with high pressure gradients, it was found that the predicted transition location can move upstream more easily than downstream during the iteration process. Thus, for this method it is essential that the initial guess is downstream of the final, predicted transition location. The method does not have any intermittency model to predict the nature of the region of transition. 11.4 LCTM Applied to a Flat Plate The LCTM method solves two additional transport equations for the intermittency (γ) and the transition momentum thickness Reynolds number ( Reθ t ). The intermittency equation can be written as ∂ D(ργ ) = Pγ − Eγ + ∂x j Dt ⎡⎛ μ ⎢⎜⎜ μ + t σf ⎢⎣⎝ ⎞ ∂γ ⎤ ⎟ ⎟ ∂x j ⎥⎥ ⎠ ⎦ (11.7) The intermittency equation source term is given by Pγ 1 = Flength ca1 ρS [γFonset ]0.5 (1 − ce1γ ) (11.8) 11-11 where S is the strain rate. This term is designed to be zero in the laminar boundary layer and active everywhere the local vorticity Reynolds number exceeds the local transition onset criteria. The vorticity Reynolds number is defined as Reν = ρy 2 Ω μ (11.9) where Ω is the vorticity magnitude. The vorticity Reynolds number is assumed to be proportional to the momentum thickness Reynolds number Reθ = max(Reν ) 2.193 (11.10) Fonset is defined as follows Fonset 1 = Reν 2.193 Reθ c (11.11) ( ( )) 4 Fonset 2 = min max Fonset1 , Fonset 1 ,2 (11.12) ⎛ ⎛ R ⎞3 ⎞ = max⎜1 − ⎜ T ⎟ ,0 ⎟ ⎜ ⎝ 2.5 ⎠ ⎟ ⎝ ⎠ (11.13) Fonset 3 ( Fonset = max Fonset 2 − Fonset 3 ,0 ) (11.14) where RT is defined as (ρk)/(μω). The two functions, Flength and Reθc must still be defined. Reθc is the critical Reynolds number where the intermittency first starts to increase in the boundary layer. This occurs upstream of the transition Reynolds number, Reθt, because the turbulence must grow to a large enough level to trigger transition before any change in the laminar profile can be seen. Hence Reθc is the location where turbulence starts to grow and Reθt is the location where the velocity profile first departs from a laminar profile. Flength is an empirical correlation that controls the length of the transition region. Both Reθc and Flength are defined as functions of a second transport equation variable Reθ t . The destruction/relaminarzation source term is defined as Eγ = ca 2 ρΩγFturb (ce 2γ − 1) where (11.15) 11-12 Fturb = e ⎛R ⎞ −⎜ T ⎟ ⎝ 4 ⎠ 4 (11.16) The constants for this equation are ce1=1.0, ce2=50.0, ca1=2.0, ca2=0.06, and σf=1.0. The transport equation for Reθ t is ( ) D ρ Reθt ∂ Reθt ⎤ ∂ ⎡ = Pθt + ⎢σ θt (μ + μ t ) ⎥ Dt ∂x j ⎣⎢ ∂x j ⎦⎥ (11.17) Outside the boundary layer, the source term Pθt is designed to force the transported scalar Reθ t to match the local value of Reθt calculated from an empirical correlation. The source term is defined as Pθt = cθt ρ t (Re θt ) − Reθt (1.0 − Fθt ) (11.18) Here t is a time scale defined as t= 500 μ ρU 2 (11.19) The blending function Fθt is defined as 4 2 ⎛ y⎞ ⎛ ⎞ ⎛ −⎜ ⎟ ⎛ γ − 1 / ce 2 ⎞ ⎞⎟ ⎜ ⎜ δ ⎝ ⎠ ⎟⎟ , 1.0 ⎟⎟ , 1.0 − ⎜⎜ Fθt = min⎜ max Fwake e ⎜ ⎟ ⎜ ⎟ ⎝ 1.0 − 1 / ce 2 ⎠ ⎠ ⎝ ⎝ ⎠ where (11.20) θ BL = Reθt μ ρU (11.21) δ BL = 15 θ BL 2 (11.22) 50Ωy δ BL U (11.23) δ= Fwake = e ⎛ Reω ⎞ −⎜ ⎟ ⎝ 1 x105 ⎠ 2 (11.24) 11-13 Re ω = ρωy 2 μ (11.25) The Fwake function ensures that the wake function is not active in wake regions downstream of an airfoil. The system of equations is closed except for the empirical relationships for Flength, Reθc, and Reθt. The first two functions are considered proprietary by the developers. The relationship for Reθt is given by ⎧⎛ 0.2196 ⎞ ⎟ F (λθ ), if Tu ≤ 1.3 ⎪⎜1173.51 − 589.428Tu + Reθt = ⎨⎝ Tu 2 ⎠ −0.671 ⎪ 331.5(Tu − 0.5658) F (λθ ), if Tu > 1.3 ⎩ (11.26) Tu is the turbulence intensity defined as Tu = 2k x100 U∞ (11.27) where k is the turbulent kinetic energy. F(λθ) is an empirical relationship for the effect of pressure gradient given by ⎛ Tu ⎞ ⎧ −⎜ ⎟ ⎪1 + (12.986λθ + 123.66λθ2 − 405.689λθ3 )e ⎝ 1.5 ⎠ , λθ ≤ 0 F (λθ ) = ⎨ ⎛ −Tu ⎞ ⎜ ⎟ ⎪ −35 λθ 1 + 0.275 1 − e e ⎝ 0.5 ⎠ , λθ >0 ⎩ 1.5 ( ) (11.28) where λθ is defined as λθ = ρθ 2 dU U ds (11.29) Note this is only valid for -0.1<λθ<0.1. For a zero pressure gradient case, λθ=0 and F(λθ)=1. The remaining two correlation functions are assumed to be Reθc = Reθt (− 4.45 x10 −4 Reθt + 0.92) Flength = 8.5 x10 7 Reθ3t (11.30) (11.31) The LCTM is used to predict the skin friction on a low speed flat plate for free stream turbulence intensities of 0.3, 0.9, 3.3, and 6.5 percent in Fig. 11.1. The 11-14 Model does an excellent job of predicting the onset and length of the transition region for all but the highest free stream turbulence level. The Tu=6.5 case may need more grid refinement near the leading edge of the plate to adequately resolve the transition process. Figure 11.1 Skin friction for a flat plate predicted using LCTM. Chapter 11 References: 1. Smith, A. M. O., and Gamberoni, N., “Transition, Pressure Gradient and Stability Theory,” Douglas Aircraft Co. Rept. ES26388, El Segundo, California, 1956. 2. Van Ingen, J. L., “A Suggested Semi-empirical Method for the Calculation of the Boundary Layer Transition Region,” University of Technology, Department of Aero. Eng., Rep. UTH-74, Delft, 1956. 3. Arthur, M. T., and Atkin, C. J., “Transition Modeling for Viscous Flow Prediction,” AIAA 2006-3052, 2006. 4. McKeel, S. A., “Numerical Simulation of the Transition Region in Hypersonic Flow,” PhD dissertation, Blacksburg, Virginia, 1996. 5. Stainback, P., Fisher, M., and Wagner, R., “Effects of Wind-Tunnel Disturbances on Hypersonic Boundary Transition,” AIAA Paper 72-181, 1972. 6. Holden, M. and Chadwick, K., “Studies of Laminar, Transitional, and Turbulent Hypersonic Flows Over Curved Compression Surfaces,” Tech. Rep. 2610-5, Calspan-UB Research Center, 1994. 7. Kimmel, R., 1993, “Experimental Transition Zone Lengths in Pressure Gradient in Hypersonic Flow,” Symposium on Transitional and Turbulent Compressible Flows (Kral, L. and Zang, T., eds.), Vol. FED 151, pp. 117-127, ASME Fluids Engineering Conference, 1993. 11-15 8. Kimmel, R., “The Effect of Pressure Gradients on Transition Zone Length in Hypersonic Boundary Layers,” Tech. Rep. WL-TR-94-3012, Wright Laboratory, Wright Patterson AFB. 9. Warren, E. S., Harris, J. E. and Hassan, H.A., “Transition Model for HighSpeed flow,” AIAA Journal, Vol. 33, No. 8, 1995, pp. 1391-1397. 10. Dhawan, S., and Narasimha, R., “Some Properties of Boundary Layer Flow During Transition from Laminar to Turbulent Motion,” Journal of Fluid Mechanics, Vol. 3, No. 4, 1958, pp. 418-436. 11. Wilcox, D., “Turbulence Modeling for CFD,” La Canada, CA, DCW Industries, 1992. 12. Wilcox, D., “Turbulence and Transition Modeling for High-Speed Flows,” NASA CR 191473, 1993. 13. Schmidt, R. and Patankar, S., “Simulating Boundary Layer Transition with Low-Reynolds Number k-e Turbulence Models: Part 2 - An Approach to Improving the Predictions,” Journal of Turbomachinery, Vol. 113, 1991, pp. 18-26. 14. Lam, C. and Bremhorst, K., “A Modified Form of the k-e Model for Predicting Wall Turbulence,” ASME Journal of Fluids Engineering, Vol. 103, pp. 456460, 1981. 15. Arnal, D., “Three-Dimensional Boundary Layers: Laminar-Turbulent Transition,” AGARD-FDP-VKI Special Course on Calcul Des Limites Tridimensionelles Avec Ou Sans Decollement, 1986. 16. Arnal, D., “Laminar-Turbulent Transition Problems in Supersonic and Hypersonic Flows,” AGARD-FDP-VKI Special Course on Aerothermodynamics of Hypersonic Vehicles, 1988. 17. Singer, B., Dinavahi, S., and Iyer, V., “Testing of Transition-Region Models: Test Cases and Data,” NASA CR 4371, 1991. 18. Singer, B., “Modeling the Transition Region,” NASA CR 4371, 1993. 19. Dey, J. and Narasimha, R., “An Integral Method for the Calculation of 2D Transitional Boundary Layers,” Journal of Aircraft, Vol.27, no.10, 1990, pp.859-865. 20. Abid, R., “A Study of Turbulence Models for Prediction of Transitional Boundary Layers,” Instability, Transition, and Turbulence (Hussaini, Kumar, and Streett, eds.), Springer-Verlag, 1992. 21. Roy, C. J., and Blottner, F. G., “Assessment of One- and Two- Equation Turbulence Models for Hypersonic Transitional Flows,” Journal of Spacecrafts and Rockets, Vol. 38, No.5, 2001, pp. 313-325. 22. Spalart, P. R., and Allmaras, S. R., “A One-Equation Turbulence Model for Aerodynamic Flows,” AIAA 92-0439, 1992. 23. Baldwin, B. S., and Barth, T. J., “A One-Equation Transport Model for High Reynolds Number Wall-Bounded Flows,” NASA TM-102847, 1990. 24. Nagano, Y., and Hishida, M., “Improved Form of the k-ε Model for Wall Turbulent Shear Flows,” J. of Fluid Eng., 109(2), 1987, pp. 156-160. 25. Menter, F. R., “Two-Equation Eddy Viscosity Turbulence Models for Engineering Applications,” AIAA J., Vol. 32, No. 8, 1994, pp. 1598-1605. 11-16 26. Wilcox, D. C, Turbulence Modeling for CFD, 2nd ed., La Canada, CA, DCW Industries, 1998, Chap.5, pp. 119-122 and 227-271. 27. Xiao, X., Edwards, J. R. and Hassan, H. A., “Transitional flow over an Elliptic Cone at Mach 8,” Journal of Spacecraft and Rockets, Vol. 38, No. 6, 2001. 28. Mokovin, M. V., “Bypass Transition to Turbulence and Research Desiderata,” Transition in Turbines, edited by R. W. Graham, NASA CP2386, 1885, pp. 161-204. 29. Papp, J. L. and Dash, S. M., “Rapid Engineering Approach to Modeling Hypersonic Laminar-to-Turbulent Transitional Flows,” Journal of Spacecraft and Rockets, Vol. 42, No. 3, 2005, pp. 467-475. 30. So, R. M. C., Sarkar, S., Gerodimos, G., and Zhang, J., “A Dissipation Rate Equation for Low-Reynolds-Number and Near-Wall Turbulence,” Theoretical Computational Fluid Dynamics, Vol. 9, 1997, pp. 47-63. 31. Warren, E. W. and Hassan, H. A., “Alternative to the en Method for Determining Onset of Transition,” AIAA Journal, Vol. 36, No. 1, 1998, pp. 111113. 32. Suzen, Y. B. and Huang, P. G., “An Intermittency Transport Equation for Modeling Flow Transition,” AIAA 2000-0287, 2000. 33. Steelant, J. and Dick, E., “Modeling of Bypass Transition with Conditioned Navier-Stokes Equations Coupled to an Intermittency Transport Equation,” International Journal for Numerical Methods in Fluids, Vol. 23, 1996, pp. 193220. 34. Cho, J. R., and Chung, M. K., “A k-e-v Equation Turbulence Model,” Journal of Fluid Mechanics, Vol. 237, 1992, pp. 301-322. 35. Huang, P. G., and Xiong, G., “Transition and Turbulence Modeling of Low Pressure Turbine Flows,” AIAA 98-0339, Reno, NV, 1998. 36. Savill, A. M., “Some Recent Progress in the Turbulence Modeling of By-pass Transition,” Near-Wall Turbulent Flows, edited by R. M. C. So, C. G. Speziale and B. E. Launder, Elsevier Science Publishers B.V., 1993, pp. 829-848. 37. Savill, A. M., “Further Progress in the Turbulence Modeling of By-pass Transition,” Engineering Turbulence Modeling and Experiments 2, edited by W. Rodi and F. Martelli, Elsevier Science Publishers B.V., 1993, pp. 583-592. 38. Walters, D. K. and Leylek, J. H., “A New Model for Boundary Layer Transition Using a Single-Point RANS Approach,” Journal of Turbomachinery, Vol. 126, Issue 1, 2004, pp. 193-202. 39. Mayle, R. E., and Schulz, A., ‘‘The Path to Predicting Bypass Transition,’’ Journal of Turbomachinery, Vol. 119, 1997, pp. 405–411. 40. Langtry, R. B. and Menter, F. R., “Transition Modeling for General CFD Applications in Aeronautics,” AIAA 2005-522, 2005. 41. Menter, F. R., Langtry, R. B., Likki, S. R., Suzen, Y. B., Huang, P. G., Volker, S., “A Correlation-based Transition Model Using Local Variables Part I Model Formulation,” Journal of Turbomachinery, Vol. 128, 2006, pp. 413-422. 42. Lodefier, K., Merci, B., De Langhe, C., and Dick, E., “Intermittency Based RANS Bypass Transition Modeling,” Progress in Computational Fluid Dynamics, Vol. 6, Nos. 1/2/3, 2006. 11-17 43. Steelant, J., and Dick, E., “Modeling of Laminar- Turbulent Transition for High Freestream Turbulence,” Journal of Fluids Engineering, Vol. 123, No. 1, 2001, pp. 22-30. 44. Lian, Y., and Shyy, W., “Laminar-Turbulent Transition of a Low Reynolds Number Rigid or Flexible Airfoil,” AIAA Journal, Vol. 45, No .7, 2007, pp. 1501-1513.
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.7 Linearized : Yes XMP Toolkit : Adobe XMP Core 4.0-c316 44.253921, Sun Oct 01 2006 17:14:39 Format : application/pdf Title : Microsoft Word - TableOfContents.DOC Creator : bnichols Create Date : 2010:09:15 16:34:19-05:00 Creator Tool : PScript5.dll Version 5.2.2 Modify Date : 2010:09:15 16:47:55-05:00 Metadata Date : 2010:09:15 16:47:55-05:00 Producer : Acrobat Distiller 8.1.0 (Windows) Document ID : uuid:7fb61758-a52c-4810-b5ac-33ddb01f5d80 Instance ID : uuid:0b16e88a-f2e4-4e3e-b3fd-990ce54ab315 Page Count : 214 Author : bnicholsEXIF Metadata provided by EXIF.tools