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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content A SUBGRID-SCALE MODEL FOR LARGE-EDDY SIMULATION BASED ON THE PHYSICS OF INTERSCALE ENERGY TRANSFER IN TURBULENCE by Brian Wayne Anderson A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (AEROSPACE ENGINEERING) May 2012 Copyright 2012 Brian Wayne Anderson Dedication To the inspirations in my life: my parents, my sisters, and my loving wife. ii Acknowledgements I am indebted to Peter Diamessis and Tawan Tantikul for their help in learning the spectral solver applied in this research. Above all, I am grateful to my advisor Andrzej Domaradzki who taught me about turbulence. Brian Wayne Anderson iii Table of Contents Dedication ii Acknowledgements iii List Of Tables vi List Of Figures vii Abstract ix Chapter 1: Introduction 1 Chapter 2: Fundamentals of Turbulence Modeling 4 2.1 Background of Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Reynolds Averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Reynolds-Averaged Navier-Stokes Modeling . . . . . . . . . . . . . . . . . 12 2.4 Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5 Large-Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Chapter 3: Subgrid-Scale Modeling in LES 23 3.1 The Smagorinsky Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 The Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 The Scale-Similarity Model . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 Mixed Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Chapter 4: Interscale Energy Transfer in Turbulence 42 4.1 The Energy Cascade in Spectral Space . . . . . . . . . . . . . . . . . . . . 42 4.2 Formulation of Partial Nonlinear Terms . . . . . . . . . . . . . . . . . . . 47 4.3 Interscale Energy Transfer Terms . . . . . . . . . . . . . . . . . . . . . . 55 Chapter 5: Development of the Interscale Transfer Model 61 5.1 Basic Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 Galilean Invariance of the Modeled SGS Stress . . . . . . . . . . . . . . . 64 5.3 Interscale Energy Transfer of the Similarity Model . . . . . . . . . . . . . 68 5.4 Test Filtering in LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.5 Filtering on a Non-uniform Mesh . . . . . . . . . . . . . . . . . . . . . . . 79 iv Chapter 6: Model Testing 81 6.1 Channel Flow Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Numerical Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.3 Test Filter Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.4 Case Descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.5 Modeling Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.5.1 Mean Velocity Results . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.5.2 RMS Velocity Results . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.5.3 Energy Dissipation Results . . . . . . . . . . . . . . . . . . . . . . 118 Chapter 7: Additional Considerations 131 7.1 Non-Symmetry of the Proposed SGS Stress Tensor . . . . . . . . . . . . . 131 7.2 Series Formulation of the Interscale Transfer Model . . . . . . . . . . . . . 142 7.3 SGS Model Computational Expense . . . . . . . . . . . . . . . . . . . . . 153 Chapter 8: Conclusions 158 References 164 v List Of Tables 6.1 LES case descriptions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.2 DNS case descriptions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.3 Integrated energy dissipation. . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.4 Decomposition of the similarity model's integrated SGS dissipation. . . . 130 7.1 Floating-point operation assignments for various mathematical operations. 154 7.2 Estimates of oating-point operation count per mesh point for various SGS models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 vi List Of Figures 4.1 Regions of wavenumber content dened by spectral ltering. . . . . . . . . 48 4.2 Normalized turbulent energy spectrum with wavenumber divisions in spec- tral space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3 Turbulent energy spectrum with accumulation in the vicinity of the LES cuto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.1 Amplitude responses of box lter approximations. . . . . . . . . . . . . . . 75 5.2 Amplitude responses of high-order lters. . . . . . . . . . . . . . . . . . . 77 6.1 Channel geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.2 Test lter amplitude responses. . . . . . . . . . . . . . . . . . . . . . . . . 91 6.3 Model test case: Re =180, Mesh: (32, 32, 33). . . . . . . . . . . . . . . . 93 6.4 Mean velocity proles: Re =180, Mesh: (32, 32, 33). . . . . . . . . . . . . 100 6.5 Mean velocity proles: Re =180, Mesh: (64, 64, 65). . . . . . . . . . . . . 102 6.6 Mean velocity proles: Re =1050, Mesh: (96, 128, 65). . . . . . . . . . . . 103 6.7 Mean velocity proles: Re =1050, Mesh: (128, 128, 97). . . . . . . . . . . 104 6.8 Mean velocity proles: Re =2000, Mesh: (256, 96, 131). . . . . . . . . . . 105 6.9 Mean velocity proles: Re =950, Mesh: (512, 384, 85). . . . . . . . . . . . 106 6.10 RMS velocity proles: Re =180, Mesh: (32, 32, 33). . . . . . . . . . . . . 111 6.11 RMS velocity proles with ltered DNS: Re =180, Mesh: (64, 64, 65). . . 112 vii 6.12 RMS velocity proles: Re =1050, Mesh: (96, 128, 65). . . . . . . . . . . . 113 6.13 RMS velocity proles: Re =1050, Mesh: (128, 128, 97). . . . . . . . . . . 114 6.14 RMS velocity proles: Re =2000, Mesh: (256, 96, 131). . . . . . . . . . . 116 6.15 RMS velocity proles: Re =950, Mesh: (512, 384, 85). . . . . . . . . . . . 117 6.16 Subgrid-scale dissipation eld: Re =180(HR). . . . . . . . . . . . . . . . . 122 6.17 Subgrid-scale dissipation eld: Re =2000. . . . . . . . . . . . . . . . . . . 123 6.18 Forward, reverse, and net subgrid-scale dissipation proles: Re =180(HR). 124 6.19 Forward, reverse, and net subgrid-scale dissipation proles: Re =2000. . . 125 6.20 Mean viscous dissipation proles: Re =180(HR). . . . . . . . . . . . . . . 126 6.21 Mean viscous dissipation proles: Re =2000. . . . . . . . . . . . . . . . . 127 7.1 Mean SGS diusion proles: Re =180(HR). . . . . . . . . . . . . . . . . . 136 7.2 Mean proles of the interscale transfer model's SGS shear stresses. . . . . 140 7.3 Interscale transfer model's predicted dierences in the non-vanishing o- diagonal values of SGS stress. . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.4 Mean velocity predictions using the original formulation and series repre- sentations of the interscale transfer model: Re =180(LR). . . . . . . . . . 149 viii Abstract In large-eddy simulation (LES) various subgrid-scale models have been proposed, all of which attempt to account for the unknown eects of the unresolved scales of turbulence on the resolved ow-eld. The scale-similarity model is one such model, which is formulated using a secondary lter applied to components of the resolved velocity and its products. The scale-similarity model is based on the assumption that the velocities associated with neighboring scales in the ow produce turbulent stresses that are similar in character. The model uses an expression that weights the scales just below the LES cuto in its approximation of the stresses of the unresolved scales. It is well-known, however, that the similarity model fails to accurately predict some of the most fundamental quantities in turbulent ows, perhaps the most important being the global energy transfer and the associated subgrid-scale dissipation. In previous research, additional dissipative terms have been added to the similarity model to improve the model's performance. In the present research, considerations of interscale energy transfer have been used to identify the deciencies of the energy transfer role of the similarity model, specically its inadequate removal of terms contributing energy to the smallest scales and its du- plication of terms producing eects in the largest scales. These considerations are then used in the development of a new model, which shows more favorable characteristics of ix energy transfer. In this approach, partial nonlinear terms are used to decompose the nonlinear transfer present in LES and to formulate an expression capable of removing small-scale production terms depositing energy near the LES cuto. The proposed model is formulated in the same vein as the scale-similarity model, consisting of test-ltered velocities and their products, but the new interscale transfer model oers improvements in its predictions of mean ow quantities and the global energy ux from the resolved to subgrid scales. The current research demonstrates that by implementing this model in a posteriori LES testing of wall-bounded ows, improved LES predictions are possible without the need for additional terms to augment subgrid-scale energy dissipation. x Chapter 1 Introduction The presence of the non-linear term in the governing equations of uid dynamics intro- duces a level of complexity that eliminates much hope of nding a general analytical solution to these equations. As a result, experiments and numerical solutions serve as valuable substitutes where analytical solutions are not available. Applications of numer- ical solution techniques have become increasingly more common and the quality of the predictions routinely obtained from these techniques has improved considerably over the last half-century. Large-eddy simulation is a technique for numerically solving the uids equations, which balances the desire for more accuracy in numerical predictions with the present and foreseeable limitations on computing power. There is a great range of scales of turbulence in most ows encountered in daily life, the smallest of these scales having the greatest demand on computational expense. A promising aspect of LES is that it is not necessary to resolve the small scales of turbulence on a computational mesh nor compute them via the governing equations. All but the small scales of turbulence are resolved and computed in LES. The eects of the unresolved small scales are approximated through the introduction of a model. 1 Appropriate treatment of the small scales of turbulence through subgrid-scale modeling has been the topic of considerable research. Subgrid-scale (SGS) models attempt to account for the unknown exact form of the SGS stress tensor, which consists of velocities and their products that are implicitly ltered in LES. It is this exact form of the SGS stresses which accounts for the unknown, unresolved eects of turbulence in the ow. Scale-similarity models make up a family of SGS models in LES. Similarity models attempt to take advantage of similarities between stresses at dierent velocity scales. These models attempt to account for the unresolved scales of turbulence through the introduction of an additional ltering applied to velocities and their products. However, it is well-known that similarity models fail to accurately predict some of the most fundamental quantities in turbulent ows. The reason for this is that these models do not properly account for the key role of the small unresolved scales of turbulence, specically their reception of kinetic energy from larger scales of turbulence and their subsequent dissipation of this energy. In previous research, additional dissipative terms have been added to scale-similarity models in order for sucient dissipation of energy to occur. In the present research an interpretation of the classic formulation of the similarity model is given with specic regard to its decient role in the energy transfer taking place between the various scales of turbulence. The present research shows that the failure of the similarity model is not in its inadequate energy removal via SGS dissipation, but rather its inadequate mechanism for energy removal in the smallest resolved scales. The goal of the present eort is to develop an alternative to the scale-similarity model that properly accounts for the energy exchange process in turbulence without requiring an additional eddy-viscosity expression. 2 In this research a new model is constructed based on considerations of the interscale energy transfer of the exact form of the SGS stress tensor. A new model is formulated which adopts more favorable characteristics of energy transfer by emulating the role of the implicit mesh lter in LES. In terms of interscale energy transfer, the implicit ltering in LES has the role of removing interscale transfer terms which contribute energy to the smallest unresolved scales of turbulence in the ow. The proposed modeling approach mimics this behavior by removing terms contributing energy to the smallest resolved scales of turbulence in LES. In this approach, an explicit lter is introduced and is used to decompose the nonlinear term in the LES equations. By introducing a SGS model that removes specic partial nonlinear terms responsible for energy production in the smallest scales of LES, the model emulates the role of the exact form of the SGS stress tensor by removing energy from the smallest scales and encouraging the most prominent exchange to occur from the largest to the smallest scales of turbulence. The proposed interscale energy transfer model has been implemented in an incom- pressible, multidomain, pseudo-spectral ow solver. The evaluation of the model has been conducted with this program through a posteriori testing of wall-bounded turbulent ow simulations. These results are compared with results from direct numerical simulations (DNS) and with predictions obtained from implementations of the Smagorinsky model, the dynamic model, and the scale-similarity model. The results for various domain sizes, resolutions, and Reynolds numbers are presented. 3 Chapter 2 Fundamentals of Turbulence Modeling For centuries, mankind has been both amazed and puzzled by turbulence. This fasci- nation has transcended the years and a complete understanding of the phenomena has continued to elude its students. Man's inability to solve this looming problem, so com- monly encountered in daily life, has fueled greater interest in the topic. And so the fascination has grown. Today, a more comprehensive understanding of the problem's complexity exists; however, the phenomenon still persists as the grand unsolved problem of classical physics. As technology has improved, ecient numerical methods for predicting the eects of turbulence in uid ow have come to the forefront of uid dynamics research. This re- search has yielded alternative techniques for attacking the problem numerically, each with its own set of strengths and weaknesses. At present, Reynolds-averaged Navier-Stokes modeling is likely the most widely used solution technique for practical applications, al- though this method can suer from considerable compromises in accuracy. On the other hand, direct numerical simulation oers considerable accuracy but is rarely used in general practice due to its taxing requirements on resolution and computing power. Large-eddy 4 simulation is a promising numerical technique which oers reasonable compromises be- tween accuracy and computing expense. 2.1 Background of Turbulence Turbulence is the irregular motion associated with the nonlinear interactions of vortical structures in uid ow. The local quantities in a turbulent ow are characteristically unsteady and seemingly random. Turbulent ows have a range of coherent structures referred to as eddies. Eddies can be thought of as any type of swirling ow structure generated by vorticity within a ow. The eddies in turbulent ows consist of any number of physical structures including randomly shaped vortex laments and other more iden- tiable structures, e.g. horseshoe vortices, pancake vortices, etc. The eddies present in a ow rapidly change size and shape and interact through stretching, coalescing, splitting, etc. These eddies are by nature time-dependent and three-dimensional and have various physical shapes and scales. There is a continual interaction and energy exchange between these structures. Turbulence stems from the presence of vorticity in a ow and can be generated in plumes, wakes, mixing layers, shear layers, etc. The size and structure of the eddies in a ow depends on various constraints such as the imposed pressure gradient, the shape and size of the ow boundaries, the uid viscosity, etc. Generally the large-scale vortical structures within a ow are directly aected by the geometry and the imposed boundary conditions. Conversely, the small-scale structures are universal in the sense that their 5 general character and time-evolving behavior is typical of small scales in other ows not sharing the same boundary conditions and physical constraints. Richardson [68] suggested that that there is a continuous spectrum of energy ranging from large scales to small scales and an energy cascade that links all eddies of intermediate sizes. This work provided a picture of turbulence indicating that the kinetic energy of turbulent ows is contained within eddies of various sizes. This idea, followed by the pioneering work by Kolmogorov [38], suggested that a turbulent ow would form a hierarchy of vortices where there is an exchange of energy between eddies that is generally directed from larger to smaller scales|the smallest eddies being where the bulk of the kinetic energy of turbulence is ultimately dissipated. The equations of uid dynamics, specically, the time-dependent, three-dimensional, Navier-Stokes equations, contain all of the practical detail describing the turbulent ow of uids. For incompressible ows with constant properties and no body forces, the continuity equation and Navier-Stokes equations take the following respective forms: @u i @x i = 0; (2.1) and @u i @t + @u i u j @x j = 1 @p @x i + @ 2 u i @x 2 j : (2.2) These equations contain the familiar dependent variables of uid ow, specically, pressure (p) and the velocity vector (u i ). The uid properties, including density (), 6 and kinematic viscosity (), also appear in these equations. These equations are well- established and have been used to study problems of varying complexity for more than a century. However, due to the nonlinearity in the equations, a general analytical theory that governs turbulent uid ow continues to elude researchers. As a result, people historically have relied heavily on experiments to extract useful information from specic types of ows. In more modern times, advances in numerical methods and computing capability have oered a viable alternative to conducting costly experiments and these simulations have supplied considerable quantitative information for many types of ows. The eld of computational uid dynamics (CFD) focuses on numerically modeling the Navier-Stokes equations in a discretized domain and numerous studies in this eld have demonstrated that accurate predictions of turbulent ows are possible. 2.2 Reynolds Averaging Turbulent ows are characteristically unsteady and the velocity and pressure variables can be decomposed into a mean component (denoted with an overbar or with a capital letter) and uctuating component (denoted with a prime). Thus, the instantaneous velocity is written as u i =u i +u 0 i =U i +u 0 i : (2.3) Similarly, pressure is denoted as p =p +p 0 =P +p 0 : (2.4) 7 This convention is referred to as the Reynolds decomposition of the ow variables where the mean components are computed as an ensemble averaging or by integrating over some adequately long time interval (T) u (~ x;t) =U i (~ x;t) = 1 T t+T=2 Z tT=2 u i ~ x;t 0 dt 0 : (2.5) The uctuating components of pressure and velocity represent the deviation between their instantaneous and mean quantities and so the time-integration of these uctuating quantities in theory is zero. By implementing the Reynolds decompositions in the governing equations and apply- ing an ensemble averaging satisfying the Reynolds conditions, one obtains the Reynolds- averaged Navier-Stokes (RANS) equations, i.e. @U i @x i = 0; (2.6) and @U i @t + @U i U j @x j = 1 @P @x i + @ 2 U i @x 2 j @u 0 i u 0 j @x j : (2.7) Clearly, the RANS equations are similar in form to the Navier-Stokes equations. The RANS equations contain mean values of dependents variables with the exception of the additional nonlinear velocity term on the right-side of these equations, specically @u 0 i u 0 j @x j : (2.8) 8 This term arises when Reynolds averaging is applied to the decomposed nonlinear term. When multiplied by density, the averaged products of the uctuating components of velocity are referred to as the Reynolds stresses, which are given by u 0 i u 0 j : (2.9) These terms exert stress on the mean ow due to turbulent uctuations. The level of these uctuations depends on the ow Reynolds number (Re =Ud=), where U and d represent mean velocity and a length scale characterizing some important dimension in the ow, respectively. The importance of the nonlinearity in a ow can be quantied by the Reynolds number, which is eectively a ratio of the ow's inertial eects to its viscous eects. In a turbulent ow the Reynolds stress acts in conjunction with the viscous stress where the viscous stress tensor ( ij ) is dened as a product of the dynamic viscosity , where ( =), and the strain-rate tensor (S ij ), thus ij = 2S ij = 2 1 2 @U i @x j + @U j @x i : (2.10) The mean strain-rate tensor is a symmetric tensor that describes the deformation of the uid, i.e. S ij = 1 2 @U i @x j + @U j @x i : (2.11) 9 The eects of pressure acting normal to the surfaces of the uid volume, and the viscous stresses caused by the action of viscosity that arises due to relative motion in a uid, are often combined to form a single stress tensor written as T ij = ij P ij ; (2.12) where ij is the Kronecker delta and is dened as ij = 8 > > > > < > > > > : 1 if i =j 0 if i6=j : (2.13) Both the viscous stress and the turbulence stress are represented by second-order ten- sors where the diagonal components represent stresses acting normal to the surfaces of a uid particle and the o-diagonals represent those acting tangential to the particle's sur- faces. The turbulent stresses account for a ow's irregular and seemingly unpredictable nature. The Reynolds stress tensor is symmetric and thus contains six unknown com- ponents. The introduction of these additional unknowns along with the mean pressure and the three components of the velocity vector exceeds the four available equations. This is known as the closure problem in turbulence. Historically, the topic of turbulence modeling has focused on the development of models capable of approximating the as- pects of turbulence through some creative treatment of the unknowns contained within the Reynolds stress tensor. This approach has taken on varying levels of complexity in turbulence modeling for the RANS equations. 10 In an incompressible ow the thermodynamic pressure cannot be dened. Practically speaking, this is not an issue since it is the gradient of pressure that has the important role in the evolution of the ow. The pressure gradient is computed from the governing equations and this along with information supplied by the boundary conditions can be used to determine the local pressure in the ow. In the equations governing incompressible ow, the dependent pressure variable is interpreted as the mean mechanical pressure (Stokes's assumption), which equals the average of the local normal stresses. This implies that the average normal viscous stress is zero. Thus, pressure can be expressed as P = 1 3 T ii : (2.14) As turbulence stresses arise the normal components of the Reynolds stress tensor can also be expressed as an isotropic part containing turbulent kinetic energy (k), specically 2 3 k = 1 3 u 0 i u 0 i , where k 1 2 u 0 i u 0 i : (2.15) In an incompressible ow, this additional term can be expressed as a contribution to the pressure term, i.e. P = 1 3 T ii u 0 i u 0 i : (2.16) In incompressible ow simulations the pressure term can absorb additional isotropic terms without aecting the pressure gradient nor the evolution of the ow. This is a 11 subtle but somewhat important concept in turbulence modeling because approximations of the Reynolds stress tensor (in RANS modeling) or of the subgrid-scale stress tensor (in LES modeling) are often expressed in such a way that the the trace of the approximate expression does not match that of the exact expression. The approximate expressions can easily be re-expressed to resolve this issue mathematically; however, the models are often left in their primitive form since in practice it is unnecessary to apply any correction to account for the mismatch of the terms since the pressure term can simply 'absorb' the inequality associated with the isotropic components with no net eect to the solution. For convenience, in the presentation of the various models, this same mathematical liberty is taken on occasion in this dissertation. 2.3 Reynolds-Averaged Navier-Stokes Modeling Despite having a set of equations that are well dened, the eld of CFD is not without its challenges. Specically, the key obstacle in the eld has been the struggle to appropriately treat the nonlinear terms in the governing equations. The earliest attempt at developing a mathematical model for the Reynolds stresses was by Boussinesq [5] who applied the eddy-viscosity concept utilizing an analogy between turbulence transport and the process of momentum transport through molecular diusion. The eddy-viscosity relation utilizes length and velocity scales of turbulence to compute an eective viscosity that accounts for the diusive nature of turbulence. The assumption that this additional viscosity can account for the eects of turbulence on the mean ow is not physically accurate, but it leads to a convenient simplication of the mathematics. The general relation showing 12 an assumed eddy-viscosity quantity ( t ) that relates the Reynolds stresses to the mean strain rate is shown below (note that the last term on the right side of Eq. (2.17) has been included to avoid the mathematical inequality of the trace terms in these tensors) u 0 i u 0 j =2 t S ij + 1 3 u 0 k u 0 k ij = t @U i @x j + @U j @x i + 1 3 u 0 k u 0 k ij : (2.17) Prandtl [66] used a technique for relating the Reynolds stress terms to the mean strain rate through an eddy-viscosity relation computed using an algebraic equation that depends on a velocity scale () and a pre-specied turbulence length scale (`), referred to as the mixing length. Models of this type are referred to as the mixing-length models, and are often expressed as t ` @U @y ` `: (2.18) Kolmogorov [38] and Prandtl [67] proposed models that related the velocity scale in the eddy-viscosity term to to the kinetic energy (k) of the turbulence t k 1=2 `: (2.19) This velocity scale can be computed through the solution of a transport equation for turbulent kinetic energy. These models, however, require that the length scale charac- terizing the ow's turbulence be specied in advance. Mixing length models result in an eddy-viscosity that can vary throughout the ow-eld since its velocity and length scales are not properties of the uid. However, arguments have been made that models of this 13 nature are only suited for ows which have single length and time scales (Tennekes and Lumley [80]). Algebraic models, such as the Baldwin-Lomax model [1], and the Cebeci-Smith model [8], implement simple expressions to compute the turbulent viscosity from the ow vari- ables. The turbulent viscosity is then used to compute the unknown Reynolds stresses. On the other hand, one-equations models (such as the Spalart-Allmaras model [75], etc.) have been introduced, which typically compute k using an additional transport equation and then estimate a length scale based on physical considerations. Kolmogorov [39] proposed a model that oered predictions of both length and time scales of turbulence and he introduced what is considered the rst \complete" model of turbulence since it required no prior knowledge of the ow's turbulence characteristics. Kolmogorov's k! model requires the solution of two dierential equations in addition to the governing equations. The k! model has undergone modications and improvements since the work by Kolmogorov, most notably by Wilcox [88]. The present form of the k! model solves transport equations for the turbulence kinetic energy (k) and the specic dissipation rate (!). Once determined, these terms are related to the velocity and length scales in the eddy-viscosity expression. The modern k! model is a popular eddy-viscosity model which allows the additional transport equations to be integrated through the viscous sublayer [8]. Thek" model formulated by Jones and Launder [37] is a related two-equation model that is also widely used in RANS calculations. This model solves transport equations for both turbulent kinetic energy and its dissipation rate ("). As in the k! model, 14 values for the velocity and length scales used for the estimate of t are obtained following the solution of the transport equations. Thus, the eddy-viscosity can be approximated directly as a function of the turbulence kinetic energy and its dissipation rate t k 2 " ; (2.20) where the dissipation rate is dened as, " 2S ij S ij : (2.21) The Reynolds stress model proposed by Launder [42] takes, perhaps, a more direct approach and solves transport equations for the Reynolds stress terms. However, addi- tional unknown terms are introduced in these transport equations and each of these terms must be treated in some manner to achieve closure. A perceived strength of the Reynolds stress model is that the solution of the Reynolds stress equations eliminates the need for an eddy-viscosity or a mixing-length relation. Eddy-viscosity models and Reynolds stress models are included in the category labeled RANS models, referring to the solution of the Reynolds-averaged Navier-Stokes equations to obtain the mean velocity eld. Many of these models were proposed well in advance of one's ability to apply them with any level of convenience because of the limited availability of the computers needed to carry out the calculations. However, with improvements in computing power, these and other models have been implemented with varying levels of success. 15 2.4 Direct Numerical Simulation Eorts in turbulence modeling have varied somewhat depending on the shortcomings of the proposed modeling technique. The key drawback in the implementation of RANS models is that the turbulence model is tasked with supplying a complete description of the eects of turbulence on the mean ow. On the other hand, direct numerical simulation has the capacity to be a much more precise method that takes a more straightforward approach to treating the governing equations. DNS oers a complete description of the turbulence by numerically integrating the exact form of the governing equations. Therefore, no modeling is required in DNS. The challenge of DNS is satisfying its great demand on mesh resolution and therefore computing resources. DNS requires that all relevant scales of turbulence be resolved both temporally and spatially. The eddies in a ow can have a great range of scales with sizes ranging from the largest length scales in the ow domain (i.e. the integral length scales) to the tiny dissipative scales of size () that were theorized by Kolmogorov and can be expressed as a function of viscosity and the dissipation rate of kinetic energy, i.e. = 3 " 1=4 : (2.22) Kolmogorov also derived a velocity scale () associated with these eddies, which takes the form = (") 1=4 : (2.23) 16 Combining these scales of length and velocity, one obtains the time scale == = (=") 1=2 ; (2.24) which is also associated with the smallest eddies in a ow. In DNS, since all relevant scales of turbulence must be resolved, the spatial and temporal resolution that is required is often based on these estimates of the smallest length and time scales. In DNS, the dynamics of the ow requires the resolution of all scales and this generally requires a domain size on the order of the integral length scale with space and time resolution capable of resolving eddies in the range of Kolmogorov's scale. For this reason, DNS generally requires considerable computational expense. To provide a physical example of the resolution expense associated with DNS, the smallest relevant length scales associated with oceanic ows are on the order of one cen- timeter. The sizes of the regions which are often modeled in these applications can vary greatly. Models of ocean current systems can extend thousands of kilometers. Assuming a domain size of just one kilometer, a mesh size consisting of 10 5 cells in each direction, or a total mesh size of roughly 10 15 cells, would be required to fully resolve all relevant structures in the ow. Simulations of this scale are very large and are considered out of reach for most all who wish to use CFD for general research and development appli- cations. Industrial ow applications do not normally involve ow domains on the order of kilometers; however, the required DNS resolution in these problems can be aected greatly by gradient resolution requirements and can scale in such a way that can keep the use of DNS in these applications similarly prohibitive. 17 Spectral methods are often implemented in DNS due to the method's eciency in spe- cic applications and its high accuracy. In spectral methods, Fourier modes associated with the primitive variables are advanced in time according to the governing equations. In DNS, the majority of the wavenumber modes are in the dissipative range (i.e. small scales), which often requires an extremely ne mesh resolution throughout the domain. For typical ow problems this resolution requirement is far too costly from a computa- tional standpoint due to the presence of the very small turbulence structures distributed throughout the ow-eld. The computing requirements of DNS increase greatly with only moderate increases in Reynolds number. The computational requirements for DNS are typically too demanding for the sake of design and development in general engineering ap- plications. Roughly speaking, the computational cost of DNS scales like Re 3 (Pope [65]), which, despite continuing rapid increases in computing capability, will likely prohibit the use of DNS in general application for many decades to come. The eld of CFD has many options for modeling turbulent ows and the best approach depends greatly on the application. DNS oers greater accuracy but is less frequently used in general practice. RANS modeling is more frequently used due to its more relaxed requirements on resolution and its robustness in various applications. LES is a numerical solution technique that bridges the gap between RANS modeling and DNS. As computing resources and capability become more readily available, LES is expected to become much more viable for general applications due to its favorable balance between accuracy and computing requirements. 18 2.5 Large-Eddy Simulation In real ows the large scales of turbulence tend to be chie y responsible for the great- est transport of momentum, heat, etc. LES resolves the largest scales of turbulence in the ow, but relaxes the resolution requirements for the smallest scales. The eects of the smallest scales of turbulence are treated through modeling. In LES, a SGS model accounting for the eects of the small unresolved scales is introduced directly into the governing equations. In a way that is analogous to the derivation of the RANS equations by the decomposition of velocity into mean and uctuating components, the LES equa- tions are obtained by decomposing the instantaneous velocity and pressure into ltered ( u i ) and residual (u 0 i ) components, specically u i = u i +u 0 i ; (2.25) p = p +p 0 : (2.26) The above convention is used throughout the remaining sections of this dissertation. In LES the governing equations undergo a ltering operation corresponding to the eect of the under-resolved mesh. The ltered portions of velocity and pressure are the components associated with the large eddies in the ow resolved by the mesh, while the residual components are those associated with small-scale motions that are unresolved by the mesh. Adopting this nomenclature, the spatial ltering operation (represented by an overbar) can be expressed as a three-dimensional convolution integral as follows: 19 f (~ x) = Z +1 1 G ~ x ~ x 0 ; f ~ x 0 d ~ x 0 : (2.27) In this expression G is a smoothing kernel representing the ltering eect that locally weights the three-dimensional elds of pressure and velocity throughout the domain. The lter width is a free parameter in this expression and is represented by the length scale . Assuming that the ltering operation commutes with dierentiation, the LES equations are obtained by applying the above ltering operation to the instantaneous continuity and momentum equations. The commutative property is satised when a spatially uniform ltering operator is used. The result of this are the equations of mass and momentum transfer (expressed here assuming no body forces) for the resolved eld in LES, which are given as @ u i @x i = 0; (2.28) @ u i @t + @u i u j @x j = 1 @ p @x i + @ 2 u i @x 2 j : (2.29) In these equations the nonlinear term contains products of both ltered and residual velocity components. It is customary to recast the nonlinear term as u i u j = u i u j + (u i u j u i u j ): (2.30) In LES, the bracketed term is referred to as the SGS stress tensor and is given the symbol ij , i.e. 20 ij =u i u j u i u j : (2.31) The SGS stress tensor contains the eects of the nonlinear quantities unresolved by an LES mesh. By decomposing the instantaneous velocity and pressure into ltered and residual components and by applying the ltering to the incompressible form of the continuity and momentum equations, the LES equations for an incompressible ow are obtained as @ u i @x i = 0; (2.32) @ u i @t + @ u i u j @x j = 1 @ p @x i + @ 2 u i @x 2 j @ ij @x j : (2.33) The exact form of the SGS stress tensor contains the unknown eects of the unre- solved scales of turbulence. Turbulence modeling in LES addresses the treatment of this unknown term. A SGS model accounting for the eects of the small unresolved scales of turbulence is typically introduced directly into the governing equations in place of ij . A key aspect of the rationale used in LES modeling relates back to Kolmogorov's theory of self-similarity [38]. Specically, the large eddies in a ow have a character that is distinct and can depend greatly on the ow type and the boundary conditions; whereas, the small scales are eectively universal between many types of turbulent ows. The large eddies, having a unique character among various ow types, are less suited for treatment through modeling and more suited for direct calculation. Conversely, the small eddies 21 in a ow have a character that is common to many dierent types of turbulent ows. Consequently, capturing the eects of small-scale structures more readily lends itself to modeling. The advantage of directly computing the largest eddies and modeling the smaller ones is further realized considering the small "universal" eddies in a turbulent ow demand far-and-away the greatest fraction of the overall computational expense due to their de- manding requirements on spatial and temporal resolution. Avoiding the heavy expense of resolving and computing the small eddies improves the viability of the LES approach. The compromises of LES are that it lacks the precision of DNS and it is usually computa- tionally more demanding than RANS modeling. Of course the benets of LES modeling are that it can oer considerably more accuracy than RANS techniques while avoiding the prohibitive computational expense of DNS. LES is thought by many to have a favorable balance of both accuracy and computational expense and for this reason is expected to see an increase in popularity. LES modeling and subgrid-scale model development have been topics of research since the early eorts in LES applications by Smagorinsky [74]. Eorts to make improvements in this area have increased considerably in recent years. The major challenge in LES modeling has been the formulation of an expression for the subgrid-scale model that captures the important eects of the unresolved scales of turbulence. The current work falls into the category of subgrid-scale model development for LES. 22 Chapter 3 Subgrid-Scale Modeling in LES 3.1 The Smagorinsky Model In LES research there are three traditional categories of SGS models: eddy-viscosity models, similarity models, and mixed-models which combine these two types. A central focus of research in the eld of LES has been the development of a model capable of properly accounting for the role of the unresolved scales of turbulence. Perhaps the most important role of the small scales of turbulence that remain unresolved in LES is their dissipation of turbulent kinetic energy. Consequently, the category of SGS models that have received the most attention are eddy-viscosity models. Eddy-viscosity models refer back to the eort by Boussinesq [5] to relate turbulent transport to molecular transport through an additional viscosity, and to Prandtl's eddy-viscosity approximation for the Reynolds stress tensor in RANS modeling. Upon assuming the SGS stress tensor follows a gradient-diusion process, the modeled SGS stress tensor takes the following form 23 ij =2 t S ij + 1 3 kk ij : (3.1) In this model, t again represents the proportionality constant referred to as the eddy- viscosity; however, in this case it relates the SGS stress tensor to the resolved strain-rate tensor S ij dened as S ij = 1 2 @ u i @x j + @ u j @x i : (3.2) Smagorinsky [74] proposed the rst SGS model for LES that assumed the SGS stresses obey a gradient-diusion relationship. Many related models of this nature have followed suit (e.g. [14], [51], [6], [89], [90]). The model is founded on the idea that the small- scales of motion in turbulence are chie y dissipative and that the energy produced in the resolved scales equals the energy dissipated by the unresolved scales [2]. The role of SGS models of this type is to remove turbulent energy from the ow that, given adequate mesh resolution, would be dissipated by small scales. In the Smagorinsky model an eddy-viscosity term is computed where the turbulent viscosity is again assumed to be a product of a length scale and velocity scale, where t = c s 2 S : (3.3) In this expression, the velocity scale is obtained by taking the norm of the ltered strain- rate tensor, i.e. 24 S = 2 S ij S ij 1=2 : (3.4) Combining this with the length scale completes the estimate of eddy-viscosity. The Smagorinsky length scale ` s is dened as ` s = c s : (3.5) In this expression c s is a constant and is the lter width, which has been estimated as a function of the mesh spacing in each direction (see Deardor [13] and Scotti et al. [73]), i.e. = ( x y z ) 1=3 : (3.6) The value of the Smagorinsky length scale is not universal for all ows and it is adjusted through the Smagorinsky constant c s , which has been reported as having an optimum value in the range (0:05c s 0:25) depending on the Reynolds number (Dear- dor [13], Schumann [72], Bardina et al. [2], Germano et al. [32], Pope [65] ). Van Driest wall damping [81] is traditionally introduced near solid boundaries to account for the reduction in the growth of small scales at the walls (Moin and Kim [57], Piomelli [62]). Thus, the corrected Smagorinsky length scale takes the form ` s =c s 1exp z + 25 ; (3.7) 25 where z + is the dimensionless coordinate indicating distance from the wall, i.e. z + = zu = (see section 6.1 for a more detailed description). The full form of the Smagorinsky model is given as ij =2 c s ( x y z ) 1=3 1exp z + 25 2 2 S mn S mn 1=2 S ij : (3.8) This model has received considerable popularity and has enjoyed some success due mainly to its ability to appropriately capture the global energy ux from resolved scales to subgrid scales (Domaradzki and Liu [22]). The Smagorinsky model has been used successfully in the modeling of transitional and turbulent ows (Germano et al. [32]); however, the model has limited generality since the Smagorinsky constant must be spec- ied with advanced knowledge of the ow characteristics. Additional criticism has come due to the model's inability to account for reverse energy transfer in turbulence and the model's requirement for ad hoc wall damping, which must be introduced as a correction to the model's calculation of a turbulence length scale (Meneveau and Katz [54]). 3.2 The Dynamic Model Germano et al. [32] proposed a model that amounts to a dynamic form of the Smagorinsky model. This model, again, uses the general form for an eddy-viscosity model, i.e. ij =2 t S ij + 1 3 kk ij : (3.9) 26 The dynamic model expands on Smagorinsky's approach by locally computing a coe- cient c appearing in a similar expression for eddy-viscosity: t =c 2 S : (3.10) The goal of this model is to improve upon the key shortcoming of the Smagorinsky model, specically its reliance on the accurate specication ofc s . The model achieves this by dynamically computing the appropriate length scale from the local ow conditions. The value of the constant c is computed and is allowed to vary in both space and time. The proposed method for dynamically computing the eddy-viscosity coecient requires the application of an explicit secondary lter (also referred to as a test lter) represented as a hat, along with the implicit grid lter represented as an overbar. Applying the test lter to the resolved velocity eld allows a sample of the smallest resolved scales to be taken for the sake of modeling the subgrid scales. In the dynamic model's approach, an expression for the resolved turbulence stress tensor is introduced as L ij d u i u j b u i b u j : (3.11) Applying the test lter on top of the implicitly ltered quantities in the exact expression for the SGS stress yields T ij = d u i u j b u i b u j : (3.12) 27 Similarly, the application of the test lter to the exact form of the SGS stress tensor gives b ij = d u i u j d u i u j : (3.13) An expression referred to as the Germano identity is obtained by taking the dierence between these terms, i.e. T ij b ij = d u i u j b u i b u j d u i u j d u i u j = d u i u j b u i b u j : (3.14) Clearly, the above dierence is equal to the expression for L ij . This expression con- tains turbulence stresses associated with the small resolved scales. These stresses are associated with the scales of turbulence resolved by the mesh that would be removed by the application of the test lter. The key to this procedure is that it assumes that the tensors ij and T ij can be accurately represented using similar forms of the eddy-viscosity expression, specically ij =2c 2 S S ij ; (3.15) and T ij =2c b 2 b S b S ij : (3.16) In LES, the width of the implicit grid lter is dened as the spatial separation of the computational grid points. The length scale corresponding to this implicit lter is represented as ; while the test lter's length scale is represented as b . These two 28 length scales are related by the type of test lter that is applied. For the dynamic model, Lund [50] gives the procedure for accurately specifying the test lter width based on the second moment of the lter transfer function. The lter shape is a function of the quadrature used in the approximation of the integration applied during ltering, and the resulting lter width varies accordingly. If quadrature using Simpson's 1/3 rule is applied, an approximation of the box lter (also referred to as the top-hat lter) in physical space is obtained. For this case, the relationship between the test lter width and the computational width is b = 2 . Thus, the length scale of the test lter is twice that of the implicit lter. If the Trapezoidal rule is used to approximate the box lter, then the relationship is b = p 6 . Applying the test ltering operation to the eddy-viscosity expression for ij and sub- tracting this from the analogous expression for T ij gives T ij b ij =2cM ij ; (3.17) where M ij is dened as M ij = b 2 b S b S ij 2\ S S ij : (3.18) By setting the Germano identity equal to the dierence in eddy-viscosity expressions one obtains L ij =T ij b ij =2cM ij : (3.19) 29 Recalling the denition of L ij , this equality places restrictions on the constant c. However, since the equation for L ij above represents multiple conditions that must be satised between tensor components, a single value of the constant c cannot satisfy these requirements without error. Lilly [47] proposed minimizing the error in c by using a least-squares approach. To keep the ratio bounded, additional M ij terms were included in the numerator and denominator of the nal expression. To avoid numerical instabil- ities associated with negative values of c, planar averaging, represented ashi, was also implemented in homogeneous directions within the ow-eld (Germano et al. [32]). Thus, the constant c takes the following nal form (note that the summation implied by the index notation is conducted separately for the numerator and the denominator): c = 1 2 hL ij M ij i hM ij M ij i : (3.20) Adopting this expression for the Smagorinsky constant and adding this to the original expression for the SGS stress tensor gives the following form for the dynamic model: ij =2 1 2 hL ij M ij i hM ij M ij i 2 S S ij ; (3.21) where L ij ;and M ij are dened above. Thus, the full detailed expression for the dynamic model is ij = D d u i u j b u i b u j b 2 b S b S ij 2\ S S ij E b 2 b S b S ij 2 \ S S ij 2 2 S S ij : (3.22) 30 The dynamic model is considered by many to be an industry standard for LES ap- plications and is often used as a benchmark for model performance in present-day LES research. The dynamic model, however, suers from some key weaknesses, specically the well-known drawback of eddy-viscosity models, which tend to predict stresses that correlate poorly with stresses in actual turbulence (Meneveau and Katz [54]). Implemen- tation of the model can also lead to numerical instabilities when reverse energy transfer is predicted (Germano [32], Park et al. [61]). Negative values of the model's constant cor- respond to negative turbulent viscosity representing energy transfer from small to large scales (referred to as backscatter). Backscatter, by denition is a non-dissipative phe- nomenon; however, the dynamic model is inherently dissipative in nature. The numerical averaging that is needed to avoid the development of numerical instabilities undermines the intended nature of the model's local length scale calculation and eliminates its ability to predict reverse energy transfer. The above restrictions can eliminate the local nature of the constant computed by the dynamic procedure and limit the model's applicability in more general geometries. An approach implemented by Carati et al. [7] uses ensemble averaging in place of aver- aging in homogeneous directions. Aside from the obvious drawback of requiring multiple realizations of the oweld, Meneveau and Katz [54] report that another disadvantage of this averaging technique is that it too can result in spatial smearing of local values of the Smagorinsky coecient. A Lagrangian approach proposed by Meneveau et al. [55], is a model that has attempted to address this concern. Nonetheless, despite its weak- nesses, the dynamic model has enjoyed much attention and some success. The model has 31 shown good predictions of transitional and fully turbulent ows and has won considerable popularity in the eld of SGS modeling for LES [65]. 3.3 The Scale-Similarity Model By numerically solving the governing equations on a mesh that is incapable of fully resolving all relevant scales of motion, an implicit ltering process becomes intrinsic to LES solutions. The result is a ltering operation applied to all terms in the governing equations. The full detail of this implicit ltering cannot be known. Thus, the fully resolved velocity eld cannot be recovered from an under-resolved simulation. On the other hand, as demonstrated in the development of the dynamic model, a prescribed ltering operation can be applied to a computed velocity eld and the remaining ltered component along with its residual counterpart can be used for the sake of modeling the subgrid scales. Much like the implicit mesh lter, a prescribed test ltering operation can also be expressed in the form of a convolution integral with a selected smoothing kernel G, i.e. b f (~ x) = Z +1 1 G ~ x ~ x 0 ; b f ~ x 0 d ~ x 0 : (3.23) In this operation, the function G denes the lter shape, and b is the lter width. In LES applications, this function is normally selected to give a low-pass eect and so typically the ltered eld is smoother than the unltered eld. The concept of the rst scale-similarity model was introduced by Bardina et al. [2]. Scale-similarity models make up an alternative family of SGS models. These models 32 attempt to account for the unresolved scales of turbulence through the introduction of an additional ltering applied to velocities and their products. This decomposition of resolved scales is used to approximate the form of the SGS stress. The formulation of this model seeks to make use of the perceived similarities associated with stresses at neighboring velocity scales, specically those scales located just below and above the LES cuto. The approach requires some method for separating the resolved velocity into components associated with large and small-scale structures. This is accomplished through the application of the test lter operation. The approach of scale-similarity models is rationalized if one assumes that the largest of the subgrid scales account for the most important SGS interactions and these scales correlate well with the smallest resolved scales. In proposing the rst of such models, Bardina et al. [2] followed Leonard's [43] de- composition of the SGS stress. Leonard decomposed the SGS stress into component con- tributions from the SGS Reynolds stress (u 0 i u 0 j ), cross stress ( u i u 0 j +u 0 i u j ), and Leonard stress ( u i u j u i u j ) terms. Bardina's model employs ltering to approximate the unknown velocity of the subgrid eld in LES, represented as u 0 i , and dened as u 0 i =u i u i ; (3.24) with the ltered uctuating velocityu i u i computed as u i u i . Approximating ltered velocity products with products of ltered velocities, Bardina obtained the approximation u i u j u i u j for the sum of the cross stress and Reynolds stress terms. By including the 33 Leonard stress and allowing for an arbitrary lter shape, the nal form of the scale- similarity model attributed to the work of Bardina et al. [2] is written as cs ij = u i u j u i u j : (3.25) Bardina's expression amounts to reconstructing the eect of the mesh lter and approx- imating the the full velocity eld (u i ) with the resolved velocity component ( u i ) that is obtained from LES, i.e. u i u i : (3.26) Adopting a Gaussian type test lter, Bardina et al. [2] demonstrated in a priori analy- ses of decaying homogeneous isotropic turbulence and in sheared homogeneous turbulence, that the similarity model was capable of producing good tensor, vector, and scalar level predictions of the SGS eld. The selected form of the Gaussian kernel is given as G ~ x ~ x 0 ; = r 6 1 ! 3 exp 6 2 ~ x ~ x 0 2 : (3.27) In these predictions the model showed improvements over the Smagorinsky model. Another noted benet of the similarity model was its ability to predict SGS stress with eigenvectors unaligned with those of the mean strain rate, thus allowing for the possibil- ity of reverse energy transfer, which is common in turbulence. It was noted by Bardina, however, that the similarity model was decient in one key aspect, its dissipation charac- teristics. For this reason a mixed model was formulated, which combined the similarity 34 model with an eddy-viscosity expression. Testing of this model showed that the good correlations predicted by the similarity model were still possible after incorporating the eddy-viscosity expression. Alternative scale-similarity models were also considered by Meneveau and Lund [56], Liu et al. [48], Meneveau [58], and Winckelmans et al. [90]. The Gaussian lter used by Bardina removes nearly all content from the largest wavenumbers and the activity in the small wavenumbers is only reduced. The research by Liu et al. [48] suggested that the Gaussian lter used by Bardina, which applied some weighting at all points of the mesh, resulted in the removal of much of the uctuating behavior of the SGS stresses. Comparisons were made with measured uctuations from actual stresses and this research showed that the similarity model in the form constructed by Bardina did not generate adequate RMS uctuation levels of the SGS stresses. Liu concluded that Bardina's formulation of the similarity model would require too great of a dependence on a constant of proportionality in order to provide acceptable uctuation levels in the predicted stresses. The initial development of the similarity model was extended by Liu et al. [48] using experimental data in the far-eld of a turbulent jet. The similarity model was refor- mulated given that the implicit lter associated with the LES mesh cannot be known explicitly. Liu et al. [48] followed a similar approach, but applied a wider test lter and also included an adjustable constant. The alternative form of the similarity model is expressed as s ij =c L d u i u j b u i b u j ; (3.28) 35 where c L is the model's constant. The research by Liu et al. [48] showed that implementing a smooth test lter gave better results than using a sharp spectral cuto. Liu's discrete approximation of the box lter was implemented with a test lter width that was twice the computational mesh spacing b = 2 . The box lter is represented with a lter function dened as G ~ x ~ x 0 ; = 8 > > < > > : 1= 0 if ~ x ~ x 0 < =2 otherwise : (3.29) Liu's form of the similarity model adopted an approximate form for this test lter. This implementation is considered the functional form of the similarity model. Citing a similar level of uctuation in the modeled and measured stresses, Liu et al. [48] elected to conduct model testing using a constant of c L =1. This form of the similarity model amounts to an approximation for the SGS stress tensor equal to the resolved turbulence stress expression popularized by Germano et al. [32] in the development of the dynamic model. Liu et al. [48] addressed the importance of the lter shape in the application of the model. A box lter, Gaussian lter, and a spectral cuto lter were implemented and the results were compared in the study. Comparisons of the SGS energy ux were also made. Testing showed that the model could produce good correlations between predicted and measured stresses when the box lter and the Gaussian lter were implemented. Much weaker correlations were reported when a sharp spectral cuto was used. Conclusions drawn from this testing were that the similarity model, when implemented with a box lter or Gaussian lter, continued to give good correlations with actual stresses and the 36 model reasonably captured the characteristics of backscatter in the ow. However, Liu et al. [48], and previously Bardina et al. [2], noted that the model in its proposed form was not adequately dissipative. In an eort to improve the model's dissipation characteristics, a backscatter control function was later added and a dynamic procedure for computing the appropriate value for c L was also tested. An extensive consideration of this model's constant was conducted by Cook [12] who used a box lter implementation. The optimum value of the constant was investigated by requiring that the model generate the correct amount of SGS turbulent kinetic energy. Optimized values of the constant were tabulated by Cook and in all cases the values were of order unity. This showed that the value ofc L =1 initially assumed by Liu et al. [48] was in the recommended range, specically c L 0.6-1.7. However, it was also shown that the optimum value for the coecient varied depending on Reynolds number, mesh resolution, and lter width. Cook also provided a formula for estimating the value of c L , but the method requires information that is typically not known in advance of a ow solution. Dynamic procedures for computing c L have also been proposed (Liu et al. [48], Cook [12]). However, the formulations of these models have typically been altered in some way in order to provide improvements in the model's ability to remove energy from the turbulence. One such modication used extensively in LES model development has been the inclusion of an additional viscosity term to account for the dissipative eects of the subgrid scales. 37 3.4 Mixed Models Perhaps the most important feature of a SGS model is that it accurately account for global SGS dissipation. In LES this implies a correct representation of the net energy ux from the resolved to the subgrid scales (Domaradzki and Saiki [29], Meneveau and Katz [54]). This feature is typically considered necessary to obtain good predictions of the important mean statistical quantities of the ow such as mean velocity and RMS values of the velocity uctuations. However, models which demonstrate good dissipative characteristics such as eddy-viscosity models, often yield poor correlations with actual SGS stresses. The SGS stresses predicted by similarity models have been found to correlate well with exact stresses in a priori analyses (McMillan and Ferziger [53], Bardina et al. [2], Liu et al. [48]). However, similarity models have been found to signicantly under-predict SGS dissipation and fail to accurately predict mean and RMS quantities in actual turbulent ows. Conversely, eddy-viscosity models have succeeded in accurately predicting SGS dissipation, but have shown poor correlations with actual stresses. A key dierence in the form of similarity models compared to many eddy-viscosity models is that similarity models do not guarantee an alignment of the principal axes of the SGS stress and the mean strain rate (Bardina [2]). Thus, local realizations of SGS dissipation (computed as " sgs = ij S ij ), when predicted by the model, can have either positive or negative values. This is a desirable characteristic; however, in practice similarity models have proved inadequate in LES predictions due to their inability to properly account for SGS dissipation (Liu et al. [48], Bardina et al. [2], Domaradzki et al. [25]). For this reason, it has been postulated that similarity models cannot account 38 entirely for the dissipative nature of turbulence. This notable deciency of similarity models has prompted researchers to combine similarity models with eddy-viscosity models (referred to as mixed models) or to introduce other modications in eorts to improve predictions in LES. Mixed models (e.g. [2], [85], [86], [91], [45], [92]), which combine similarity and eddy- viscosity models, have been proposed and implemented with the intent of beneting from the positive attributes of each class of model, namely the dissipation of eddy-viscosity models along with the similarity model's ability to accurately predict real stresses. Mixed models have been formulated by simply combining an eddy-viscosity expression with the similarity model giving the following form m ij = d u i u j b u i b u j 2 t S ij : (3.30) Specic variations of the mixed-model formulation have been proposed and tested (Bar- dina et al. [3], Piomelli et al. [64], Zang et al. [94], Vreman et al. [84]). The results of testing have shown that mixed models can demonstrate the favorable dissipative eects shown by eddy-viscosity models along with the good predictive ability of similarity models. These simulations have also demonstrated some capacity to pre- dict the presence of backscatter. Backscatter is the energy exchange between scales of turbulence that is directed from small to large scales of turbulence, which is in reverse of the general direction of energy exchange. This is a well-known feature in turbulence and its complete elimination can be a detriment to an LES model's predictive capability. Backscatter predictions are expected from similarity models, but are typically problematic 39 for eddy-viscosity models. Similarity models can generally exhibit backscatter without issue. Research conducted by Kornev et al. [40] has shown that, like eddy-viscosity mod- els, predictions of backscatter while using mixed-models can also lead to issues associated with numerical instabilities. In mixed-model simulations of homogeneous turbulence with a mean shear, the cor- relations between real and predicted stresses reported by Bardina, et al. [2] were equal to, and sometimes slightly higher than, those reported for separate implementations of the similarity model and the Smagorinsky model. However, in mixed models it is well- known that the eddy-viscosity expression provides the energy dissipation necessary in LES. Eddy-viscosity models alone have shown comparable levels of success with regard to predicting bulk quantities of the ow and so the rationale for having the additional complexity of a similarity expression is not well supported. Mixed models suer from the same issues as traditional eddy-viscosity models (i.e. spatial or temporal averaging for the sake of numerical stabilization, advanced specication of model constants, ad hoc treatment near walls, etc.) and are complicated further by the introduction of an ad- ditional length scale from test ltering. Furthermore, as discussed by Domaradzki [19], similarity models can disrupt the natural energy exchange expected in turbulence. Thus, for a mixed model to have success the action of the eddy-viscosity term must overcome this detriment. For these reasons, mixed models have not been fully embraced as a viable alternative for SGS modeling. The shortcomings of both similarity models and eddy- viscosity models has motivated further work in the eld of SGS model development for LES. 40 The failure of the similarity model without an additional eddy-viscosity term under- scores its inadequacy for use in LES. Considerations of interscale energy transfer will show that this failure is due to the model's inability to properly account for the key role of the small unresolved scales of turbulence, specically the reception of kinetic energy from resolved scales and the subsequent dissipation of this energy. The goal of the current eort has been to develop an alternative model, created in the same vein as the similarity model, that properly accounts for the energy exchange process in turbulence and elimi- nates any need for an additional dissipative term. It will be shown that the new model has a functionally simple formulation, and much like the similarity model, is constructed by applying a test-lter to the resolved velocity components and their products. How- ever, the proposed model is formulated based on considerations of the interscale energy transfer associated with the exact form of the SGS stress tensor. The model has been implemented and a posteriori testing has been conducted in wall-bounded turbulence. Predictions by this model are compared with those obtained by implementations of the dynamic model and the similarity model. 41 Chapter 4 Interscale Energy Transfer in Turbulence 4.1 The Energy Cascade in Spectral Space In the present research an alternative to the similarity model is developed. Like previously discussed similarity models, the proposed model consists solely of products of ltered velocity components; however, the present model draws upon similarities associated with energy transfer occurring near the LES cuto rather than on velocity scales and SGS stresses. The present model is derived from these considerations and its properties of energy transfer are compared with those of the scale-similarity model. To foster the discussion of this topic, interscale energy transfer is introduced in the present section. Energy transfer between the various scales of turbulence is a fundamental process in turbulent ows that can be traced back to the presence of the nonlinearity in the governing equations. Richardson's theory [68] of the energy cascade where energy is transferred from the large eddies to small eddies, is a well-known concept. As a result of this process, turbulent ows consist of a hierarchy of eddies that can have great ranges of length scales and time scales. The energy exchange between turbulent eddies theorized 42 by Richardson is often approached in spectral space. The energy spectrum expressed in terms of wavenumbers gives a global description of the energetics of turbulent structures. Localized eddies of varying sizes contribute to this global representation of turbulence in wavenumber space. The general direction of energy transfer in turbulence is from small to large wavenum- bers so this most prominent exchange process is referred to as forward transfer. However, reverse transfer also occurs when the energy of large wavenumber modes causes increases in the energy of small wavenumber modes. Studies of the energy exchange between scales of turbulence in shear ows have been carried out by Domaradzki [24] and Domaradzki et al. [28]. Related studies in isotropic turbulence were conducted by Domaradzki and Rogallo [26], [27], Yeung and Brasseur [93], and Domaradzki et al. [23]. Included in this research have been studies into the exchange of energy across various length scales. The methodology used for the analysis of interscale energy transfer is outlined by Domaradzki and Rogallo [27], Domaradzki et al. [24], and Domaradzki and Liu [22]. The concept of interscale transfer provides a method for both evaluating existing models and developing new ones. This is a key consideration in the development of the present model. In three-dimensional space the Cartesian coordinate system is represented as x i = (x 1 ; x 2 ; x 3 ): (4.1) The Fourier decomposition of a ow-eld leads to a description of the scales of turbulence in terms of wavenumbers. In three-dimensional space a wavevector is expressed as 43 k i = (k 1 ; k 2 ; k 3 ): (4.2) The components of the vector are computed as follows, where L i represents the corre- sponding length scale in physical space: k i = 2=L i : (4.3) A turbulent structure has a wavenumber corresponding to its average size in physical space. Thus, a wavenumber in three-dimensional space is represented as the norm of the wavevector, specically k = k 2 1 +k 2 2 +k 2 3 1=2 : (4.4) A velocity eld in physical space can be represented as a Fourier expansion where (c m ) represents Fourier coecients corresponding to each wavenumber mode, i.e. u j (x;t) = X m c m exp (ik m x): (4.5) The Navier-Stokes equations can be expressed in spectral space by taking the Fourier transform of these equations. This form is written as @ @t u n ~ k;t =k 2 u n ~ k;t +N n ~ k;t ; (4.6) 44 whereN n ~ k;t contains both the pressure and nonlinear eects in the ow and is rep- resented as N n ~ k;t =i P njm 2 ~ k Z u j (~ p;t)u m ~ k~ p;t d~ p; (4.7) and P njm is a given wavenumber function (see Batchelor [4] and Lesieur [46]), and is dened as P njm =k j nm k n k m k 2 +k m nj k n k j k 2 : (4.8) Note that variables expressed in spectral space are distinguished from those in physical space by their arguments (e.g. ~ x or ~ k). The energy amplitude is the contribution to turbulent kinetic energy from specic wavevectors, and is expressed as 1 2 u ~ k;t 2 : (4.9) The equation governing spectral energy content in turbulence is obtained from the spectral form of the Navier-Stokes equations and can be written as @ 1 2 u ~ k;t 2 @t =2k 2 1 2 u ~ k;t 2 +T ~ k;t : (4.10) This equation describes the transfer of kinetic energy associated with mode ~ k. In this equation T ( ~ k;t) is the nonlinear energy transfer and is dened by 45 T ~ k;t =Im P njm 2 ~ k u n ~ k Z u j (~ p)u m ~ k~ p d~ p ; (4.11) where () represents a complex conjugate. This term originates from the nonlinear term N n ~ k;t in the Navier-Stokes equations and embodies the important role of the con- vective term and its eects on the evolution of energy in turbulence. Specically, this term describes the advection of energy between dierent modes that occurs due to the nonlinear interaction of eddies of various wavenumbers in the ow. Energy contributions to an arbitrary mode ~ k arise from the various interactions involving modes~ p and~ q, where ~ q = ~ k~ p. These interactions are restricted to modes~ p, which form a triangle with modes ~ k and~ q in vector space. The integration in the triad interactions accounts for the eects of all possible interactions in uencing mode ~ k, and these interactions account for the continuous coupling of all modes in turbulence. The importance of the locality of the interacting modes becomes apparent in these considerations. In spectral space the energy equation can be expressed in such a way that kinetic energy transfer is associated with wavenumber ranges. Given that the energy transfer terms are rooted in the nonlinear term in the momentum equation, the origin of these interscale transfer terms can be traced directly to dierent components of the nonlinear term. These individual nonlinear components have specic physical interpretations and are referred to as partial nonlinear terms. 46 4.2 Formulation of Partial Nonlinear Terms In the study of interscale energy transfer, the nonlinear term is decomposed to explicitly show the energy transfer eects from nonlinear interactions between dened wavenumber regions. The kinetic energy transfer associated with the various scales of turbulence can be identied by considering discrete wavenumber ranges (e.g. of size k) contributing to the nonlinear term. For illustration, a ow's energy spectrum can be divided into three discrete wavenumber regions denoted by R 1 ; R 2 ; and R 3 , and dened as: R 1 : 0k<k 1 ; R 2 : k 1 k<k 2 ; (4.12) R 3 : k 2 k<1: These regions can be described as containing the smallest, intermediate, and largest wavenumbers, respectively. For the present work, R 3 is taken to represent the range of scales that would remain unresolved by an LES mesh, which generally includes the smallest most dissipative scales of turbulence. Thus, only the content of R 1 and R 2 is computed in LES. The intermediate region, R 2 , represents the smallest resolved scales of turbulence, specically a band of wavenumbers residing adjacent to, but below, band R 3 . In spectral space the truncations needed to construct these regions can be carried out through treatment of the Fourier coecients. 47 Expressed in terms of a ltering operation, the transfer function associated with sharp spectral ltering has the following relationship to the lter's cuto wavenumber (k c ): G (k) = 8 > > < > > : 1 if kk c 0 if k>k c : (4.13) Plots of amplitude response using sharp spectral ltering with two dierent cuto wavenumbers is shown in the gure; the hypothetical regions (R 1 ;R 2 ;R 3 ) constructed in wavenumber space are also shown. Figure 4.1: Regions of wavenumber content dened by a spectral lter's amplitude re- sponse in wavenumber space. If a fully resolved ow-eld is assumed, then through treatment of the Fourier coe- cients thei th andj th components of the velocity eld can be truncated to contain regions R p and R q , respectively. In this case, R p and R q are taken to represent any one of the 48 three dened regions in wavenumber space. Velocity scales belonging to these regions are denoted simply with superscripts 1, 2, or 3. These truncations of the velocity elds are expressed as u p i ~ k 0 ;t = 8 > > < > > : u i ~ k 0 ;t if ~ k 0 2R p 0 otherwise ; (4.14) and u q j ~ k 00 ;t = 8 > > < > > : u j ~ k 00 ;t if ~ k 00 2R q 0 otherwise : (4.15) Partial nonlinear terms represented asN pq n ~ k;t , which account for separate contri- butions to the full nonlinear termN n ~ k;t , are introduced to describe the interactions between modes in the specic wavenumber bands. The divisions of these nonlinear terms are represented as N pq n ~ k;t = 8 > > > > < > > > > : N n u i ~ k 0 ;t ; u j ~ k 00 ;t o if ~ k 0 2R p and ~ k 00 2R q 0 otherwise ; (4.16) where the spectral representation of the nonlinear term is shown asN (), and its argu- ments are the velocities (taken here as truncations in spectral space) indicated on the right side of Eq. (4.16). Thus, the full nonlinear term has been decomposed according to its contributing velocity components. 49 In the momentum equation, these partial nonlinear terms describe the eects of the velocity eld interacting with itself. In fully resolved turbulence, these interactions can generate eects across the entire wavenumber spectrum. Restrictions on the contributing interactions require that the lengths of the interacting wavevectors ~ k; ~ k 0 ; ~ k 00 form a triangle in spectral space. For this reason the partial nonlinear terms and their corre- sponding terms in the energy equations are often referred to as \triad interactions." An additional superscript m can be added to identify the region of wavenumbers represented by (R m ) that is aected by these triad interactions, i.e. N pqm n ~ k;t = 8 > > > > < > > > > : N pq n ~ k;t if ~ k 2R m 0 otherwise : (4.17) The parameters represented by R p ; R q ; R m appearing in the denitions above can correspond to any of the three wavenumber regions described previously (R 1 ; R 2 ; R 3 ). Individually, each partial nonlinear term represents the eects on modes in R m due to the complex nonlinear interactions between modes in R p with modes in R q . Assuming the ow is fully resolved numerically, the sum of all the partial nonlinear terms accounts for all nonlinear eects in the ow. The full list of partial nonlinear terms that gives a complete description of the nonlinear eects in the governing equations is given as N 111 ;N 121 ;N 131 ;N 211 ;N 221 ;N 231 ;N 311 ;N 321 ;N 331 ; N 112 ;N 122 ;N 132 ;N 212 ;N 222 ;N 232 ;N 312 ;N 322 ;N 332 ; N 113 ;N 123 ;N 133 ;N 213 ;N 223 ;N 233 ;N 313 ;N 323 ;N 333 : (4.18) 50 Note, these terms have been organized so that the rst row includes terms aecting modes in R 1 , the second row includes terms aecting R 2 , and similarly for R 3 . For a full DNS solution we have theorized thatR 3 represents the band of wavenumbers that would remain unresolved in an LES simulation, specically those scales above the LES cuto that are removed by the implicit ltering eect in LES. Consequently, in LES all terms involving nonlinear interactions with modes from R 3 , or having an eect in R 3 , are eliminated. Thus, of the 27 terms shown, the application of the implicit ltering from the LES mesh and the subsequent removal of band 3, removes all but eight 2 3 of these terms. Assuming all eects from the third band of wavenumbers are removed by the implicit mesh lter, the eight partial nonlinear terms that remain embody all the nonlinear eects in LES (assuming the absence of a SGS model). The eight remaining terms are given as N 111 ;N 121 ;N 211 ;N 221 ; N 112 ;N 122 ;N 212 ;N 222 : (4.19) An analog of these terms derived in spectral space can be given in physical space assuming an appropriate lter can be used to distinguish between physical scales. Again, by denoting the implicit ltering applied in LES with an overbar and the test ltering operation with a hat, one obtains the following velocity scales and their corresponding wavenumber regions: 51 R 1 : u 1 i = b u i ; R 2 : u 2 i =u 0 i =u i b u i ; R 3 : u 3 i =u 0 i =u i u i : (4.20) These velocity scales are decomposed with a progressive level of ltering and corre- spond to the three regions discussed above. Given these divisions of the velocity eld's wavenumber content, the nonlinear term is decomposed as @ @x j u i u j = @ @x j X m X p u p i X q u q j ! m : (4.21) The above expression represents the full nonlinear term where the superscriptsp,q, andm identify the regions in spectral space and their corresponding translation to physical scale. Assuming the same lter is applied to the velocities and their products, the individual components that this expression is comprised of are written simply as N pqm = @ @x j u p i u q j m : (4.22) Expressions of the above form are the physical representations of the partial non- linear terms decomposed through a general ltering operation. These terms shall be distinguished from N pqm , which will be used to identify the components of stress con- tained inN pqm . Specically, the decomposed stress components present inN pqm are given as N pqm = u p i u q j m : (4.23) 52 The quantitiesN pqm and N pqm are of course related through the divergence operator present in the governing equations. In the present work, specic stress components are incorporated into a SGS model and due to the presence of the divergence operator in the LES equations the implementation of the model is equivalent to also introducing the partial nonlinear terms directly into the LES equations. Likewise, the energy transfer terms associated with these stresses are introduced in the energy equation. For this reason, any reference to adopting partial com- ponents of stress into a SGS model is taken to interchangeably describe the introduction of the corresponding partial nonlinear terms in the LES equations, and the correspond- ing interscale energy transfer terms in the energy equation. In the present context, the introduction of one of these terms also implies the presence of the other two. The various contributions from the nonlinear terms to the three regions of interest are given as R 1 : u p i u q j 1 = d u p i u q j ; R 2 : u p i u q j 2 =u p i u q j d u p i u q j ; R 3 : u p i u q j 3 =u p i u q j u p i u q j : (4.24) Again, if a simulation is fully resolved numerically, the sum of all 27 nonlinear terms would contain all the nonlinear eects in the ow. Individually, the partial nonlinear terms represent the eects on modes in R m due to modes in R p interacting with those in R q . However, all terms involving interactions with modes from R 3 , or producing an eect in R 3 , are removed by the grid lter. In the presence of the implicit lter in LES and the subsequent removal of information in R 3 , a fraction of these terms are known in 53 LES. The small-scale velocity components u 3 k are unknown, as are the partial nonlinear terms contributing to R 3 . Only the partial nonlinear terms combining the explicit test lter operation with the known test-ltered components ( b u i ) and residual components (u 0 i ) of velocity are known. These terms can be computed directly from the resolved eld once a test ltering operation has been dened. Omitting the divergence operator, the eight terms that can be computed directly from the resolved eld in LES are given as N 111 = d b u i b u j ; N 112 = b u i b u j d b u i b u j ; N 121 = d b u i u 0 j ; N 122 = b u i u 0 j d b u i u 0 j ; N 211 = d u 0 i b u j ; N 212 =u 0 i b u j d u 0 i b u j ; N 221 = d u 0 i u 0 j ; N 222 =u 0 i u 0 j d u 0 i u 0 j : (4.25) It can be shown that the sum of these partial components of stress equals the resolved term present in the LES equations (i.e. u i u j ). By constructing these terms with physical lters, the presence of their partial nonlinear terms in the momentum equation has a direct role in the nonlinear transfer of momentum and energy in LES, thus representing the physical analog of the spectral space variablesN ( ~ k;t) andT ( ~ k;t). It is the consideration of these terms that is emphasized in the formulation of a new SGS model, specically, considerations of the eects of these terms and their interpretation pertaining to energy transfer. Specic nonlinear terms have been incorporated into the proposed model based on their role in energy transfer. While the energy transfer eects are emphasized in the formulation of the model, the implementation of the model is of course carried out in the momentum equation through the introduction of specic partial nonlinear terms 54 formulated in physical space. The role of these terms in energy transfer is discussed in the following section. 4.3 Interscale Energy Transfer Terms The partial nonlinear terms describe the complex nonlinear interactions between wavenum- ber modes, but limited information is known about how local ow conditions respond to these interactions. It is well known, however, that the net result of these interactions is an extraction of energy from the mean ow. It is also known that larger scales of turbu- lence generally evolve independent of viscous eects while the smaller scales of turbulence tend to extract energy from the larger scales and the energy passed along the cascade is eventually dissipated by viscosity (Batchelor [4]). The energy exchange resulting from the interactions described by the partial nonlinear terms are expressed as a series of interscale energy transfer terms (represented simply as T pqm ). Thus, each partial nonlinear term N pqm ( ~ k;t) appearing in the momentum equation has a corresponding term representing interscale energy transfer T pqm ( ~ k;t) present in the energy equation. Expressing the in- terscale energy transfer term as a function of its corresponding partial nonlinear term, one obtains T pqm ~ k;t =Re h u n ~ k;t N pqm n ~ k;t i : (4.26) The eects of the nonlinear interactions can conveniently be expressed in terms of their role in the energy equation. The interscale energy transfer terms dene the energy exchange between modes. Specically, these triads describe the rate of energy change in 55 modes residing in R m due to interactions of modes in R p with those in R q . Given this breakdown of activity, each of these terms represents a rate of energy exchange from the complex triad interactions between the discrete wavenumber bands. Studies of locality conducted by Domaradzki and Rogallo [27], Zhou [95], and Domaradzki and Carati [20] have shown that the interactions of eddies occur most prominently between neighbor- ing bands. Thus, energy exchange is more prevalent for local interactions, specically, interactions between modes having similar scales. The full energy spectrum of a turbulent ow contains the kinetic energy of all wavenum- ber modes of turbulence present in the ow. A sample schematic of a turbulent energy spectrum is shown in Fig. 4.2. This spectrum depicts the distinct wavenumber regions dened by ideal sharp spectral ltering. Figure 4.2: Normalized turbulent energy spectrum with wavenumber divisions in spectral space. 56 The interscale energy transfer terms describe the exchange occurring between the dened regions of the spectrum. If these wavenumber regions are considered then the triad expressions represent a total of 27 dierent terms, each corresponding to one of the partial nonlinear terms given previously. Each of these terms has some unique contribution to the evolution of the kinetic energy of turbulence and the list is given as: T 111 ; T 121 ; T 131 ; T 211 ; T 221 ; T 231 ; T 311 ; T 321 ; T 331 ; T 112 ; T 122 ; T 132 ; T 212 ; T 222 ; T 232 ; T 312 ; T 322 ; T 332 ; T 113 ; T 123 ; T 133 ; T 213 ; T 223 ; T 233 ; T 313 ; T 323 ; T 333 : (4.27) These terms embody the eects of all nonlinear energy transfer originating from the full nonlinear term (u i u j ). Each term represents a rate of energy exchange from triad interactions between the discrete bands of wavenumbers. Given that the energy exchange is more prevalent for local interactions, specically interactions between modes not having a wide separation of scales, terms such asT 331 ,T 133 , etc. are often small relative to terms such as T 122 . However, this is not always the case since the relative magnitudes of these terms can depend greatly on the size and arrangement of the discrete wavenumber regions. An important nding with regard to LES was that of Domaradzki et al. [23] who showed that the majority of the energy that is transferred to the subgrid scales originates from the length scales just above the mesh cuto, specically length scales within one oc- tave above (<`< 2). Studies of this nature have supported the picture of the energy transfer process provided by Richardson [68], and despite the presence of backscatter, the transfer of energy is generally from modes with small wavenumbers to those with large wavenumbers. 57 The energy transfer terms known in LES and associated with ( u i u j ) are T 111 ; T 121 ; T 211 ; T 221 ; T 112 ; T 122 ; T 212 ; T 222 : (4.28) Conversely, the interscale energy transfer terms associated with the unknown nonlinear term (u i u j ) are given as T 111 ; T 121 ; T 131 ; T 211 ; T 221 ; T 231 ; T 311 ; T 321 ; T 331 ; T 112 ; T 122 ; T 132 ; T 212 ; T 222 ; T 232 ; T 312 ; T 322 ; T 332 : (4.29) In a fully-resolved simulation (i.e. DNS), the evolution of turbulent energy in the ow leads to a net energy exchange occurring generally from R 1 and R 2 to R 3 , where it is ultimately dissipated. However, due to the implicit ltering in LES the terms contributing energy toR 3 are not present and as a result are absent in the list above. Thus, the eect of the implicit lter in LES is to remove those terms contributing energy to R 3 . It is the terms contributing to R 3 and the action of terms exchanging energy with R 3 that is absent in LES. In DNS the tiny dissipative scales are resolved and these scales adjust to dissipate energy according to the Reynolds number. However, in under-resolved simulations the smallest scales are by denition not properly resolved and there is no adequate means for energy removal by these scales. There is eectively a discontinuity in the energy cascade in under-resolved simulations. In LES, the use of SGS models, such as the similarity model, which fail to account for the absence of the dissipative scales (e.g. by adequately removing energy near the mesh cuto), can lead to an accumulation of energy near the 58 cuto wavenumber. In the context of discrete wavenumber regions this energy tends to accumulate inR 2 . The accumulation of energy in LES due to inadequate energy removal (or dissipation) is depicted in the energy spectrum in Fig. 4.3. Figure 4.3: Turbulent energy spectrum with accumulation in the vicinity of the LES cuto. Given that the general ow of energy in turbulence is from the large energy-containing scales to the smallest scales and that the eect of the implicit ltering of the full nonlinear term in LES (i:e: u i u j ) is the removal of terms contributing energy transferred to the smallest physical scales, it is feasible that a SGS model can be formulated to properly account for the general role of energy transfer of the exact form of the SGS stress tensor, i.e. the removal of terms contributing energy to the smallest resolved scales in LES. This concept can be accomplished by removing specic interscale transfer terms contributing energy to R 2 in a manner that is analogous to the implicit lter which removes terms 59 contributing toR 3 . Alternatively, since in a fully-resolved simulation the actual subgrid- scale terms (unresolved in LES) act to remove energy from region R 2 , the proposed approach is to approximate this eect by preventing the energy from being deposited in R 2 by the action of the resolved energy transfer terms (shown in Eq. 4.28). This idea is the general concept behind the development of the current model. Specically, the goal is to formulate an alternative similarity-type SGS model that properly accounts for the natural process of interscale energy transfer in turbulence. 60 Chapter 5 Development of the Interscale Transfer Model 5.1 Basic Model Formulation The central goal of this research is to formulate a new SGS model that addresses the deciencies of the similarity model, specically this model's eect on SGS dissipation in the vicinity of the mesh cuto. The concept of interscale energy transfer oers a tool to parse the nonlinear transfer present in LES and to formulate expressions capable of removing (or reducing) terms that deposit energy near the LES cuto. The nonlinear term in the governing equations results in a purely advective transfer of energy between wavenumber modes that occurs due to the nonlinear interaction of eddies. In a turbulent ow there is no net energy change as a result of these interactions; however, if the turbulent energy spectrum is categorized into bands of wavenumbers, there can be a net exchange between these discrete regions. The decomposition of velocity into scales allows the interscale energy transfer terms to be categorized as production or redistribution terms based on the form of their corresponding partial nonlinear terms in 61 the momentum equation. The production of turbulent energy relies on the presence of a mean velocity gradient. Expressed with the additional restriction on the interacting scales dened by super- script m, the partial nonlinear terms have the form u q j @u p i @x j m : (5.1) Thus, if velocity component u p i present in these terms consists of the large-scale compo- nent of the LES velocity (i.e. u 1 i ), then this component along with the spatial derivative operator gives an approximation of the mean shear. In the energy equation terms of this form are taken to represent production type eects. In terms of this form,u q j identies the velocity component interacting with the large-scale shear. In these terms the lter level represented by superscript m identies the wavenumber region experiencing the eects of this interaction. Terms of this form with the additional restriction of m = 2 represent energy production eects in the smallest resolved scales. Conversely, terms containing the large-scale component of shear with the additional restriction of m = 1 are best de- scribed as advective terms. Despite an extraction of energy from the large-scale shear, these terms are chie y restricted to energy redistribution among scales in R 1 . In the above convention, the terms N 112 and N 122 are responsible for production eects in the smallest resolved scales of turbulence in LES. The term N 112 represents production in R 2 due to interactions among scales only in R 1 , while N 122 represents production in R 2 due to interactions involving scales in R 1 with those in R 2 . Both of these terms are expected to have an important role in the overall energy transfer from 62 R 1 to R 2 . A SGS model incorporating some (full or partial) combination of these terms can alter their role in the LES equations. A key aspect of selecting these terms is that by including both, the SGS model will have in uence over energy production in R 2 due to interactions from all resolved scales (i.e. R 1 and R 2 ). Specically, a SGS model that incorporates the negative of some approximation for the resolved stress components associated with termsN 112 andN 122 , has the capability of canceling some fraction of the energy production in R 2 caused by interactions from all resolved wavenumber modes. The complete removal of both the production and redistribution eects from a specic region risks freezing those scales and decoupling their activity from the cascade dynamics. Thus, it is the removal of only the production type terms that has been emphasized here. Incorporation ofN 212 was, however, considered for the sake of achieving model symmetry. Including this term in a model would not have the same interpretation as the inclusion of the production type terms discussed above; however, this term could still contribute to the model's SGS dissipation. Initial testing has shown that this form for a model, i.e. III ij = N 112 N 122 N 212 = d u i u j d u 0 i u 0 j u i u j + u 0 i u 0 j ; (5.2) can perform similarly to one removing only the production terms. However, the above form does not satisfy Galilean invariance and an acceptable surrogate to this model has not been found. For this reason, a model incorporating only the two key production terms in R 2 is being proposed. 63 In the eort to formulate a model capable of removing the eects of energy accumula- tion in the smallest resolved scales, reducing or removing terms T 112 andT 122 was of key interest. A SGS model formulated to remove the eects of the nonlinear energy transfer in R 2 associated with these terms takes the following general form II ij = N 112 N 122 = b u i b u j d b u i b u j b u i u 0 j d b u i u 0 j : (5.3) Combining terms, one obtains the simplied model expression II ij = d b u i u j b u i u j : (5.4) This is the primitive formulation of the interscale transfer model. A modication to this form of the model for the sake of Galilean invariance is proposed in the following section. 5.2 Galilean Invariance of the Modeled SGS Stress Galilean invariance is a property satised by the governing equations of uid dynamics. This property guarantees that the underlying physics in the governing equations remain unchanged under a transformation to an arbitrary inertial reference frame (Speziale [76]). Frame invariance is not, however, satised by the primitive form of the two-term interscale transfer model reported above. This violation can be demonstrated by comparing the 64 general form of the model to a form undergoing a Galilean transformation to an arbitrary inertial reference frame moving at xed velocity V k , specically ij = \ \ (u i +V i ) (u j +V j ) \ (u i +V i ) (u j +V j ): (5.5) Carrying out the ltering operations and simplifying one obtains ij = d b u i u j b u i u j +V i b u j u j +V j b b u i b u i : (5.6) Thus, it has been shown that the transformed expression for the model diers from the original by factors involving the reference velocity V k , and therefore violates frame invariance. Depending on a model's expression, Galilean invariance can be achieved by various means. In the present research, one such method is to remove the (implicitly ltered) mean velocity component from all velocity terms in the model, thus giving a model of the following form: ij = \ \ u i U i u j U j \ u i U i u j U j : (5.7) All V k terms will cancel in this expression and clearly this form of the two-term model would eliminate any extraneous terms introduced by a Galilean transformation. This model formulation was tested, however, and yielded relatively poor results when compared with results of the primitive form of the model. The likely reason for the poor performance is that the foundation of the proposed model relies on the presence of terms approximating 65 the mean velocity gradient. Removing the mean velocity component from these terms tends to undermine the objective of this approach. Another alternative form that also satises Galilean invariance is obtained by intro- ducing a double ltering operation in the N 112 term, specically N 112 = b b u i b b u j d b u i b u j : (5.8) When implemented using a physical lter, this modied form of theN 112 term more heav- ily weights the interscale transfer contributions toR 2 from the largest scales contained in R 1 . This modication is not expected to be necessarily detrimental to the intended role of the model. It can be shown that the above modications completes the Galilean invariant formu- lation of the model since the termN 122 satises invariance when aided by the divergence operator present in the governing equations. For incompressible ows, the imbalanced ltering of the velocity terms present in N 122 does not result in a violation of Galilean invariance because the divergence of the velocity in the momentum equation vanishes. This eectively removes any V k terms that remain. Specically, this is demonstrated below: 66 @ @x j N 122 = @ @x j b u i u 0 j d b u i u 0 j = @ @x j ( \ u i +V i ) (u j +V j ) ( \ u j +V j ) @ @x j \ ( \ u i +V i ) (u j +V j ) ( \ u j +V j ) = @ @x j b u i u 0 j d b u i u 0 j +V i @ @x j u 0 j b u 0 j = @ @x j b u i u 0 j d b u i u 0 j : (5.9) The divergence property of incompressible ows allows the expression for N 122 to remain in its original form. Adopting this term along with the the modied term N 112 and simplifying this expression, the modied Galilean invariant form of the model is obtained as II ij = N 112 N 122 = d b u i u j b b u i b b u j b u i u 0 j : (5.10) In testing, the double ltering applied in the modied form of the model had only a minor impact on the model's predictions. Thus, it is the above formulation, represented as ( II ij ) and referred to as the "interscale transfer model", that is the adopted form of the proposed model and has been implemented in the present research. 67 5.3 Interscale Energy Transfer of the Similarity Model Given that the similarity model is also constructed from test-ltered quantities, it too can be considered in terms of its role in interscale energy transfer. As outlined by Domaradzki [19], adopting a model constant of unity, the decomposition of this model into its partial components of stress gives d u i u j b u i b u j = N 111 +N 121 +N 211 +N 221 N 111 +N 112 : (5.11) The impact of this model is the subtraction of N 112 and the addition of terms N 121 ; N 211 , andN 221 to the nonlinear term on the left side of the LES equations. Thus, the similarity model is related to the present model in that both forms remove N 112 , and so these models remove the eects of the energy production in R 2 due to its presence. This aspect of the similarity model is expected to be benecial with respect to the model's ability to account for the removal of energy near the LES cuto. However, the similarity model also has the eect of duplicating three of the terms representing transfer in R 1 . Clearly, the duplication of these terms will alter the character of the nonlinear interactions beyond simply aecting energy production in the smallest scales. It is this modication to the dynamics of interscale energy transfer that is in contrast to the key role of the exact form of the subgrid-scale stress tensor in fully resolved turbulence, specically its removal of terms contributing to the small scales. In terms of energy transfer there seems to be little physical rationale for the duplication of these terms. Furthermore, some of the benet of removingN 112 is likely negated by the duplication of those terms aecting energy transfer in R 1 . It is expected that it is this feature that has earned the category 68 of similarity models the general reputation of providing inadequate energy dissipation. Through testing, however, it will be shown that the similarity model may produce con- siderable levels of SGS dissipation, although in contrast to a SGS model's expected role, the terms accounting for the majority of this are those responsible for energy transfer in the largest resolved scales. The present interscale transfer model has been proposed as an alternative to the scale- similarity model. This presently proposed model and the classic form of the similarity model have been implemented in a numerical solver and the results of this testing are compared in this research. 5.4 Test Filtering in LES By introducing a SGS model that removes partial nonlinear terms responsible for energy production in the smallest scales of LES, it emulates the role of the exact form of the SGS stress tensor by encouraging the most prominent exchange to occur from the largest to the smallest scales of turbulence. The eects of nonlinear production fromT 112 andT 122 can be removed by constructing their resolved stress components and incorporating them into a SGS model. In spectral space, discrete wavenumber bands can be selected through corresponding Fourier coecients giving the separate wavenumber regions. Construction of these terms with spectral ltering gives sharp divisions between the high and low wavenumber regions. However, spectral ltering is not convenient for many types of CFD codes. On the other hand, ltering in physical space is more amenable to a broad range of LES solvers. Physical ltering, however, does not typically give a sharp division 69 between wavenumbers and the denitive nature of the wavenumber bands becomes less apparent. Intuitively it may seem ideal then to use a lter with a sharp cuto to limit removal only to the smallest scales. However, using a smooth lter allows for some energy removal at larger scales and relaxes the requirement that energy rst be passed to the smallest scales before action is taken. In practice, the ltering of velocity in physical space is typically carried out through the application of a lter kernel G, which weights neighboring values on a mesh. Again, this operation is expressed as a convolution integral of the form b u i (~ x) = Z G ~ x~ x 0 ; u i ~ x 0 d~ x 0 : (5.12) The operation generates a spatially-weighted average of the velocity u i (x), where G sets the weighting applied at physical points on the mesh. In LES, a lter function is typically selected that satises the normalization constraint Z +1 1 G ~ x ~ x 0 ; d ~ x 0 = 1: (5.13) In physical terms, this expression states that the full eect of the weighting applied by the lter function sums to unity. Each of the lters commonly used in LES have specic characteristics and a unique convolution kernel. For illustration, the sharp spectral lter satisfying the above normal- ization constraint has the following lter kernel in physical space 70 G ~ x ~ x 0 ; = sin j~ x ~ x 0 j ~ x ~ x 0 : (5.14) The Gaussian lter has a lter shape corresponding to a symmetric distribution function with zero mean. The Gaussian lter in physical space has the form G ~ x ~ x 0 ; = 6 2 1=2 exp 0 B @ 6 ~ x ~ x 0 2 2 1 C A: (5.15) Another lter commonly used in LES is the box lter. The box ltering operation in physical space is represented as G ~ x ~ x 0 ; = 8 > > < > > : 1= 0 if ~ x ~ x 0 < =2 otherwise : (5.16) Work by Stolz et al. [77] has shown that physical approximations of a large variety of lter shapes are possible by applying adaptations to the discrete implementation of the box lter. These n-point ltering techniques oer greater control over the lter shape. The above expressions represent a few of the spatial lters encountered in LES. Typ- ically a lter is chosen for its shape and cuto characteristics. A discrete representation of a ltering operation can be obtained by numerically approximating the integration. For instance, a three-point lter kernel can be formulated as an approximation of the box lter. Two common approximations of the box lter are obtained by applying Simpson's 1/3 rule for integration and the trapezoidal rule. The discrete smoothing kernel that results corresponds to a spatial weighting of a function f(x) at neighboring values on a mesh. 71 Recalling that Simpson's 1/3 rule for integration states b Z a f (x)dx ba 6 f (a) + 4f a +b 2 +f (b) : (5.17) Therefore, if ltering is applied to a function over an interval corresponding to two com- putational cell widths, Simpson's rule for integration gives an operation of the form 1 2x +x Z x f x 0 dx 0 (2x) (12x) [f x + 4f 0 +f +x ]: (5.18) This approximation of the integral is equivalent to multiplying the length of the integra- tion interval (ba) by a function b f(x), which corresponds to f(x) with some smoothing kernel applied. Applying the test ltering operation at discrete points on an LES mesh gives a lter kernel 1 6 ; 2 3 ; 1 6 , and a discrete ltering operation of the form b f (x n ) = 1 6 f x n1 + 2 3 f (x n ) + 1 6 f x n+1 : (5.19) Similarly, the trapezoidal rule for integration is dened as b Z a f (x)dx ba 2 [f (a) +f (b)]: (5.20) 72 If ltering is again taken over an interval of two cell widths, the trapezoidal rule is applied in a similar manner. Choosing again an expression for G corresponding to the box lter, one obtains 1 2x +x Z x f x 0 dx 0 = 1 2x 2 4 0 Z x f x 0 dx 0 + +x Z 0 f x 0 dx 0 3 5 1 2x x 2 [f x + 2f 0 +f +x ]: (5.21) Thus, an approximation of the box lter using the trapezoidal rule results in a lter kernel of 1 4 ; 1 2 ; 1 4 . For simplicity, the ltering operations have been shown in one-dimension; however, in LES, ltering is often applied in two or more directions on a Cartesian mesh. Multidi- mensional ltering is carried out by applying the lter kernel sequentially in each spatial direction. The eect of ltering is often characterized by considering its eect on a function of the form e ikx , where k represents wavenumber. The ltering operations represented by the transfer functions are often plotted versus normalized wavenumber showing the lter's amplitude response, or specically its attenuation eects across the resolved range of wavenumbers. Phase response is not typically a concern since the ltering is applied at discrete time steps with a spatially symmetric lter kernel. Spatial symmetry may not be guaranteed for ltering on a non-isotropic mesh, although this eect can be reduced through interpolation. Applying the lter kernel derived from the trapezoidal rule to the test function gives the expression 73 d e ikx j =e ikx j 1 4 e ikx + 1 2 + 1 4 e ikx : (5.22) The term in brackets is the lter transfer function and it represents the eect of the smoothing kernel applied to the functione ikx . Making use of trigonometric identities, the transfer functions corresponding to the discrete approximations of the box lter using the trapezoidal rule G T and Simpson's 1/3 rule G S are expressed as: G T (k) = 1 2 cos (kx) + 1 2 ; (5.23) G S (k) = 1 3 cos (kx) + 2 3 : (5.24) The transfer functions induced by these ltering operations are expressed above in wavenumber space. The eect of the transfer functions are often considered up to a maximum wavenumber of k max = 2= (2x), where 2x corresponds to the smallest wavelength that can be resolved by the mesh. The attenuation eects of a transfer func- tion can be plotted versus normalized wavenumber, where wavenumber k is normalized by 1=x. Thus, the normalized wavenumber is considered over the range 0. Typical test lters for LES remove content from the high wavenumber modes present in the resolved eld. As a result, the test-ltered LES velocity b u i is generally smoother than its unltered counterpart u i . As an illustration of the eect of the lter, the am- plitude response of the lter's transfer function is plotted versus normalized wavenumber 74 to identify the eect of the lter across the wavenumber spectrum. Comparisons of the amplitude response for the two approximations of the box lter are shown in Fig. 5.1. Figure 5.1: Amplitude responses of box lter approximations: Trapezoidal rule 1 4 ; 1 2 ; 1 4 , Simpson's rule 1 6 ; 2 3 ; 1 6 These plots of amplitude response illustrate how varying the test lter coecients can control the cuto and lter shape. The behavior of the two representations of the box lter, particularly at high wavenumbers, are clearly dierent. Specically, the cuto of the lter obtained from Simpson's rule does not fully remove scales at the highest wavenumber, whereas, the trapezoidal lter does. The dierence in this behavior is controlled by the simple relation between ltering coecients (a 1 ;a 0 ;a 1 ), specically, the value ofa 0 relative toa 1 . In general, increasinga 0 relative toa 1 results in a weaker lter, 75 although limits to this relation exist. The dierent weighting of the coecients translates to a change in the amplitude response given by the representation of the lter. There are limitations to the lter shapes that can be obtained from the three-point kernel. Physical lters with sharper cutos can, however, be formulated using straight- forward adaptations of the box lter implementation. These techniques oer considerable control over the lter shape. However, the resulting lter kernels must operate on more than three points and some approaches require solving a system of equations (Stolz et al. [77], Jeanmart and Winckelmans [36]). The ltering procedure outlined by Jeanmart and Winckelmans is implemented in a manner similar to the three-point implementation of physical lters; however, in this case the function operating on velocity is altered as follows: b u i (~ x) = Z [I [IG] n ]u i ~ x 0 d ~ x 0 : (5.25) In this expressionI represents the identity operator and the operatorG can take the form of some lter function corresponding to a box lter, Gaussian lter, etc. The exponent n is an integer value that, along withG, controls the lter shape. ChoosingG to represent a discrete approximation of the box lter using the trapezoidal rule, the amplitude response for various values of n are shown in Fig. 5.2. Increasing the exponent n results in a steepening of the cuto and a broadening of the eective lter kernel. This implementation of ltering in physical space is capable of giving a better representation of the sharp cuto associated with spectral ltering. 76 Figure 5.2: Amplitude responses of high-order lters: n = 1; 2; 8 Of course limitations in this technique are reached when the ltering stencil reaches the extent of the spatial resolution in the domain. The eect of applying a sharper lter in the proposed SGS model is that more well- dened wavenumber regions are represented for the sake of interscale energy transfer considerations. However, the above ltering implementation was tested with the present model and was not found to improve upon results obtained using simple three-point kernels in physical space. Instead, it was found that a more smooth lter with a relatively weak cuto at the highest wavenumbers was necessary to obtain good results for the present model. As a result, a three-point lter resembling the discrete approximations of the box lter in physical space has been applied in the current work. 77 The test lter shape and the control of the cuto level are clearly important in inter- preting the role of the interscale energy transfer terms. Test ltering applied in physical space can generally give a smooth lter shape; however, the divisions between wavenum- ber regions are typically not as well dened and so the resulting partial nonlinear terms are only approximate representations of those terms idealized in spectral space. In testing, the smoothing characteristics of physical ltering have actually shown benecial eects in the present model. In addition to selecting the type of test lter, the lter's eect on the partial nonlinear terms was taken a step further in the model's initial implementation. The test lter's kernel was modied initially to tune the lter strength, thus dening the eect of the terms across the energy spectrum. The lter strength was adjusted by changing the shape of the lter function G applied in the convolution, i.e. Eq. (5.12). Specically, the lter strength was controlled using the relative magnitudes of the coecients in the three-point kernel (i.e. a 0 versus a 1 ). At a given wavenumber, this adjustment controls the fraction of the partial nonlinear terms removed by the model. The optimum setting of the lter kernel was found to be somewhat weaker that what is given by the standard approximations of the box lter. This represents a reduced eect of the terms removed by the model. The partial removal of the energy transfer terms conceptualized in spectral space means that some fraction of these terms must remain in the equations. Results suggests that this feature is actually helpful in maintaining the role of the smallest scales in the ow. Specically, a smooth lter likely helps to avoid decoupling those scales in the largest wavenumbers by eliminating the strict wavenumber dependence of the energy removal mechanism and by relaxing the requirement that modal 78 activity be passed across a xed cut-o wavenumber before it can be removed. Thus, the role of the smooth lter shape implemented with the terms adopted in the present model gives an energy removal mechanism that simply becomes increasingly more active at higher wavenumbers. 5.5 Filtering on a Non-uniform Mesh In the present LES the test-ltering was carried out in all coordinate directions. In all cases a uniform mesh spacing was used with the exception of the wall-normal direction where no-slip boundary conditions were applied and the mesh was rened near the walls in order to better resolve the local velocity gradients and the near-wall structures of turbu- lence. Quadratic interpolation was implemented for velocity ltering in the wall-normal direction in order to reduce integration errors that would arise from equally weighting neighboring points that have dierent spatial separations. This interpolation was carried out using the approach outlined by Loh and Domaradzki [49]. The modied expression for the test lter operation distinguishes between the sizes of adjacent cells ( and + ) and takes the form b f (x n ) =a 1 f x n1 +a 0 f (x n ) +a 1 f x n+1 ; (5.26) where the lter kernel is represented as (a 1 ;a 0 ;a 1 ), and f (x) represents a dependent variable modied as follows to account for the unequal mesh spacing: 79 f x n+1 = + + + f x n1 + 2 1 + f (x n ) +2 + ( + + ) f x n+1 : (5.27) This modication is applied to those points having the larger of the two neighboring cell spacings in the three-point stencil. The correction acts to minimize the phase-shift associated with ltering using a uniform lter kernel on a nonuniform mesh. The result is a three-point lter kernel that is adjusted for mesh nonuniformity and is applied up to the wall when ltering in the vertical direction. This modied ltering formulation was used for all test ltering carried out in the present eort. 80 Chapter 6 Model Testing All numerical simulations were conducted using the pseudo-spectral ow solver developed by Diamessis et al. [18], and modied and implemented by Tantikul and Domaradzki [78], [79]. The details of the numerical solver are summarized in section 6.2. The proposed SGS model has been implemented in the ow solver and has been used to simulate channel ows. The similarity model and the dynamic model have also been implemented and these results have been reported for comparison. The predictions from these simulations are compared against DNS carried out in the present research and also DNS results provided by by del Alamo and Jimenez [15], [16] and del Alamo et al. [17]. To provide a background, some basic concepts along with a brief summary of the theory associated with wall-bounded turbulent ows are presented in the following section. 6.1 Channel Flow Fundamentals Shear ows make up an important class of turbulent ows due to their presence in various applications. In channel ow and in many other types of ows the mean shear provides the energy for turbulent motion. The bulk characteristics of free shear ows evolve 81 independent of viscosity, while wall-bounded shear ows see the important eects of viscosity at locations close to the wall (Kundu [41]). Physical considerations have been used to arrive at expressions for the mean velocity prole in wall-bounded ows. The so-called \law of the wall" expresses the velocity as a functional relationship with distance from the wall. The velocity in the thin region adja- cent to the wall, referred to as the viscous sublayer, is expressed as a linear relationship with distance. Farther away from the wall a logarithmic relation between velocity and distance has been demonstrated. Thus, a typical boundary layer can be divided into an inner viscous sublayer, an outer logarithmic layer, and a region referred to as the buer layer that serves as a transition region between these two. When expressed in terms of ap- propriate non-dimensional variables, plots of velocity proles versus wall-normal distance fall along the same curve for ows at varying Reynolds numbers. Writing an expression for the shear stress at the wall one obtains wall dU dz z=0 : (6.1) Combining this expression with density gives a quantity u having dimensions matching that of velocity, and is expressed as u r wall : (6.2) This expression is referred to as the friction velocity. In channel ows, the friction Reynolds number formulated from friction velocity, chan- nel half-width (), and kinematic viscosity, has the form 82 Re u : (6.3) By combining the friction velocity with viscosity and a length scale measuring the physical distance from the wall, one obtains a dimensionless group, which gives a descrip- tion of separation distance from the wall, i.e. z + zu : (6.4) By combining the friction velocity with the mean velocity (in a pipe, channel, etc.), another important dimensionless quantity is obtained, i.e. u + U u : (6.5) The dimensionless groupsu + andz + have been found to share a relationship of a general form that describes how velocity varies in a turbulent boundary layer. This general expression is given as u + =f z + : (6.6) Near the wall, the dimensionless velocity has been found to have the following rela- tionship: u + =z + : (6.7) 83 This is the linear relationship in the viscous sublayer and it has been found to hold true up to values of around z + 7. Farther from the wall, beginning at z + 30, the velocity prole takes on the relationship which behaves logarithmically with the non-dimensional distance. The velocity-distribution law for large Reynolds numbers (Schlichting [70]) takes the form u + = 2:5ln z + + 5:5: (6.8) Clearly the velocity proles in the viscous sublayer and in the logarithmic layer exhibit dierent behaviors; however, the velocity versus distance relations in each region merely assume dierent functional forms of the same dimensionless variables. The expressions for wall separation distance and physical velocity in each region can be scaled identically for various Reynolds numbers. Numerous experimental measurements have validated these relationships. In the outer region of the boundary layer the logarithmic prole eventually breaks down and additional theories have been developed to describe the velocity in this region. The extent that the logarithmic prole is valid depends on the Reynolds number. A higher Reynolds number extends the non-dimensional thickness of a fully developed turbulent boundary layer, thus giving a thicker log-layer. In LES a model's ability to correctly predict statistical quantities of both the mean and uctuating components of velocity is critical. A calculation of turbulence intensity is often used to quantify the level of turbulence in a ow-eld. This quantity is expressed as the following function of the uctuating component of velocity: 84 1 3 D (u i U i ) 2 E : (6.9) This is eectively an averaged quantity over all components of velocity. The turbulence corresponding to the individual velocity components normalized by the friction velocity in the viscous wall region can be expressed similarly in the form of a standard deviation (Mathieu [52]) * u i U i u 2 + 1=2 : (6.10) The quantities obtained from this expression are referred to as the root-mean-square of the velocity uctuations. These RMS velocities provide a statistical representation of the turbulence in each coordinate direction. The study of the predicted distributions of RMS velocities in shear ows is useful because they re ect the combined eects of production, diusion, and dissipation of turbulence in the ow. These quantities help describe how the turbulence evolves and how it distributes itself once it has been generated. Analytical and experimental predictions of wall-bounded ows have been used exten- sively for anchoring numerical ow solvers. Wall-bounded ows have also been used to anchor and test turbulence models and have been a useful testbed for SGS model devel- opment in LES. The present SGS model has been implemented in channel ow LES. Its ability to predict mean and RMS velocity proles at various resolutions and Reynolds numbers has been evaluated. The results of such tests are presented in following sections. 85 6.2 Numerical Solver An incompressible multidomain pseudo-spectral solver has been used to evaluate the performance of the proposed SGS model. Spectral solvers are not typically considered ideal for general CFD applications due to their limitations in eciently treating problems with complex geometries or non-traditional boundary conditions. However, spectral codes are ideal for SGS model testing and validation because they oer the highest level of numerical accuracy relative to other techniques (Fletcher [30]). This high level of accuracy does not guarantee good ow predictions by the code, but rather it allows the user to more directly associate errors in the simulations with the performance of the model rather than numerical errors generated by the solver. The incompressible ow solver developed by Diamessis et al. [18] and modied by Tan- tikul and Domaradzki [78], [79], was used in the present eort. This code has been used to simulate wall-bounded ows at various Reynolds numbers. The solver uses a multidomain implementation of the spectral method. Temporal discretization of the governing equa- tions is achieved using a high-accuracy pressure projection scheme. The advancement of the velocity solution is decomposed into three fractional steps. The computational domain is divided into multiple subdomains extending spanwise across the channel. In this approach a patching condition along with a penalty method is used to treat values across the subdomain interfaces. For reference, the resolution for all cases is reported in terms of its global discretization, i.e., N x , N y , N z . Fourier spectral discretization is used withN x andN y modes specied in the streamwise and spanwise directions, respectively. 86 In the vertical direction spectral discretization is implemented using Legendre polynomi- als dened on N z points, the order of the polynomial approximation being (N z 1) in the vertical direction. Point-wise interpolation is implemented in nodal form (i.e. with collocation points) using Lagrange polynomial shape functions. The domain used in the simulations was constructed from a Cartesian grid where the streamwise direction is oriented along the x-axis. No-slip boundary conditions were applied at the vertical boundaries, while periodic boundary conditions were applied at the streamwise and spanwise borders of the domain. The domain represents a channel bounded vertically by solid surfaces and having innite width and length. The geometry is illustrated in Fig. 6.1. Figure 6.1: Channel geometry. Three dierent channel sizes were used in the simulations. The geometry chosen for the simulations corresponds to cases which have been analyzed extensively in previous studies, or for which sucient data was available for comparison. In all cases the selected Reynolds number (Re ) closely coincides with that for which DNS results were available. 87 A rough guideline for choosing mesh resolution in LES is to resolve at least 80% of the kinetic energy contained in the turbulence (Pope [65]). This requires ne near-wall mesh resolution in order to capture the near-wall structures of turbulence. More quantita- tive guidelines for selecting mesh resolution have been established for wall-bounded ow simulations (Sagaut [69], Piomelli [63]). Assuming the mesh resolution is uniform in the streamwise and spanwise directions and cell clustering is applied near the walls, the guidelines for mesh spacing can be expressed spatially in wall units as x + i = x i u : (6.11) By entering the physical cell size in each coordinate direction into the formula the mesh resolution in terms of wall units is obtained. In the wall-normal direction where mesh spacing varies, the size of the rst cell adjacent to the wall was used. Uniform spacing was used in the streamwise and spanwise directions. Given that the x, y, and z directions correspond to the streamwise, spanwise, and wall- normal directions, respectively, the rough quantitative guidelines for mesh resolution are given as x + 40 60; y + 20 30; z + 1: (6.12) 88 These guidelines help ensure that the near-wall structures of turbulence and the energy- containing scales of the ow are adequately resolved and thus reasonable ow predictions can be achieved when using typical LES models such as the Smagorinsky model or the dynamic model. However, regardless of these guidelines, the results of such simulations of course depend on the predictive quality of the SGS model. A poorly formulated SGS model will generally give poor results regardless of resolution. 6.3 Test Filter Optimization The eects of nonlinear production from T 122 and T 112 can conceptually be removed by constructing the appropriate partial components of stress using physical ltering and incorporating them into a SGS model. Interscale energy transfer conceptualized with sharp spectral ltering can only be approximated using ltering in physical space. In the current work, the smooth aspect of physical ltering has proved benecial in the sense that the full representation of the partial nonlinear terms would result in the complete removal of the eects of selected energy transfer terms. This full cancelation of terms is not likely desired, in fact testing has suggested the model's performance currently relies on only a partial cancelation of those terms incorporated into the model. The implementation of a naturally smooth test lter accomplishes this. Thus, the current research introduces a SGS model with a physical lter that results in what can be interpreted as the partial removal of specic energy transfer terms idealized in spectral space. The lter used in testing the model was a three-point stencil with quadratic inter- polation used to account for non-isotropy of the mesh. Initially this testing was carried 89 out in wall-bounded turbulence simulations using the standard three-point trapezoidal approximation of the box lter in physical space. The initial results showed that only marginal predictions of the mean velocity prole could be obtained. However, in testing the model the quality of these predictions exhibited a noted dependence on lter strength. Testing showed that the predictions could be improved by modifying the strength of the test lter. Specically, a weaker lter was implemented and the model's predictions of mean velocity were improved considerably. Adjustments to the lter strength were made by modifying the lter kernel (a 1 ;a 0 ;a 1 ) where the sum of the coecients satises the normalization constraint and a higher or lower weighting of the middle coecient (within certain limits) generally represents a weaker or stronger lter, respectively. An optimum lter setting was established by running a low Reynolds number case (Re = 180) and comparing with predictions from DNS results. The lter setting was optimized for this test case using a trial-and-error modication of the lter coecients. Comparing the model's predictions with the ideal mean velocity prole, the optimum test lter for the Galilean invariant version of the two term model was found to be 1 12 ; 10 12 ; 1 12 . After optimizing the lter setting additional Reynolds numbers were simulated and no additional changes to the lter were needed. The proposed model has continued to give good results in subsequent testing. Comparisons of the amplitude response for the selected lter versus two standard approximations of the box lter: 1 4 ; 1 2 ; 1 4 and 1 6 ; 2 3 ; 1 6 , along with a much weaker lter 1 32 ; 30 32 ; 1 32 , are shown in the Fig 6.2. It should be noted here that the lter coecients 90 were selected solely for the sake of modifying the lter strength and this was done with- out consideration of the predictive quality of the lter kernel for general applications in numerical integration, etc. Figure 6.2: Test lter amplitude responses. Filter kernels: 1 4 ; 2 4 ; 1 4 , 1 6 ; 4 6 ; 1 6 , 1 12 ; 10 12 ; 1 12 ; 1 32 ; 30 32 ; 1 32 . In LES, changes in the test lter strength correspond to changes in the relative sizes of the resolved wavenumber regions (i.e. R 1 versusR 2 ). In the interscale transfer model, changes to the test lter strength correspond to changes in the fraction of the energy pro- duction terms removed by the model. The present model's dependence on lter strength is analogous to the dependence of dissipative models on the eddy-viscosity term. This relationship can be demonstrated in the Smagorinsky model by altering the value of the 91 Smagorinsky coecient c s in channel ow simulations. A greater value of the Smagorin- sky coecient represents an increase in the model's eddy-viscosity, which ultimately alters the diusive transport of momentum across the boundary layer. A similar relation exists in the interscale transfer model. In Fig. 6.3 this relation is displayed by altering the lter strength in the present model. In the process of optimizing the proposed model's lter strength, the lter coe- cients were selected which gave the best results for predictions of mean velocity for a low Reynolds number case (Re = 180) at a relatively low mesh resolution (N x ;N y ;N z ) = (32; 32; 33). To illustrate the model's dependence on lter selection, mean velocity re- sults for wall-bounded ow with this Reynolds number and resolution are shown for the four discrete lter kernels discussed above. Also shown are results for the Smagorinsky model (implemented with wall damping [57]) with diering values of c s . Approaching the lower extremes of the Smagorinsky constant (c s =0) and the lter strength (0; 1; 0), the results simply approach under-resolved DNS (i.e. no SGS model). Near the upper extremes of the Smagorinsky coecient and lter strength, both of these models have the eect of overly damping the turbulence in the ow. For comparison, the ideal log-linear velocity prole for ow in a turbulent boundary layer is also plotted. The ideal law of the wall type prole that is shown in the gure (labeled as \Linear/Log- Law") closely matches results obtained from the DNS conducted in the present eort using the identical domain size and ow conditions. Clearly the present model exhibits a strong sensitivity to the lter strength. It is clear that for the Galilean invariant version of the proposed interscale transfer model, selection of the lter kernel 1 12 ; 10 12 ; 1 12 gives a result that best matches the actual velocity prole in the turbulent boundary layer. 92 Figure 6.3: Model test case: Re =180, Mesh: (32, 32, 33). Log-Linear Law (dotted); Present model (solid) with lter kernels: (a 1 ;a 0 ;a 1 )= 1 4 ; 1 2 ; 1 4 , 1 6 ; 2 3 ; 1 6 , 1 12 ; 10 12 ; 1 12 , 1 32 ; 30 32 ; 1 32 ; Smagorinsky model (dashed) with constants: c s =0.02, 0.06, 0.1, 0.13. Mean velocity increases with both c s and a 1 . Filter kernels 1 4 ; 2 4 ; 1 4 , and 1 6 ; 4 6 ; 1 6 result in an upward shift of the velocity prole, specically in the logarithmic layer. The weakest lter with kernel 1 32 ; 30 32 ; 1 32 , gives a result that considerably underpredicts the theoretical mean velocity prole in this region. Since the partial nonlinear terms in the model remove terms aecting energy pro- duction in R 2 , an increase in lter strength results in an increase in the size of R 2 and thus a general increase in the energy transfer terms removed by the model. Clearly the model's eect on the small scales is analogous to the intended role of the dissipative eects 93 of eddy-viscosity models, specically the suppression of activity in the smallest resolved ow structures. The well known drawback of the Smagorinsky model is the need to change the model's coecient depending on Reynolds number or ow type. Testing at dierent Reynolds numbers, resolutions, and domain sizes has shown that the present model does not demon- strate this same dependence on lter strength. The same lter kernel has been used for all results presented for the interscale transfer model. The present testing, however, has been limited to wall-bounded ow simulations. It should be noted that adjustment of the test lter in the similarity model does not yield the same behavior as that exhibited above by the interscale transfer model. This is most likely because the similarity model has a mixed interpretation of its energy transfer eects given this model's adopted nonlinear terms. 6.4 Case Descriptions In LES, the SGS model directly impacts the exchange and dissipation of turbulent energy in the ow. The nature of this dissipation aects the turbulent velocities and ultimately the mean velocity of the ow. A fundamental feature of a SGS model is its ability to accurately predict mean and RMS velocities at various ow conditions. Predictions of these statistical quantities were emphasized in testing of the present model. Since the ow dynamics can change considerably with test conditions it is useful to demonstrate the model's predictability over a sizeable test space. The ow solver was used to simulate several cases of wall-bounded turbulence. Plots of mean and RMS velocities for various 94 domain sizes, resolutions, and Reynolds numbers are presented. The present model along with the similarity model and the dynamic model were implemented and these results have been included for comparison. The predictions from these simulations are compared against DNS carried out in the present eort and DNS results provided by del Alamo and Jimenez [15], [16], and del Alamo et al. [17]. The domains used in the simulations were discretized using a Cartesian mesh where the streamwise direction is oriented along the x-axis. Uniform spacing was used in the streamwise and spanwise (y) directions while cell-clustering was used in the wall-normal (z) direction. The domain represents a rectangular channel bounded vertically by solid surfaces and having innite width and length. Wall boundary conditions were applied at the vertical boundaries, while periodic boundary conditions were applied at the stream- wise and spanwise borders of the domain. In industrial CFD applications, the geometry is often very complex and strict control over mesh resolution is not always possible. Thus, it is worthwhile to evaluate a SGS model's robustness by testing its performance outside the range of recommended mesh resolution. In choosing the resolution for each case there were two objectives; the primary objective was to evaluate the performance of the model using the resolution guidelines, and the second was to evaluate the model's robustness by running simulations that de- viated from these guidelines. In setting the resolution for all cases, equal partitioning across parallel processors was required by the solver. The mesh resolution was selected according to this requirement and the aforementioned objectives. Table 6.1 summarizes the LES cases considered and shows the Reynolds number, channel dimensions, global discretization, and mesh resolution expressed in wall units. 95 The vertical resolution in wall units is shown for the rst cell height adjacent to the wall. Also shown is the ratio of LES to DNS resolution normalized for domain size. These ratios underscore the fact that the computational expense of LES is typically a very small fraction of DNS. Re L x ;L y ;L z N x ;N y ;N z x + ; y + ; z + N LES =N DNS 180 4; 2; 2 32; 32; 33 71; 35; 1:2 0.008 180 4; 2; 2 64; 64; 65 35; 18; 0:32 0.063 950 8; 3; 2 512; 384; 85 47; 22; 0:99 0.0061 1050 5=2;=2; 2 96; 128; 65 86; 13; 1:8 0.0056 1050 5=2;=2; 2 128; 128; 97 64; 13; 0:83 0.011 2000 5=2;=2; 2 256; 96; 131 61; 33; 0:86 0.0034 Table 6.1: LES case descriptions: Reynolds number, domain dimensions, mesh resolution, wall-unit resolution, ratio of LES to DNS mesh resolution normalized for domain size. Generally, the domain sizes used in LES were smaller that those used in the DNS; however, the exceptions to this are theRe =180 domain, which matches the DNS domain used in the present work, and theRe =950 case, which implemented the identical domain size used in the DNS supplied by del Alamo et al. [16]. Note that since DNS results for Re =1050 were not available, the results atRe =950 have been used for comparison with these cases. Details of all the DNS cases are given in Table 6.2. Re L x ;L y ;L z N x ;N y ;N z x + ; y + ; z + 180* 4; 2; 2 256; 256; 131 8:8; 4:4; 0:077 180** 12; 4; 2 762; 342; 97 8:9; 6:1; 0:096 950** 8; 3; 2 3072; 2304; 385 7:8; 3:9; 0:032 2000** 8; 3; 2 6144; 4608; 633 8:2; 4:1; 0:32 Table 6.2: DNS case descriptions: Reynolds number, domain dimensions, mesh resolution, wall-unit resolution. *Results provided by the author. **Results provided by del Alamo and Jimenez [15], [16], and del Alamo et al. [17]. 96 The DNS conducted in the present eort matches the low Reynolds number DNS supplied by del Alamo and Jimenez [15]. These two simulations dier in their domain size but only slightly in their eective resolution. The DNS conducted in the present eort uses channel dimensions identically matching those used in the LES, while the domain used by del Alamo and Jimenez [15] was considerably larger. For this reason, the present DNS is expected to give a more direct comparison with the LES. The two DNS cases gave quite similar results, the only notable deviation was in their prediction of spanwise RMS velocity. The degree of deviation, while not unusual, is possibly explained by the dierences in spanwise resolution. The present DNS used slightly higher spanwise resolution perhaps better capturing the ne-scale structures of turbulence accounting for the spanwise uctuations. Hence the slight increase in spanwise activity predicted by the present results. Dierences in the DNS channel dimensions also have the capability of impacting the RMS activity. Regardless of the cause, this dierence in the DNS results was acceptably small and in general the two DNS cases provided ow statistics that were quite similar. In the current work, all LES with a SGS model enabled were started from an initial condition obtained by running the code to statistically steady-state with no SGS model. The data used for initialization is best described as under-resolved DNS (URDNS) since all settings used in these simulations match those of the DNS with the exception of the ne mesh resolution. In all cases, the resolution used for the initialization run was identical to that used in the LES. Comparing the under-resolved DNS results to those obtained with a SGS model provides a method to isolate the eects of the SGS model in the solution. For reference, DNS results are included in all gures. Unless otherwise stated, the DNS 97 results shown in the gures are those provided by del Alamo and Jimenez [15], [16], or del Alamo et al. [17]. All LES cases were advanced in time until a statistically steady-state solution was reached. Steady-state conditions were veried by plotting resolved and turbulence shear stresses and other statistical quantities of the ow. Statistical samples were collected by averaging over homogeneous planes. The planar averages were computed over many realizations. The simulation washout time scales are estimated by computing the ratio of the streamwise channel dimension and the bulk ow velocity (see del Alamo and Jimenez [15]). The LES cases corresponding to Re =180, 950, 1050, 2000 were computed with the model of interest enabled for approximately 25; 30; 110, and 25 washout time scales, respectively, prior to collecting data. The data was collected over another 30; 10; 60, and 15 washout times, respectively. In these cases the data collection involved averaging over 60-240 ow realizations. 6.5 Modeling Results A mean pressure gradient between the in ow and out ow boundaries was applied in the simulations. The presence of this gradient gives rise to a corresponding rate of mean ow and a generation of shear stress at the walls. In this study, the interest is in evaluating characteristics of the mean velocity and the RMS quantities of turbulent uctuations after the eects of the friction forces exerted at the walls have diused to the core of the channel and the ow has reached a fully-developed state of turbulence. The non-dimensional mean velocity proles are plotted on semi-logarithmic scales for all cases shown in Table 6.1. 98 With the exception of the DNS data supplied by del Alamo and Jimenez [15], [16], and del Alamo et al. [17], all results presented in this section were obtained during the current research eort. Included in each plot of mean velocity are DNS results, which generally reproduce the theoretical linear/logarithmic prole expected of turbulent ows. Results obtained from implementing the present model are shown, and for comparison additional LES results are included. Specically, results from URDNS, the dynamic model, the similarity, and in some cases the Smagorinsky model are also shown. In plots of RMS velocity, the DNS data is again plotted along with the LES results. 6.5.1 Mean Velocity Results Channel ow simulations have been a mainstay in turbulence model development because the characteristics present in these ows lend themselves to many other types of internal and external ows. The logarithmic velocity distribution is a fundamental property of wall-bounded turbulent ows. Thus, an important feature of a SGS model in LES is its ability to predict this common prole for various ows. Since the dynamics of turbulence can change considerably with ow conditions it is informative to demonstrate the model's predictability at various Reynolds numbers. TheRe =180 case has a physical Reynolds number based on mean centerline velocity (Re c ) and the channel half width of approximately 3,300. This Reynolds number, which has a relatively a coarse resolution (32, 32, 33), served as the initial test case for which the lter kernel implemented in the present model 1 12 ; 10 12 ; 1 12 was established. The results obtained by the proposed interscale transfer model (denoted in these results as \model") are naturally good as a result. The present model performed reasonably well in both the 99 inner and outer regions of the boundary layer. The underprediction of mean velocity by the URDNS is seen in this case and was fairly consistent among all cases. The similarity model also greatly underpredicted the DNS results and is actually slightly less than the URDNS results over part of the boundary layer. The dynamic model performed rather poorly in this case, giving an overprediction of the mean velocity in the logarithmic region. For this case the dynamic model's predictions also coincide with results published by Park and Mahesh [61]. Figure 6.4: Mean velocity proles: Re =180, Mesh: (32, 32, 33). For comparison, a repeat of theRe =180 case using a higher resolution is also shown. This case was conducted with resolution nearly doubled in all coordinate directions (64, 100 64, 65). As expected, the URDNS results improved somewhat with this increase in resolution, although these results are still well below the DNS curve. Despite the higher resolution, the similarity model underpredicted the mean velocity and generally did not oer an improvement over the results obtained by simply running the code with no model enabled. The predictions from the similarity model show a marginal improvement with the resolution increase that was in line with what was seen in the URDNS results. At this higher resolution the dynamic model shows noticeable improvement. The dynamic model still exhibits a considerable overprediction of the mean velocity; however, it is not as marked as in the previous case. Again, the interscale transfer model gave a good prediction of the mean velocity using the identical lter setting established in the low- resolution case. Aside from an improvement in the boundary layer's outer region, this model shows relatively little change in comparison to its coarse resolution results. At both resolutions the proposed model reasonably reproduced the DNS prole for the Re =180 case. This consistency in the predictions of the proposed model (despite a doubling of the mesh resolution) provides the rst evidence that the model's performance does not have an excessively high sensitivity to mesh resolution. In fact, at this Reynolds number and beyond, the present model demonstrates a noted lack of excessive sensitivity to changes in resolution, aside from what is expected from typical SGS model behavior. This point is made due to concerns that the model's inherent sensitivity to the test lter strength and the test lter's implicit dependence on resolution would lead to excessive deterioration of the model's results under changes to mesh resolution. This concern has not proven to be an issue. 101 Figure 6.5: Mean velocity proles: Re =180, Mesh: (64, 64, 65). Simulations were carried out using a smaller domain size at Re =1050, which corre- sponds to a physical Reynolds number Re c 24; 000. Again, the URDNS results gave a large underprediction of the mean velocity prole. In this case the similarity model continued to underpredict the mean velocity prole given by the DNS data. Despite this considerable error, the similarity model appeared to perform its best at this Reynolds number. The dynamic model gave an overprediction of the mean velocity. In this case the present model performed reasonably well, overlapping with the DNS results, and performing considerably better than both the dynamic and similarity models. 102 Figure 6.6: Mean velocity proles: Re =1050, Mesh: (96, 128, 65). The Re =1050 case was repeated with resolution (128, 128, 97), which more closely falls into the recommended resolution ranges. The resolution in both the streamwise and wall-normal directions was increased. The URDNS results changed only marginally with this change in resolution. The similarity model returned reasonable results and gave an improvement over the previous resolution and oered a signicant improvement over the URDNS values. However, this model continued to underpredict the DNS data. The dynamic model returned results that were nearly identical to those found with the coarse resolution. The dynamic model's level of overprediction was on par with the similarity 103 model's level of underprediction. At this resolution the proposed model returned values that were a very good match with the DNS. Figure 6.7: Mean velocity proles: Re =1050, Mesh: (128, 128, 97). TheRe =2000 case (corresponding toRe c 50,000) was the highest Reynolds number considered. The LES cases used the same relatively small computational domain as in the Re =1050 cases, although the resolution expressed in wall units was somewhat dif- ferent. The similarity model again fell short of the DNS proles returning a considerable underprediction of the mean velocity. The dynamic model again gave an overprediction of the DNS results although the degree of this overprediction was slight. Overall, the 104 model's prediction was reasonably good throughout the entire boundary layer. The pre- dictions obtained from the present model were also quite good and the results overlapped considerably with the DNS deviating only slightly in the buer region. Figure 6.8: Mean velocity proles: Re =2000, Mesh: (256, 96, 131). The small domains used in the high Reynolds number cases are of course capable of limiting the ow in its low wavenumber content and can eliminate the action of those modes approaching or exceeding the size of the domain. This limitation is in contrast to the direct numerical simulations. The Re =950 case (Re c 21,500) is of interest because it considers a domain size with a notable increase over the other cases and identically matches the large domain used in the DNS case [17]. 105 Figure 6.9: Mean velocity proles: Re =950, Mesh: (512, 384, 85). Given the similarity model's poor predictions in previous cases, the model was not implemented in theRe =950 case. In this case the dynamic model gave an improvement over the URDNS results, but again somewhat overpredicted the mean velocity. The present model performed well in this simulation providing a prediction of mean velocity that compared well with the DNS results. In a simulation of this size, DNS is often a prohibitively expensive alternative. However, the LES resolution in this case resolved only a fraction of the modes captured in the DNS ( 0:6%). At a fraction of the computational expense the simulation using the proposed interscale transfer model has provided mean velocity predictions that very nearly reproduce those obtained from the DNS. 106 As expected, in all Reynolds number considered, practically all of the LES models accurately captured the linear velocity prole in the viscous sublayer. However, in the buer layer and throughout the logarithmic layer the models diered considerably in their results. The similarity model generally did a poor job of predicting the mean velocity throughout the outer portion of the boundary layer for all cases in which it was implemented. In some cases the similarity model did not provide any improvement over the predictions obtained by running the code without a SGS model enabled. The consistent underprediction of the mean velocities by the similarity model in the present wall-bounded ow simulations support the previous claims of the model's weaknesses and underscores its general inadequacy for use in LES applications. In contrast to the similarity model, the dynamic model tended to slightly overpredict the results of the direct numerical simulations. Predictions for the low resolution case withRe =180 proved to be problematic for the dynamic model, which for this case, gave a severe upward shift of the mean velocity prole in the logarithmic region. The dynamic model tended to overpredict the DNS proles and in at least one case the magnitude of the model's error was on the order of that obtained by the similarity model. However, with the exception of the low Reynolds number case, the dynamic model generally gave reasonable predictions of the mean velocities. The proposed model performed relatively well for all cases considered. Its performance was found to be consistently superior to the similarity model and at least similar in quality and often better than that of the dynamic model. 107 6.5.2 RMS Velocity Results Due to their random nature the uctuating components of velocity in a turbulent bound- ary layer are best represented as averages in space and/or time (see Eq. 6.10). For wall- bounded ows the character of the RMS values of the uctuations have been researched extensively both experimentally and numerically. Plots of mean turbulence intensities across an inner section of the boundary layer are shown for each case discussed above. The DNS data is also plotted. In these ows turbulent kinetic energy is generally concen- trated inside the boundary layer where production due to mean shear is most prevalent. The presence of the wall and the alignment of the freestream gives distinct character to the components of these velocities and their combined contribution to the the regions of elevated turbulent kinetic energy inside the channel. The streamwise component of uctuating velocity u 0 maintains an elevated level of intensity relative to the other com- ponents because it readily receives energy production from the freestream kinetic energy. Thev 0 andw 0 components cannot extract energy directly from the mean ow, but rather they obtain energy from u 0 through pressure-velocity interaction terms present in the energy budget equation (Tennekes and Lumley [80]). The locations of peak production of energy in these components occurs near the location of maximum u 0 . This energy, however, is not shared equally among the v 0 and w 0 components. The w 0 component is aligned with the wall-normal direction and the production of this component is impeded somewhat by the in uence of the wall (White [87]). Lacking this direct in uence of the wall, the v 0 component is somewhat higher than w 0 . Analyzing the kinetic energy and 108 Reynolds stress budgets, along with turbulence shear stress u 0 w 0 , etc. can give further insight into the behavior of various models. The RMS velocities were computed using averaging over horizontal planes and were then normalized by the mean friction velocity. As with mean streamwise velocity, the RMS velocities were averaged by sampling over several time steps. The following plots show these terms plotted across a section of the channel closest to the wall. The general shape of the LES and DNS curves are expected to be similar. Their magnitudes should also be close; however, it should be noted that LES predictions of RMS values of velocity are expected to slightly underpredict DNS calculations of these same quantities. The reason is that direct numerical simulations fully resolve the ow-eld and so all relevant contributions to the uctuating components of velocity are accounted for. In LES, the ow-eld is under-resolved and so only a portion of the real ow uctuations are present. However, since LES resolves the large energy-carrying eddies which account for the large contributions to the RMS velocities, the LES and DNS curves are expected to be relatively close in magnitude. With the exception of the high resolution Re =180 case, the DNS results have not been ltered to the LES resolution; therefore, the DNS results provide a qualitative comparison with results from the simulations and should realistically give an upper bound for the intensities obtained in the LES. In terms of the magnitude of the RMS plots, a better agreement of the magnitude of the LES and the DNS proles would be possible if an appropriate lter could be applied to the DNS data and these results were then sampled at the LES resolution. However, the approach often taken in literature is to compare the LES results to the unltered DNS data (Vreman [83]). Nonetheless, the general shape of the LES proles are expected 109 to match that given by the DNS, and the magnitudes of these curves can reasonably be expected to fall within approximately 15% of the DNS proles. In any case, the LES predictions of RMS velocities should be slightly lower than the DNS data. The eect of mesh resolution on RMS velocities becomes more clear by viewing results from the rst two cases presented. These cases compare RMS velocities at identical Reynolds numbers for a mesh resolution doubled in all coordinate directions. Along with unltered DNS results, the high resolution Re =180 case also shows DNS results that have been ltered and sampled at the LES resolution. The RMS velocities for theRe =180 case with resolution (32; 32; 33) is shown in Fig. 6.10. In this case the similarity model considerably underpredicts the DNS proles for all components of RMS velocity. The dynamic model gives a substantial overprediction of the streamwise component which extends into the core of the channel. The present model gives only a slight underprediction of the streamwise component. An underprediction of the DNS data at this level is expected from LES. All models capture the general shape of the spanwise and wall-normal curves. A notable feature in the spanwise and wall-normal plots is the relatively low predictions given by the dynamic model in regions near the peaks of these curves. Given that the dynamic model overpredicted the streamwise RMS velocity and underpredicted the spanwise and wall-normal components, suggests that the model's energy exchange process between component directions may be decient at this resolution. To provide a more realistic comparison with the RMS predictions by the LES, the DNS results computed in the present eort were ltered (using a trapezoidal approximation of the box lter) and sampled at a resolution matching the LES. The ltered DNS predictions 110 Figure 6.10: RMS velocity proles: Re =180, Mesh: (32, 32, 33). (referred to here as FDNS) of the RMS velocities are plotted along with their unltered counterpart for the Re =180 case with resolution (64; 64; 65) in Fig. 6.11. The ltering and sampling operations account for the noticeable reduction in the FDNS proles. Clearly, in both the Re =180 cases the dynamic model exhibits its known weakness in that it tends to overpredict the streamwise velocity uctuations [61]. In the low resolution case this overprediction is accompanied by somewhat lower activity in the spanwise and wall-normal directions. The dynamic model's level of discrepancy was 111 Figure 6.11: RMS velocity proles with ltered DNS: Re = 180, Mesh: (64; 64; 65). All LES, DNS, and FDNS results provided by the author. reduced considerably as the resolution was increased, although its streamwise prediction is perhaps still unusually high given its proximity to the unltered DNS curve. Doubling the resolution generally improved the dynamic model's predictions of spanwise and wall- normal RMS velocity components. With the exception of the dynamic model's tendency to overpredict the RMS velocity in the streamwise direction, both the present model and the dynamic model performed relatively well for resolution (64, 64, 65). The present model provided perhaps the best match with the FDNS results. 112 In both the high and low resolution cases at this Reynolds number, the similarity model gives comparatively low predictions of nearly all component activity, most notably in the streamwise velocity. A doubling of the mesh resolution appeared to have little eect on this model's performance. This consistent underprediction of RMS velocities was typical of the similarity model for practically all cases considered. Figure 6.12: RMS velocity proles: Re =1050, Mesh: (96, 128, 65). 113 Since DNS results at Re =1050 were not available, DNS statistics at Re =950 are used for comparison with the LES at Re =1050. This discrepancy in Reynolds numbers is not an issue when comparing mean velocities since the boundary layer proles are self- similar. However, in comparing RMS velocities the LES results at Re =1050 should be expected to have a slightly higher magnitude than what would be present at Re =950. Thus, the discrepancy in DNS Reynolds numbers manifests itself as a reduced magnitude of RMS velocities in the DNS results. Figure 6.13: RMS velocity proles: Re =1050, Mesh: (128, 128, 97). 114 For the Re =1050 cases, the results at higher Reynolds numbers are not plotted across the entire half-width of the channel since these proles tend to be fairly at as they approach the center of the channel. In the low resolution case, the streamwise RMS results for the present model generally fell between those of the similarity model and the dynamic model, and were often higher than both in the predictions of the wall-normal component. In the streamwise component, the present model closely matches the DNS results over much of the boundary layer. Given that in this case the DNS data corresponds to the lower Reynolds number (Re =950) simulation, the LES curve's close proximity to the unltered DNS curve is somewhat expected. Comparing the dynamic model's predictions of streamwise RMS velocity, again suggests that there is an overprediction; however, this point cannot be made with certainty for this case given the LES/DNS Reynolds number discrepancy. The performance of the proposed model is again consistent with previous cases. Overall the interscale transfer model has generally predicted the magnitude and shape of the RMS curves reasonably well. On the other hand, at this Reynolds number the similarity model underpredicts the peaks of all RMS components. The proposed model and the dynamic model performed similarly after increasing the streamwise and wall-normal resolution. In general, at this resolution the similarity model continued to underpredict the magnitudes of the DNS data. In the Re =2000 case, the dynamic model's prediction of streamwise RMS values overlapped the DNS results near the wall and out to z + 70. Beyond this the dynamic model tended to underpredict the DNS prole. Conversely, the similarity model un- derpredicted the streamwise DNS in the near-wall region and gave results that closely 115 Figure 6.14: RMS velocity proles: Re =2000, Mesh: (256, 96, 131). overlapped the DNS prole away from the wall. Again, the similarity model consider- ably underpredicted the spanwise and wall-normal activity throughout the channel. The present model's predictions of the wall-normal and spanwise components closely matched the results obtained from the dynamic model, and in general the present model gave a consistent level of underprediction of the DNS for all RMS components throughout the boundary layer. This model's underpredictions were generally less exaggerated compared 116 with the similarity model, and its proles tended to be within a range that is expected from LES when compared to unltered DNS. Figure 6.15: RMS velocity proles: Re =950, Mesh: (512, 384, 85). For the Re =950 case, the dynamic model's prediction of streamwise RMS velocity again appear slightly high near the peak; whereas, the present model resides slightly below the DNS results. With the exception of the slight overprediction of streamwise activity by the dynamic model, the present model and the dynamic model performed similarly in all coordinate directions. 117 As anticipated, the results from the proposed model tended to slightly underpredict the unltered DNS values of RMS velocities. At times the dynamic model tended to predict unrealistically high RMS values that matched or even exceeded the DNS results. Aside from this discrepancy, the results of these two models were often similar in their proles. The proposed model and the dynamic model generally outperformed the scale- similarity model in these simulations. The important conclusion is that the presently proposed formulation of the alternative similarity-type model oers obvious improvements over the predictions of mean and RMS velocities obtained from the classic similarity model. 6.5.3 Energy Dissipation Results Similar to the LES equations which are obtained by ltering the Navier-Stokes equations, a ltered version of the energy equation can be obtained by applying a lter to the inner product of velocity with the Navier-Stokes equations, i.e. u i @u i @t +u i @u i u j @x j = 1 u i @p @x i +u i @ 2 u i @x 2 j : (6.13) Assuming ltering and dierentiation commute one obtains 1 2 @ (u i u i ) @t + 1 2 @ (u i u i u j ) @x j = 1 @ (pu i ) @x i +u i @ 2 u i @x 2 j : (6.14) A straightforward method for obtaining an equation governing the transport of re- solved kinetic energy is to take the inner product of resolved velocity with the LES equations. The resulting equation governs resolved energy transfer and can be written as 118 u i @ (u i ) @t +u i @ (u i u j ) @x j = 1 u i @ (p) @x i +u i @ 2 u i @x 2 j u i @ ij @x j : (6.15) Reformulating terms, the following equation is obtained 1 2 @ (u i u i ) @t + 1 2 @ (u i u i u j ) @x j = 1 @ (pu i ) @x i + 2 @ 2 (u i u i ) @x 2 j @u i @x j @u i @x j @ (u i ij ) @x j + ij S ij ; (6.16) where S ij is the strain-rate tensor constructed from resolved velocities. Subtracting the resolved energy equation (6.16) from the ltered version of the full energy equation (6.14) gives the transport equation for subgrid-scale kinetic energy, i.e. 1 2 @ ( ii ) @t + 1 2 @ ( ii )u j @x j = 1 2 @ (u i u i u j u i u i u j ) @x j 1 @ (pu i pu i ) @x i + 2 @ 2 ( ii ) @x 2 j @u i @x j @u i @x j @u i @x j @u i @x j + @ (u i ij ) @x j ij S ij : (6.17) In this equation the last two terms contain the eects of the SGS stress tensor on the transport of resolved and SGS energy. These terms appear (with opposite signs) in both the resolved and SGS energy equations and they represent the diusion and dissipation of energy due to the SGS stresses, respectively. The SGS diusion term is expressed in ux-divergence form limiting its energy contributions to the domain boundaries. Globally 119 this term does not account for any net energy transfer in the equations for the resolved and subgrid scales. On the other hand, local non-zero values of SGS dissipation can result in a net energy exchange between the resolved and subgrid scales. Energy dissipation is naturally a topic of consideration in LES model development. A SGS model can aect activity across all resolved scales so a model's full impact on dissipation occurs directly through SGS dissipation and indirectly through its eects on scales undergoing the eects of viscosity. A SGS model inevitably impacts both the subgrid-scale dissipation and resolved scale dissipation in a posteriori analyses. The value of a simulation's resolved scale dissipation is represented as " visc and is dened as " visc 2S ij S ij : (6.18) This term, along with the subgrid-scale dissipation dened as " sgs ij S ij ; (6.19) are considered in the present eort. The values of viscous and SGS dissipation were ob- tained by integrating their planar averaged proles across the channel and were computed for all LES cases discussed previously. For reference, predictions of dissipation obtained by employing the Smagorinsky model have also been included for the high resolution Re =180 case and for the low resolution Re =1050 case. Viscous dissipation is the irreversible conversion of a ow's kinetic energy to inter- nal energy caused by uid deformation (Panton [59], Kundu [41]). This shows up as a 120 negative value in the energy equation representing a continual energy sink. The SGS dissipation represents the energy exchange between the resolved and subgrid scales in the ow. Locally this term can assume either negative or positive values. The general direction of energy transfer in turbulence is from small to large wavenumbers so this most prominent exchange process is referred to as forward transfer. However, reverse transfer (backscatter) also occurs when the energy of large wavenumber modes causes increases in the energy of modes with small wavenumbers. Negative values of SGS dissipation correspond to forward transfer representing a loss of energy from the resolved scales to the subgrid scales, while positive values represent reverse transfer and a gain of energy in the resolved scales. The dynamic model, which is often implemented using a computed coecient for- mulated using planar-averaged terms in its numerator and denominator [32], is generally restricted to giving only forward transfer predictions of SGS dissipation. The Smagorin- sky model, which has a tensor formulation giving an eigenvector alignment with that of the resolved strain rate tensor is also limited to predictions of forward scatter. A per- ceived strength of the similarity model, however, is its ability to permit backscatter. This property is also exhibited by the present model. To give examples of the present model's ability to predict both forward and reverse transfer, snapshots showing the SGS dissipation eld predicted by the interscale transfer model in the high resolutionRe =180 case and in theRe =2000 case are shown in Figures 6.16 and 6.17, respectively. Regions of forward and reverse SGS energy transfer show up as blue and red regions, respectively. Similar contour plots were made for the other LES 121 cases, the main dierence being the change in the prevalence of scale and magnitude of the dissipation regions with Reynolds number. Figure 6.16: Instantaneous subgrid-scale dissipation eld: Re =180 (HR). Walls are at top and bottom of image. Flow is from left to right. In the resolved energy equation, the SGS diusion term represents a redistribution of energy due to the SGS stresses. In LES, this term locally redistributes the resolved scale energy but does not alter the total energy content in the ow. On the other hand, the SGS dissipation term accounts for the net energy transfer between the resolved and subgrid scales. Given that the small scales are not resolved in LES, this term is capable of altering the total resolved energy contained in the ow. Arguably the most important role of a SGS model is the dissipation of resolved scale energy accounting for the absence of the viscous dissipation of the smallest scales in LES. For this reason, it is a model's energy dissipation eects that are of major interest. For illustration, sample planar-averaged proles of SGS 122 Figure 6.17: Instantaneous subgrid-scale dissipation eld: Re =2000. Walls are at top and bottom of image. Flow is from left to right. and resolved scale dissipation corresponding to the high resolutionRe =180 case and the Re =2000 case are shown in Figures 6.18, 6.19, 6.20, 6.21. These proles show that generally the planar-averaged proles of net SGS dissipation were comparable for the dynamic model and the interscale transfer model. This similarity is worthy of comment given that these two models take a considerably dierent approach in arriving at this result. The forward transfer predicted by the similarity model was, in comparison, quite high and resulted in a net SGS energy transfer that was considerably greater than the other models. In 3-D simulations with no SGS model there is a gradual build-up of activity in the vicinity of the mesh cuto. The accumulation of energy at these scales occurs until the action of viscosity across all scales adequately balances the dissipation requirements of the 123 Figure 6.18: Forward, reverse, and net subgrid-scale dissipation proles: Re =180 (HR). energy being generated in the turbulence. The integrated viscous dissipation in URDNS can provide a practical estimate for the upper bound of total dissipation expected from LES. In Figures 6.20 and 6.21, the similarity model's peak values in viscous dissipation rival those predicted by the under-resolved DNS. The similarity model's viscous dissipa- tion prole also tends to increase somewhat more rapidly than that of the dynamic and interscale transfer models suggesting a dierent resolved scale dissipation trend in the vicinity of the wall. In the LES results, single quantities summarizing reverse and forward SGS energy transfer were computed for all cases by individually integrating the positive and negative proles of planar-averaged SGS dissipation. A resulting net forward transfer is expected in these calculations. The net SGS dissipation and total dissipation (SGS plus viscous) 124 Figure 6.19: Forward, reverse, and net subgrid-scale dissipation proles: Re =2000. values were also computed. The ratios of net SGS dissipation to total dissipation have also been calculated. These values are shown in Table 6.3 for all LES cases considered. In these comparisons the forward transfer predicted by both the interscale transfer model and the similarity model tended to be quite high. However, the present model had an accompanying high level of reverse transfer so it generally had somewhat lower values of net SGS dissipation. The SGS dissipation predicted by the present model was also generally less than that of the dynamic model, although the total dissipation (resolved plus net subgrid) for these models tended to be close given their respective levels of resolved scale dissipation. 125 Figure 6.20: Mean viscous dissipation proles: Re =180 (HR). The Smagorinsky model has often been criticized for having too much dissipation (Vreman et al. [86]), although its results depend strongly on the model's selected coe- cient (c s ). Here,c s values of 0.1 and 0.06 have been used in the Re =180 andRe =1050 cases, respectively. The SGS dissipation predicted by the Smagorinsky model tended to be in the same range as the other models tested here. On the other hand, the similarity model has earned a reputation of providing an inadequate amount of SGS dissipation (Bardina et al. [2], Liu et al. [48]). However, Vreman et al. [86] and Vreman [82] observed that the similarity model is inadequate, specically in its predictions of SGS dissipation 126 Figure 6.21: Mean viscous dissipation proles: Re =2000. in the smallest resolved scales. Based on the present testing, the latter of these state- ments seems to be more accurate. In all cases, the similarity model's SGS dissipation exceeded the values of the dynamic model and the present model. Also, in nearly all cases the similarity model provided a total dissipation that rivaled that of the URDNS results. The exception to this is in the Re =1050 cases where the similarity model's total dissipation was considerably less than the URDNS. Perhaps not coincidentally, this is Reynolds number in which the similarity model provided its best predictions of mean velocity. A model's eect on energy dissipation occurs directly through its presence in the SGS dissipation term, but also indirectly impacts viscous dissipation through its role in altering 127 the activity of the resolved scales. Generally, it is not possible to decouple a model's role in resolved scale and subgrid-scale dissipation. Although intuitively, as mesh resolution is increased (holding Reynolds number xed) a model's SGS dissipation is expected to have a diminishing role in the total energy dissipation. By comparing values of net SGS dissipation (" net sgs ) to total dissipation (" visc +" net sgs ), it is evident that for the Re =180 case that this trend holds for both the dynamic and interscale transfer models. However, the reverse of this was observed for the scale-similarity model. On the other hand, this expected behavior was observed for all LES cases in the Re =1050 simulations. Liu et al. [48] claims that a strength of the similarity model is its applicability across dierent ow regimes. This statement is based on the idea that the similarity model, formulated from the tensor L ij , tends to give smaller values of its tensor components when the cell Reynolds number is decreased in a priori testing. Thus, the role of the similarity model should tend to diminish as a ow's Reynolds number is decreased or a simulation's resolution is increased. This feature is desirable and one that some eddy- viscosity models do not benet from|the dynamic model being an exception since it also contains a dependency on L ij . Despite not having a direct dependence on L ij , the interscale transfer model also exhibits this property. The present model simply contains partial components of stress that vanish in a fashion similar to those contained in L ij . Models not satisfying this condition, for instance, can tend to overpredict SGS energy dissipation at low Reynolds numbers. As discussed above, a posteriori testing of the similarity model has shown that this condition was not consistently exhibited for all cases considered. This is another example where the similarity model has shown a tendency 128 to deviate from the expected behavior of the subgrid scales. On the other hand, the interscale transfer model and the dynamic model have consistently shown this behavior. It is possible to decompose the similarity model's integrated SGS dissipation into components, each owing to the model's incorporated nonlinear terms. Parsing the model's aggregate SGS dissipation is useful given that the dissipation originating from the model's terms aectingR 1 orR 2 can be identied. This breakdown of the similarity model's SGS dissipation is shown in Table 6.4 for the various simulations. Decomposing the SGS dissipation in this way shows a split where approximately 65- 80% of the model's net SGS dissipation is credited to terms aecting energy change in the largest resolved scales (i.e R 1 ). Relative to other models, the similarity model provides dissipation levels that are substantial; however, it is the model's excessive contribution to SGS dissipation from terms aecting the largest resolved scales that likely drives the relationship between SGS dissipation and resolved scale dissipation out of its correct balance. In contrast to this, the present model is constructed entirely from terms aecting energy transfer in R 2 . Thus, all of this model's SGS dissipation can be associated with energy removal from terms aecting the smallest resolved scales|the implementation of a smooth lter being the only factor that relaxes this principle somewhat. Thus, the proposed arrangement of terms seems to more realistically match the natural energy exchange process observed in real turbulence; specically, the greatest contribution to the model's SGS dissipation is associated with energy removal from small scales in the vicinity of the LES cuto. 129 Model " visc " fwd sgs " rev sgs " net sgs " visc +" net sgs " net sgs =(" visc +" net sgs ) Re =180 (LR): None -37.16 0 0 0 -37.16 0 Similarity -30.25 -10.43 3.78 -6.66 -36.91 0.18 Dynamic -24.52 -4.74 0 -4.74 -29.26 0.162 Interscale Transfer -26.67 -9.99 7.03 -2.97 -29.64 0.10 Re =180 (HR): None -38.89 0 0 0 -38.89 0 Similarity -31.02 -9.24 1.61 -7.63 -38.65 0.197 Dynamic -29.26 -3.82 0 -3.82 -33.08 0.115 Interscale Transfer -31.82 -7.17 3.88 -3.29 -35.11 0.0937 Smagorinsky -29.05 -4.19 0 -4.19 -33.24 0.126 Re =950: None -37.75 0 0 0 -37.75 0 Dynamic -25.25 -7.17 0 -7.17 -32.43 0.221 Interscale Transfer -28.62 -16.53 11.76 -4.77 -33.38 0.143 Re =1050 (LR): None -40.25 0 0 0 -40.25 0 Similarity -26.29 -14.44 4.62 -9.82 -36.11 0.272 Dynamic -24.09 -7.62 0 -7.62 -31.71 0.24 Interscale Transfer -28.59 -17.77 13.25 -4.51 -33.11 0.136 Smagorinsky -31.86 -5.1 0 -5.1 -36.96 0.138 Re =1050 (HR): None -38.0 0 0 0 -38.0 0 Similarity -25.57 -12.23 3.01 -9.22 -34.79 0.265 Dynamic -25.20 -6.15 0 -6.15 -31.35 0.196 Interscale Transfer -28.82 -14.26 10.0 -4.26 -33.08 0.129 Re =2000: None -31.81 0 0 0 -31.81 0 Similarity -22.66 -12.67 3.44 -9.22 -31.88 0.289 Dynamic -19.97 -6.87 0 -6.87 -26.84 0.256 Interscale Transfer -23.43 -17.55 13.18 -4.37 -27.79 0.157 Table 6.3: Integrated values of resolved scale dissipation, forward, reverse, and net SGS dissipation, total dissipation, ratio of net SGS dissipation to total dissipation. Case " net sgs " R1 sgs " R2 sgs " R2 sgs =" net sgs Re =180 (LR) -6.66 -5.24 -1.42 0.21 Re =180 (HR) -7.63 -5.17 -2.46 0.32 Re =1050 (LR) -9.82 -6.66 -3.16 0.32 Re =1050 (HR) -9.22 -6.01 -3.21 0.35 Re =2000 -9.22 -6.61 -2.61 0.28 Table 6.4: Decomposition of the similarity model's integrated SGS dissipation. Total net SGS dissipation, contribution to net SGS dissipation associated with terms aecting interscale energy transfer in R 1 , and in R 2 , ratio of net sgs dissipation associated with terms aecting R 2 to total net SGS dissipation. 130 Chapter 7 Additional Considerations 7.1 Non-Symmetry of the Proposed SGS Stress Tensor As previously mentioned, the proposed model is not symmetric in its present form. This is of course in contrast to the exact form of the SGS stress tensor. Non-symmetric models are not foreign to the eld of LES and have been proposed in previous research including the work of Domaradzki and Holm [21] who formulated the alpha model, and Geurts and Holm [33] who implemented a regularization modeling approach implementing the non-symmetric Leray model. In the LES equations, the exact form of the SGS stress tensor has the role of imparting local stresses representing the eects of the unresolved scales of turbulence. Similarly, the viscous stress tensor formulated from the resolved strain rate, accounts for viscous stresses of the resolved scales. As a fundamental property, when a uid is at rest the shearing components of stress must vanish (Kundu [41]). This is satised by the exact form of the SGS stress tensor given its dependence on velocity, and the viscous stress tensor given its relation to the strain rate. The proposed SGS stress tensor also satises this fundamental 131 requirement (including vanishing normal stresses) given that it is constructed solely from the resolved velocity components. As uid moves it develops stresses due to the presence of viscosity and in most cases turbulence is generated. However, in simple non-turbulent ows such as solid body ro- tation or simple shear ow (having a constant velocity gradient), viscous stresses result in pure normal stress or pure shear stress, respectively. In laminar ows such as these the relevant physical scales are usually on the order of those dened by the geometry and typical LES models will predict SGS stress components that are very nearly zero under these conditions. The exact form of the SGS stress tensor again does not introduce any articial forces under these conditions. This vanishing contribution is not a trivial feature of all SGS models, although this property has been veried for the current model in laminar wall-bounded ow. The Smagorinsky model in particular has been criticized for violating this condition (Park et al. [60]). A uid volume of nite size can experience a change in its rotation rate due to the presence of a gradient in the local shear stress. The gradient results in a force imbalance across the particle that imparts an angular acceleration. As the size of the uid volume is allowed to shrink the eect of the gradient on the uid particle approaches zero and a force balance is attained that must result in a zero net angular acceleration of the particle. The viscous stress tensor describes the viscous stresses at a point of innitesimal size in the ow and is thus a symmetric tensor. Note that this is the case for isotropic and Newtonian uids, which accounts for most uids in general applications (White [87]). Given that the exact form of the SGS stress tensor is similarly a symmetric expression, this symmetry guarantees that the SGS stresses acting on a uid particle of innitesimal 132 size will not experience angular acceleration at the hand of SGS forces. The modeled SGS stress tensor presently proposed is not, however, a symmetric expression. This non- symmetry can be displayed by expressing the full form of the model's stress tensor and comparing its o-diagonal components, i.e. II ij = 2 6 6 6 6 6 6 6 4 [ b u 1 u 1 b b u 1 b b u 1 b u 1 u 0 1 [ b u 1 u 2 b b u 1 b b u 2 b u 1 u 0 2 [ b u 1 u 3 b b u 1 b b u 3 b u 1 u 0 3 [ b u 2 u 1 b b u 2 b b u 1 b u 2 u 0 1 [ b u 2 u 2 b b u 2 b b u 2 b u 2 u 0 2 [ b u 2 u 3 b b u 2 b b u 3 b u 2 u 0 3 [ b u 3 u 1 b b u 3 b b u 1 b u 3 u 0 1 [ b u 3 u 2 b b u 3 b b u 2 b u 3 u 0 2 [ b u 3 u 3 b b u 3 b b u 3 b u 3 u 0 3 3 7 7 7 7 7 7 7 5 : (7.1) The non-symmetry of this tensor suggests that the use of the model can lead to the generation of angular acceleration of a uid volume of nite size in the absence of a gradient in shear stress. This is generally not the case since the model's non-symmetry reduces to a dierence in ltered and unltered quantities, which would tend to vanish in the absence of a gradient. However, the SGS expression does permit contributions to angular momentum due to local ow properties. It has been shown that these eects cannot contribute to the net angular momentum of the ow and steps have been taken to demonstrate that these local redistribution eects are small and are of negligible im- portance relative to the numerical accuracy of typical CFD codes. The implications of the non-symmetry have been quantied both analytically and numerically. The results of this are reported in the present section. In LES, the SGS stress along with the resolved strain rate combine to contribute to SGS dissipation. In the resolved energy equation, the non-symmetric part of the proposed 133 model's SGS stress tensor does not contribute to SGS dissipation given that this quantity is computed as a contraction with the resolved strain rate, which is a symmetric tensor by denition. The non-symmetric part of SGS stress tensor can, however, aect SGS diusion. Through model testing and scaling considerations, it has been shown that the role of the model's non-symmetry in SGS diusion is quite small in magnitude and tends to be distributed randomly inside the computational domain. Analytical and numerical considerations of the non-symmetry quantify the imbalance in the model's o-diagonal components and help describe any resulting local introduction of angular momentum due to this condition. A series expansion of the terms in the model has been used to determine the comparative dierences in the o-diagonal terms. The order of the dierence of these terms is shown to be quite small relative to those in the model itself. To analytically quantify the imbalance in the SGS shear stresses, the dierences in the o-diagonal components are computed for each term in the model's expression. This dierence is expressed as II ij II ji = h d b u i u j b b u i b b u j + b u i u 0 j i h d u i b u j b b u i b b u j +u 0 i b u j i : (7.2) Noting that the N 112 term is symmetric and does not contribute, this dierence can be expressed more conveniently as II ij II ji = h u i b u j d u i b u j i h b u i u j d b u i u j i : (7.3) 134 From the expression above it is clear that the model's non-symmetry has contributions only from the dierences of ltered and unltered quantities. Given that the test ltering operation is a convolution of the form given by Eq. (3.23), specically b f (~ x) = Z +1 1 G ~ x ~ x 0 ; b f ~ x 0 d ~ x 0 ; (7.4) integrating the ltered expression over the domain and changing the order of integration taken on the right side of this expression gives Z +1 1 b f (~ x)d~ x = Z +1 1 Z +1 1 G ~ x ~ x 0 ; b f ~ x 0 d~ xd ~ x 0 : (7.5) If f and G are integrable functions, a property of convolution integration permits the integral of the convolution to be expressed as the product of two integrals Z +1 1 b f (~ x)d~ x = Z +1 1 G ~ x ~ x 0 ; b d~ x Z +1 1 f ~ x 0 d ~ x 0 : (7.6) From Eq. (5.13), the integral of G on the right side is equal to unity. Therefore, the following relation is obtained for the integration of a ltered quantity Z +1 1 b f (~ x)d~ x = Z +1 1 f (~ x)d~ x: (7.7) The above equality shows that the global integration of ltered and unltered quan- tities gives an identical result. This has implications with regard to the present model's impact on angular momentum. Specically, the imbalance in the model's o-diagonal 135 components of stress, consisting of dierences in ltered and unltered quantities, are re- stricted to having only local eects on angular momentum. Globally angular momentum is conserved and there is no net eect on this due to an imbalance in the model's SGS shear stress. Figure 7.1: Mean SGS diusion proles: Re =180 (HR). All LES and DNS results pro- vided by the author. The present model is also limited to having only local eects on energy transfer due to SGS diusion since this quantity appears in ux-divergence form in the resolved energy equation. A sample case showing the planar-averaged SGS diusion proles predicted by the various SGS models considered in this research is shown in Fig. 7.1. Also shown in this gure is a SGS diusion prole computed from the DNS. This prole was generated 136 by constructing the expression for SGS stress from the DNS results using a trapezoidal approximation of the box lter in place of the implicit LES lter. The mean proles in this gure all integrate to zero. From these results it is evident that the present model's formulation does not produce eects on SGS diusion that are out of family from the other SGS models considered. The eect of the model's non-symmetry in LES is limited to a redistribution of energy and angular momentum. The local impact of this can be quantied by carrying out series expansions of terms in the model. Applying the discrete test ltering operation to products of ltered velocity components introduces ltered velocities located at a central point n and at neighboring mesh points, i.e. b u n1 i =a 1 u n2 i +a 0 u n1 i +a 1 u n i ; (7.8) b u n i =a 1 u n1 i +a 0 u n i +a 1 u n+1 i ; (7.9) b u n+1 i =a 1 u n i +a 0 u n+1 i +a 1 u n+2 i : (7.10) For simplicity, ltering in only one spatial direction is shown. Thus, the term d b u i u j is expressed as d b u i u j =a 1 b u n1 i u n1 j +a o b u n i u n j +a 1 b u n+1 i u n+1 j : (7.11) 137 The other terms included in Eq. (7.3) are expressed in a similar way. The Taylor series expansions of the velocities present in these expressions are written as u n2 i = 1 X k=0 (2x) k k! @ k u n i @x k ; (7.12) u n1 i = 1 X k=0 (x) k k! @ k u n i @x k ; (7.13) u n+1 i = 1 X k=0 (x) k k! @ k u n i @x k ; (7.14) u n+2 i = 1 X k=0 (2x) k k! @ k u n i @x k : (7.15) Expanding all terms in the dierence expression, one obtains II ij II ji = 2a 2 1 x 4 @ 3 u n i @x 3 @u n j @x @u n i @x @ 3 u n j @x 3 ! +a 2 1 x 4 @ 4 u n i @x 4 u n j u n i @ 4 u n j @x 4 ! +O x 6 : (7.16) The leading-order terms that remain after computing this dierence are fourth-order in space, where the coecienta 1 corresponds to the lter coecient 1 12 selected for the present model. Thus, the leading-order coecient in this dierence expression is quite small 1 72 . On the other hand, a similar procedure shows the leading terms in the model itself are second-order with a leading-order coecient of 1 6 . For this reason, one would 138 expect that simulations conducted with numerical solvers accurate up to O x 4 should not realize eects introduced by the model's non-symmetry. Based on this comparison and considering the relatively small leading-order coecient, one would expect only small dierences in the computed o-diagonal SGS stress terms in real LES. Filtering in only one spatial direction has been assumed in these considerations, al- though similar terms are obtained upon ltering in all coordinate directions. A uniform mesh was also assumed in this expression; additional contributions to non-symmetry due to mesh non-uniformity can also be estimated with series expansions. A non-uniform mesh would add contributions from terms multiplied by dierences in adjacent cell sizes, i.e. (x + x ). Accounting for mesh nonuniformity in Eq. (7.16), it can be shown that this introduces second-order derivative terms proportional to (x + x ) 2 =144. In testing a series representation of the ltering operation that removes the eects of these terms, it was observed that their contribution to tensor non-symmetry was gener- ally quite small compared with theO x 4 terms discussed above. It is expected that the quadratic interpolation implemented to account for non-isotropy of the mesh also serves to reduce the eects of these terms. For this reason, the dierence in the o-diagonal terms is expected to be quite small and generally fourth-order in nature. Given the solver's spectral-order accuracy it is well-suited to numerically quantify the comparative dierences in the o-diagonal components of the SGS stress. The order of the non-symmetry has been quantied through the model's implementation in the numerical solver. The analytical estimate of the eects of non-symmetry suggests a small magnitude for the modeled stress imbalance. This appears to also coincide with what has been found numerically. 139 In test simulations of wall-bounded ows, planar-averaged values of the o-diagonal components of the SGS stress have been computed. This gives a description of the local introduction of angular momentum. The results have shown that the magnitude of these dierences is quite small. The planar-averaged proles of the non-vanishing o-diagonal terms (i.e. 13 and 31 ) for the case Re = 180, with resolution (32; 32; 33), are compared at a single time-slice. Due to the alignment of the mean ow direction, the ensemble averaging of the other o-diagonal terms in the SGS expression approach zero. Visually, the plotted proles of these shear stresses appear practically identical in Fig. 7.2. A more illustrative plot showing the percent dierence in the planar-averaged values of these o-diagonal terms is included in Fig. 7.3. Figure 7.2: Mean proles of the interscale transfer model's non-vanishing SGS shear stress predictions (single time-slice). 140 Figure 7.3: Interscale transfer model's predicted dierences in the non-vanishing o- diagonal values of SGS stress as a percentage of peak SGS stress (single time-slice). In these results, the dierences had a maximum of approximately 0:6% of the peak values of the predicted SGS shear stress. For comparison, in this simulation the peak values of SGS shear stress were approximately 10% of the peak resolved shear stress. In general, the dierences in the o-diagonal SGS stress values represent very small fractions of the actual peak values in stress. The dierence comparison in Fig. 7.3 is made at a single time slice and these dierences tend to diminish further upon averaging over multiple realizations. Furthermore, the deviations in the o-diagonal SGS shear stresses tend to be randomly distributed throughout the domain; thus, the eect of these dierences cancel and have no net eect on angular momentum. For these reasons the 141 eect of the SGS model's non-symmetry on the angular momentum of the resolved ow is expected to be insignicant. 7.2 Series Formulation of the Interscale Transfer Model The preceding analysis prompted the consideration of a series formulation of the proposed model. This form is obtained by implementing a series representation of the test ltering operation in the model's expression, i.e Eq. (5.10). Series representations of ltering expressions have been used previously in LES applications (e.g. [43], [31], [44], [35], [92]). This alternative to using the discrete three-point test lter oers a simple way to imple- ment the model and also allows a convenient way for terms contributing to the model's non-symmetry to be identied (dissipative or dispersive type terms in the formulation can also be identied given the order of the derivatives). It is important to note that the non- symmetry in the present model cannot be removed by simply averaging the o-diagonal components of the tensor. This gives rise to terms which violate frame invariance. Expressing the discrete lter as a 1-D operation dened as b u n i =a 0 u n i +a 1 u n1 i +u n+1 i ; (7.17) and combining the series expressions for the neighboring velocities, and recalling the normalization constraint (a 0 + 2a 1 = 1), this 1-D ltering operation takes the form b u n i =u n i +D x (u n i ): (7.18) 142 This is an expression for thex-direction ltering operation where theD x operator contains allx-derivative contributions (arising from the series expansions ofu n1 i andu n+1 i ) to the ltered term, and is dened as D x (u n i ) = a 1 (u n1 i +u n+1 i ) 2a 1 u n i = 2a 1 (x) 2 2! @ 2 u n i @x 2 + (x) 4 4! @ 4 u n i @x 4 +::: ! : (7.19) By conducting a sequence of the x, y, and z derivative operations analogous to ap- plying the three-dimensional lter kernel on a computational mesh, an approximation for the full test lter operation is obtained as b f = [(1 +D z ) (1 +D y ) (1 +D x )]f; (7.20) where similar derivative operators for the y and z directions have been introduced. Selecting a desired level of precision in the approximation of the directional derivative operators (e.g. 7.19), the above sequence of operations can give an ecient and accurate method for approximating the discrete three-point ltering operation. Applying the above formulation in the present model, the coecient a 1 in the operators above corresponds to the coecient selected in the optimized lter kernel (i.e. 1 12 ). Given that this value remains xed, the practical form of the model is easy to implement using this formulation. Also, since most all numerical ow solvers are equipped to handle derivative calculations of this type, this ltering formulation provides a convenient alternative for implementing the present model. Perhaps more importantly, this series formulation oers a clear way 143 to identify terms that contribute to non-symmetry. Non-symmetry contributions arising due to mesh non-uniformity can also be identied if one distinguishes between + and in the above series expansions. Noting that derivatives can remove terms violating frame invariance, it is convenient to represent this test lter operation as b f =f +D(f); (7.21) where theD operator represents all derivative operations that arise from the sequence of directional operators in the test lter expression above, specically D() = [(1 +D z ) (1 +D y ) (1 +D x )] 1: (7.22) For the sake of model development it is convenient to group the derivative operations together in this way to allow one to easily decompose the expression, determine the respective scaling of terms, and to identify terms giving rise to tensor non-symmetry or even Galilean invariance (if further modications to the model were considered). The discrete representation of the D operator has the form (a 1 ;2a 1 ;a 1 ). Given that the three-point lter kernel G has the relation G=I+D, or alternatively I-G=-D (where I is the identity operator), similar forms of this operator have been used in previous developments (Stolz et al. [77], Jeanmart and Winckelmans [36]) of discrete n-order lters of the form b f = [I (IG) n ]f: (7.23) 144 Recall, in this expression n is a parameter that sets the integer-valued iteration count applied in the ltering operation. This parameter sets the amplitude response dened by the lter function. The analog of this lter operation in terms of a series expansion is expressed simply as b f = [1 (D) n ]f: (7.24) It has been shown that this type of lter can give enhanced control over the lter shape. A lter of this type was not used in the present model, although it is clear that the present series approximation of the three-point lter kernel is easily adaptable for use with other lter types. Written explicitly, the comprehensive derivative operator has the form D() =D x () +D y () +D z () +D x D y () +D x D z () +D y D z () +D x D y D z (): (7.25) Also, given the relations D x () +D y () +D z () = X m a 1 x 2 m @ 2 () @x 2 m + a 1 12 x 4 m @ 4 () @x 4 m +::: ; (7.26) and 145 D x D y () +D x D z () +D y D z () = X m X n a 2 1 2 x 2 m x 2 n @ 4 () @x 2 m @x 2 n X m a 2 1 2 x 4 m @ 4 () @x 4 m +:::; (7.27) this operator has the series representation given below D () = a 1 X m x 2 m @ 2 () @x 2 m + 1 12 a 1 2 x 4 m @ 4 () @x 4 m + a 2 1 2 X m X n x 2 m x 2 n @ 4 () @x 2 m @x 2 n +::: (7.28) In these expressions summation is not implied over the repeated indices. Note that the detailed expansion of the D operator is shown here only for instruction and clarity. In practice, rather than implementing an expanded form of this operator, the application of the series formulation of the test lter can most eciently be carried out using the procedure that sequentially applies the directional derivative operators, i.e. Eq. (7.20). From the expanded form of the D operator in Eq. (7.25) it is clear that even low- order approximations for the directional derivatives can produce a nal ltering eect that contains contributions from high-order terms. High-order terms are introduced given the sequence of directional operators in D and given that the nonlinear terms in the model involve multiple applications of the test lter. It has been found that the retention of these terms is critical in order to obtain an accurate representation of the eects of the test lter. 146 Note that the lter expression in Eq. (7.21) contains only even-ordered derivative terms (small contributions from odd-ordered derivatives can be introduced when account- ing for mesh non-uniformity). However, both dissipative and dispersive terms can be generated in the representation of the model's nonlinear terms since (due to the chain rule) applying higher-order derivative operations to velocity products can produce both even and odd ordered derivatives. From the expressions above, it is clear that the application of the comprehensive derivative operator gives a representation of the small-scale eld, i.e. u i b u i = u 0 i = D(u i ): (7.29) Thus, the term N 122 has the expanded form: N 122 =D[u i D(u j )] +D[D(u i )D(u j )]: (7.30) Similarly, the term N 112 takes the form: N 112 =D [(u i +D(u i ))(u j +D(u j ))]: (7.31) Also, noting that b b u i =u i +D(u i ) +D[u i +D(u i )]; (7.32) 147 the Galilean invariant form of N 112 can be expressed as N 112 = D [(u i +D(u i ))(u j +D(u j ))] +D[u i +D(u i )]D[u j +D(u j )] +D[u i +D(u i )](u j +D(u j )) +(u i +D(u i ))D[u j +D(u j )]: (7.33) Taking the negative of the expressions forN 112 andN 122 above, and given the series representation of the D operator in Eq. (7.28), the following is obtained for the series expression of the interscale transfer model II ij = a 1 X m 2x 2 m @u i @x m @u j @x m + 1 12 a 1 2 x 4 m @ 4 (u i u j ) @x 4 m u i @ 4 u j @x 4 m @ 4 u i @x 4 m u j a 2 1 X m X n x 2 m x 2 n 2 @ 2 u i @x 2 m @ 2 u j @x 2 n 1 2 @ 4 (u i u j ) @x 2 m @x 2 n u i @ 4 u j @x 2 m @x 2 n @ 4 u i @x 2 m @x 2 n u j +u i @ 4 (u j ) @x 2 m @x 2 n 2 @ 3 (u i ) @x m @x 2 n @ (u j ) @x m +O x 6 ; (7.34) where summation over repeated indices is not taken here. In addition to being a compli- cated expression to implement, it is not obvious that retaining the above leading-order terms in the series expression is adequate for fully capturing the intended eects of the model. For simplicity, to implement the series form of the model it is instead advisable to utilize the proposed sequence of derivative operations (shown in Eq. 7.20) to replace the discrete test lter in the model's original expression (Eq. 5.10). 148 The present model has been tested using series approximations for the ltering proce- dure. It has been veried that representing terms N 122 andN 112 using a test lter con- structed from second, fourth, and sixth-order representations of the directional derivative operators gives an increasingly better match with the original model formulation. Figure 7.4 shows this comparison where mean velocity predictions obtained from these series approximations of the model are shown. The original prole obtained from the discrete ltering formulation of the present model is also shown. Figure 7.4: Mean velocity predictions using the original formulation and series represen- tations of the interscale transfer model: Re =180 (LR). If low-order approximations of the ltering operation are used it is expected that the model's predictability will be aected. However, it has been shown that these predictions 149 and those from alternative formulations of the present model can be improved by adjusting the lter strength. Conducting the ltering with the series formulation provides an alternative view for the role of the lter coecients in the proposed model. To a rst-order approximation the model is proportional to the coecienta 1 (or equivalentlya 0 ), which is used to control the lter strength. However, the higher order terms in the model have a stronger dependence on this coecient. This is in contrast to the linear relation that would be obtained if a constant of proportionality had been introduced for the sake of the model's optimization. The expanded expressions for N 112 and N 122 above allow one to identify leading- order terms and modications required for the sake of symmetrizing the model. It is the rst term in the expanded expression for N 122 , which does not permit the model's o-diagonal terms to be averaged for the sake of achieving tensor symmetry. In addition to partially replacing the production term N 122 with an advection term, averaging the o-diagonals produces a condition that violates Galilean invariance. Specically, N 122 containsu i D 2 (u j ), which satises Galilean invariance only when the divergence of veloc- ity is zero. Upon averaging the model's o-diagonals, the term D 2 (u i )u j is introduced, which is not Galilean invariant. It is of course possible to fully remove the problematic term present inN 122 ; however, initial testing has indicated that the removal of this term coupled with symmetrizing the model results in a deterioration of the model's perfor- mance. An important item of interest is that the leading-order term in the series expression for the proposed model identically matches the tensor-diusivity model (or nonlinear model) proposed by Leonard [43]: 150 nl ij =c nl 2 @u i @x m @u j @x m ; (7.35) where a value of c nl = 1 12 for this model has been proposed. In the above expression, corresponds to the mesh spacing in the direction indicated by the derivatives. It is clear that the leading-order coecients in the present model and in the tensor diusivity model dier only by a factor of two, although more recent testing by Winckelmans et al. [92] has shown c nl 1. Variants of this model have also been considered by Clark et al. [11] and more recently by Liu et al. [48], Vreman et al. [86], Leonard [44], and Winckelmans et al. [92]. This model has given reasonable results in testing but leads to instability issues due to its eventual prediction of negative viscosity (Chumakov [10], Winckelmans et al. [92]). Like the similarity model, the tensor-diusivity model also suers from an inability to properly account for aspects of energy transfer. In an attempt to address this issue, Bardina et al. [2] and others (e.g. [85], [86], [91], [45], [92]) have altered the formulation of this model to include eddy-viscosity terms (i.e. mixed models). It should be noted that an expansion of the similarity model also shares the same leading-order term as that included in the tensor diusivity model [85]. This is due to the similarity model's incorporation of N 112 . In theory, the tensor-diusivity model amounts to a series-expansion alternative to the similarity model. However, the tensor diusivity model's incorporation of the leading-order term alone does not provide an adequate amount of detail needed to reproduce the eects of the similarity model's test- ltered quantities. This is unimportant given that the similarity model tends to fail 151 regardless of its implementation, so practically speaking, the two models can be thought of as being independent. In the present research it was veried that making modications to the nonlinear model's coecient and the similarity model's coecient and test lter strength could not emulate the important role of the test lter in the interscale transfer model. For the sake of the present model, retaining the formal role of the test lter is critical to maintaining the model's intended eect on energy transfer. Thus, in the proposed series formulation, particular attention has been paid to the expanded form of the lter and its use in constructing the model's partial nonlinear terms. In the interest of simplicity, the subtle but important point to be made is that choosing an Nth-order approximation of the directional derivative operators and then applying the constructed test lter accord- ing to the model's expression, is superior to merely representing the model using series expansions and retaining its terms up through order N. These two approaches can give considerably dierent levels of accuracy, the former being more accurate and quite simple to implement. Clearly, expressing the interscale transfer model as a series expansion provides ad- ditional options for implementing the model and/or altering its formulation. The series formulation of this model is likely more amenable to a broader range of solvers, and also more readily lends itself to implementation on unstructured grids. Future eorts to achieve tensor symmetry will remain under consideration; however, given the relatively small magnitude of the non-symmetry and its minor role in the model's performance, this issue is not presently an item of major concern. 152 7.3 SGS Model Computational Expense Among the limiting factors in applications of computational physics is the computing expense associated with modeling the governing equations and properly resolving a prob- lem's important features. During a solver's development stage, the modeling approach and the type of numerical algorithms to be used are among the items that must be con- sidered with regard to the desired level of accuracy and program eciency. The ratios of LES to DNS resolution given in Table 6.1 show that the mesh resolution (and therefore computational work) of LES is generally a very small fraction of what is required in di- rect numerical simulations. In LES, both the accuracy and eciency of a SGS model are among the items that must be considered. It has been shown that the interscale transfer model can predict statistical quantities closely matching those obtained from DNS. Esti- mates of the computing expense associated with using this model are considered in the present section. In regard to run-time eciency, most numerical applications are limited by various factors including memory allocation, data access, and algorithm eciency. Many mod- ern CFD applications are carried out across multiple compute nodes so scalability and processor communication are also important. In CFD applications, computationally in- tensive models can greatly aect run-time eciency. For this reason, it is instructive to quantify the computational load associated with the proposed interscale transfer model and compare this to other more established SGS models in the eld. Basic integer or oating-point operations such as addition and subtraction are carried out by a computer's central processing unit. A computer's performance is generally 153 measured in terms of its ability to complete specic oating-point operations over a block of time. Likewise, a program's eciency can be measured by the number of equivalent oating-point operations required to complete a specic task (see Chapra and Canale [9]). Based on estimates of required compute time, calculations of various mathematical functions can be expressed in terms of an equivalent oating-point operation count. By using a set of equivalent operations for these basic functions, a standardized estimate of the computational expense associated with the various SGS models can be obtained. This interpretation of computing expense allows one to estimate the relative compute time, or the comparative computational load, associated with using these models in LES. Estimates of the equivalent oating-point operations supplied by Schorghofer [71], and used in the present analysis, are given in Table 7.1. integer addition/subtraction/multiplication 1 oat addition/subtraction/multiplication 1 integer division 8 oat division 4 square root 13 exponent/logarithm/trigonometric functions 25 Table 7.1: Floating-point operation assignments for various mathematical operations. Given that the computational work associated with calculating derivatives (required in the Smagorinsky, dynamic, and the series implementation of the present model) can vary greatly with solver type, in this study the operation count for taking derivatives was associated with fourth-order accurate nite-dierence expressions and was computed us- ing the listed operation assignments. This simplication gives more generalized estimates of the equivalent number of oating-point operations per mesh point required for these models. 154 The total operation count for each SGS model has been computed assuming an ecient implementation of the model. Assuming the subroutine for each SGS model is called only once per time step, the relative computational expense at each time step can be estimated using the model's total operation count and number of cells in the computational mesh. Approximations of the oating-point operation count per mesh point for the similarity model, Smagorinsky model, dynamic model, and the original and series formulations of the interscale transfer model are given in Table 7.2. Similarity Model 230 Smagorinsky Model 260 Dynamic Model 1280 Interscale Transfer Model 360 Interscale Transfer Model (Series Form): O(x 4 ) 1240 Table 7.2: Estimates of oating-point operation count per mesh point for various SGS models. These estimates provide a relative picture of the computational expense associated with each SGS model. Of course, not accounted for in these estimates are items such as the model's memory utilization and the communication required between processors. For instance, it is known that the planar-averaging operation implemented by the dynamic model can require a considerable amount of information transfer across computational subdomains. This communication time is not accounted for in the above estimates. However, given that each model has the same general purpose within the solver and utilizes similar information stored within the primitive variables, concerns about memory and inter-processor communication are assumed to be comparable and have not been addressed here. 155 The values in Table 7.2 give standardized estimates of operation count at each point on an LES mesh. Thus, the total computational expense can be estimated by taking the product of these values with the mesh size. These results show that the similarity model, the Smagorinsky model, and the discrete formulation of the interscale transfer model are generally the least computationally intensive options. The similarity model is only slightly more ecient than the Smagorinsky model and both of these are somewhat more ecient than the interscale transfer model. Compared with the present model, the similarity and Smagorinsky models represent about a 35% decrease in computational expense depending on whether mesh interpolation is used. The additional expense of the present model is generally quite small given this model's improvement in performance relative to the similarity and Smagorinsky models. Given that the interpolation implemented in the discrete ltering procedure used in the similarity and interscale transfer models is optional and in some cases should be eliminated for eciency, the operation count for these two models can be reduced from 230 and 360 to approximately 160 and 270, respectively. In this case, this puts the computational expense of the interscale transfer model on par with the Smagorinsky model while the similarity model remains somewhat less than both of these models. In general, the computational expense of the Smagorinsky model is considered to be quite low compared with the dynamic model, which was the most expensive SGS model tested. Despite the dynamic model's additional expense relative to other options, it remains a popular SGS model in LES. Using derivative operations associated with fourth-order nite dierence expressions, the computational expense associated with the 156 series formulation of the interscale transfer model is slightly less than that of the dy- namic model. This also holds if sixth-order dierence expressions are used. In the series formulation of the present model and in the dynamic model, somewhat more ecient formulations of both of these models could feasibly be constructed by utilizing velocity gradient information already computed in the viscous term in the governing equations. By today's standards a computing cluster can complete on the order of 510 9 oating- point operations per processor, per second [34]. Using this estimate of performance, the compute time associated with each SGS model's computational load can be estimated. For instance, an implementation of the dynamic model in the Re = 950 case having resolution (512, 384, 85), would ideally require about 4.3 seconds per time step (1280 512 384 85=5 10 9 4:3), while the interscale transfer model would require less than one-third of this time. This savings is not trivial given that many thousands of time steps are required to reach a statistically steady-state solution. Given that the workload can be distributed over multiple machines, this compute time, however, can be reduced considerably. This estimate is of course associated with the ideal compute time of the SGS term alone; typically the other terms treated by the numerical solver account for the greatest fraction of the computational expense. The actual compute time varies with the computing hardware, code eciency, and the distribution of the computational load across multiple processors. However, regardless of the hardware and parallelization, the general trends in terms of relative computing expense between SGS models should remain somewhat consistent with the values given above. 157 Chapter 8 Conclusions The under-resolved nature of LES reduces the dissipative role of molecular viscosity that acts largely on the smallest scales of turbulence. The absence of this dissipative mechanism in LES can lead to numerical instabilities and generally poor predictions of actual ow quantities. A feature of SGS models that is of great importance is the ability to properly account for energy dissipation of the subgrid scales and the resulting energy ux from the resolved to the unresolved scales that is absent in LES. Generally, a model that properly accounts for this key deciency also provides good predictions of mean velocity and RMS values of velocity uctuations. A model's predictions of the mean velocity prole in a boundary layer is also key to predicting such things as wall shear stress, skin friction, etc. Eddy-viscosity models are known to account for the dissipative nature of the small scales of turbulence and the majority of SGS models adopt some naturally dissipative relationship to account for the absence of adequate dissipation in large-eddy simulation. The family of SGS models referred to as scale-similarity models use a secondary lter to construct a model with the assertion that the smallest resolved scales of turbulence 158 in LES can be used to approximate the stresses of the unresolved scales. This model employs arguments based on similarities of stresses generated by velocity components having neighboring scales. In previous applications and testing, similarity models have been unable to correctly capture the important dissipative nature of the small unresolved scales of turbulence. The present research has shown that the majority of the similarity model's SGS dis- sipation is associated with the largest resolved scales. The similarity model's inadequate prediction of SGS dissipation in the smallest scales has generally been addressed indi- rectly with the addition of ad-hoc dissipative terms. The central goal of this research has been to formulate a new SGS model that addresses the deciencies of the similarity model directly. The proposed model uses a similarity-type expression that accounts for the dissipative nature of the subgrid scales of turbulence. The present interscale transfer model has been formulated using concepts of interscale energy transfer, which are rooted in spectral space. These concepts oer a tool to parse the nonlinear transfer present in LES and formulate expressions capable of removing (or reducing) terms that deposit en- ergy near the LES cuto. The present formulation is such that a fraction of the nonlinear terms leading to production in the the smallest resolved scales of turbulence are removed by the model. The strength of the test lter has been used to dene the proportion of these terms removed by the model. The lter strength combined with implicit factors such as mesh resolution and Reynolds number determine the model's collective eect on energy transfer. Selecting a smooth lter allows for some energy removal at larger scales and relaxes the requirement that energy rst be passed to the smallest scales before SGS 159 dissipation takes eect. Tuning of the test lter strength was required initially and numer- ical simulations at various Reynolds numbers and mesh resolutions have demonstrated that good performance can be obtained without additional changes to the lter. In testing, the model has continued to provide good mean and RMS velocity predic- tions using this identical lter setting. The model's eectiveness over a range of Reynolds numbers, each having some unique level of scale separation, suggests a exibility and robustness important for more general applications in LES. The model's consistent pre- dictions at dierent wall-unit resolutions is also reassuring given the model's intrinsic sensitivity to the lter setting and the inevitable link between resolution and the test lter cuto. Modeling over various Reynolds numbers illustrates its ability to capture the mean ow activity with changes in the turbulent scales prevalent throughout the boundary layer. The proposed model has demonstrated considerable improvements when compared with predictions of the similarity model. Through the considerations of interscale energy transfer, it has been shown that the present model is comparatively better suited to account for the dissipative eects of the unresolved scales of turbulence. The present model emulates the role of the exact form of the SGS stress tensor by encouraging the most prominent exchange process to occur from the largest to the smallest scales of turbulence. The model has also demonstrated a predictive capability that rivals (and often surpasses) that of the dynamic model without suering from the same drawbacks such as the dynamic model's elimination of backscatter or its limiting requirements of spatial or ensemble averaging. 160 The dynamic model has an inherent sensitivity to backscatter that forces special treatment of the modeled expression, in this case averaging along homogeneous planes. This spatial averaging limits the applicability of the model for more complex geometries. On the other hand, the proposed model is relatively simple in its formulation and lacks any obvious drawbacks in terms of its general applicability. The proposed model has given predictions comparable to the dynamic model's in a series of LES simulations of wall-bounded turbulent ows. Future research is likely to be conducted with regard to this model's performance under dierent ow conditions. A SGS model's direct role in resolved energy transfer is limited to SGS dissipation and SGS diusion. It is known that the present model's non-symmetry cannot produce a net eect on the energy content in LES given that SGS dissipation is computed as a contraction with the SGS stress tensor and the (symmetric) resolved strain-rate tensor, and also given that SGS diusion appears in ux-divergence form in the resolved energy equation and is conservative in LES. Also, given that the dierences in the model's o- diagonal terms can be expressed as dierences in ltered and non-ltered quantities, the model's non-symmetry is limited to having only local eects on the ow's angular momentum. Given the convolution properties of the test ltering, no net eect on angular momentum is produced by these dierences. The non-symmetry of the present model's SGS stress has been considered both ana- lytically and numerically. These considerations have shown that the model's local eects due to non-symmetry are small relative to the model's leading-order terms. It has been shown analytically that the order of the non-symmetry is best approximated as being O( 4 ) and has been characterized as having a minor role relative to the leading-order 161 stress components produced by the model itself. In the test simulations of wall-bounded ows, the planar-averaged values of the o-diagonal components of the SGS stress were computed. The results showed that the magnitude of these dierences was quite small. In these results, the dierences had maximum values that represented a fraction of the peak value of the predicted SGS shear stress. Thus, this eect is quite small relative to the resolved scale and SGS stress terms. In comparing values of the planar-averaged, non-vanishing, SGS shear stresses, the proles of these plots appear practically identical. Furthermore, the minor deviations tended to be randomly distributed throughout the domain and the resulting eects are expected to be negligible. A series representation of the model was derived by expanding the model's nonlinear terms. This form is obtained by representing the model's discrete test ltering operation as a sequence of operations involving second and higher order derivative operators in each coordinate direction. The series formulation oers an alternative form for the present model where the role of the incorporated terms can be identied more directly. This formulation of the model has been implemented and the results obtained from the model's series representation compared well with the original formulation. It is possible that this new formulation will be a more desirable form of the model given the ease of implementing it in unstructured solvers and in complex geometries. Finally, estimates of the computational expense of various SGS models were given and these were compared with similar estimates for the discrete and series formulations of the proposed interscale transfer model. The discrete formulation of the present model is only slightly more costly than the similarity and Smagorinsky models, but is considerably less computationally intensive than the dynamic model. The series formulation of the present 162 model is slightly less than the dynamic model in terms of its computational expense. Relative to the other SGS models considered, the computational expense of the present model is quite advantageous given this model's improvements in performance. 163 References [1] B. S. Baldwin and H. Lomax. Thin-layer approximation and algebraic model for separated turbulent ows. AIAA J., 16, 78-257, 1978. [2] J. Bardina, J. H. Ferziger, and W. C. Reynolds. Improved subgrid scale models for large eddy simulation. AIAA J., 80-1357, 1980. [3] J. Bardina, J. H. Ferziger, and W. C. Reynolds. Improved turbulence models based on large eddy simulation of homogeneous, incompressible, turbulent ows, volume TF-19 of Stanford University Technical Report. 1983. [4] G. K. Batchelor. The theory of homogeneous turbulence, volume 1. Cambridge University Press, 1953. [5] J. Boussinesq. Theorie d l'ecoulelement tourbillant, volume 1. Mem. Pres., 1877. [6] D. Carati, K. Jansen, and T. Lund. A family of dynamic models for large eddy simulation, volume 237-48 of Center for Turbulence Research, Proceedings of the Summer Program. 1995. [7] D. Carati, A. Wray, and W. Cabot. Ensemble averaged dynamic modeling, volume 237-48 of Center for Turbulence Research, Proceedings of the Summer Program. 1996. [8] T. Cebeci. Turbulence models and their application: ecient numerical methods with computer programs. Springer-Verlag, 2004. [9] C. Chapra and R. Canale. Numerical Methods for Engineers. McGraw-Hill, 3rd edition, 1998. [10] S. G. Chumakov. Statistics of subgrid-scale stress states in homogeneous isotropic turbulence. J. Fluid Mech., 562:405{414, 2006. [11] R. A. Clark, J. H. Ferziger, and W. C. Reynolds. Evaluation of subgrid-scale models using an accurately simulated turbulent ow. J. Fluid Mech., 91:1{16, 1979. [12] A. W. Cook. Determination of the constant coecient in scale similarity models of turbulence. Phys. Fluids, 9:1485{1487, 1997. [13] J. W. Deardor. A numerical study of three-dimensional turbulent channel ow at large Reynolds numbers. J. Fluid Mech., 41:453{480, 1970. 164 [14] J. W. Deardor. Stratocumulus-capped mixed layers derived from a three- dimensional model. Boundary-Layer Meteorology, 18:495{527, 1980. [15] J. C. Del Alamo and J. Jimenez. Direct numerical simulation of the very large anisotropic scales in a turbulent channel. Center for Turbulence Research, Annual Research Briefs, pages 329{341, 2001. [16] J. C. Del Alamo and J. Jimenez. Spectra of the very large anisotropic scales in turbulent channels. Phys. Fluids, 15, 6:41{44, 2003. [17] J. C. Del Alamo, J. Jimenez, P. Zandonade, and R. D. Moser. Scaling of the energy spectra of turbulent channels. J. Fluid Mech., 500:135{144, 2004. [18] P. J. Diamessis, J. A. Domaradzki, and J. S. Hesthaven. A spectral multidomain penalty method model for the simulation of high Reynolds number localized incom- pressible stratied turbulence. J. Comp. Phys., 202:298{322, 2005. [19] J. A. Domaradzki. Large eddy simulations without explicit eddy viscosity models. Int. J. of Comput. Fluid D., 24 10:435{447, 2010. [20] J. A. Domaradzki and D. Carati. A comparison of spectral sharp and smooth lters in the analysis of nonlinear interactions and energy transfer in turbulence. Phys. Fluids, 19, 085111, 2007. [21] J. A. Domaradzki and D. D. Holm. Navier-Stokes-alpha model: LES equations with nonlinear dispersion. Modern Simulation Strategies for Turbulent Flow. R.T. Ed- wards, Inc., 2001. [22] J. A. Domaradzki and W. Liu. Approximation of subgrid-scale energy transfer based on the dynamics of resolved scales of turbulence. Phys. Fluids, 7:2025{2035, 1995. [23] J. A. Domaradzki, W. Liu, and M. E. Brachet. An analysis of subgrid-scale inter- actions in numerically simulated isotropic turbulence. Phys. Fluids, A 5:1747{1759, 1993. [24] J. A. Domaradzki, W. Liu, Hartel C., and L. Kleiser. Energy transfer in numerically simulated wall-bounded ows. Phys. Fluids, A 6:1583{1599, 1994. [25] J. A. Domaradzki, K. Loh, and P. P. Yee. Large eddy simulations using the subgrid- scale estimation model and truncated Navier-Stokes dynamics. Theor. Comput. Fluid Dyn., 15:421{450, 2002. [26] J. A. Domaradzki and R. S. Rogallo. Energy transfer in isotropic turbulence at low Reynolds numbers. Center for Turbulence Research, Proceedings of the Summer Program. 1988. [27] J. A. Domaradzki and R. S. Rogallo. Local energy transfer and nonlocal interactions in homogeneous, isotropic turbulence. Phys. Fluids, A 2:413{426, 1990. 165 [28] J. A. Domaradzki, R. S. Rogallo, and A. A. Wray. Interscale energy transfer in numerically simulated homogeneous turbulence. Center for Turbulence Research, Proceedings of the Summer Program. 1990. [29] J. A. Domaradzki and E. M. Saiki. A subgrid-scale model based on the estimation of unresolved scales of turbulence. Phys. Fluids, 9:2148{2164, 1997. [30] C. A. J. Fletcher. Computational Techniques for Fluid Dynamics, volume 1. Springer-Verlag, 2nd edition, 1991. [31] M. Germano. Dierential lters for the large eddy simulation of turbulence. Phys. Fluids, 29:1755{1758, 1986. [32] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgrid-scale eddy viscosity model. 3:1760{1765, 1991. [33] B. J. Geurts and D. D. Holm. Regularization modeling for large-eddy simulation. Phys. Fluids, 15:13{16, 2003. [34] W. Group. USC Center for High-Performance Computing and Communications. http:==www:usc:edu=hpcc=images=HPCC brochure 2011 FINAL:pdf, 2011. [35] K. Horiuti. A new dynamic two-parameter mixed model for large-eddy simulation. Phys. Fluids, 29:3443{3464, 1997. [36] H. Jeanmart and G. S. Winckelmans. Investigation of eddy-viscosity models modied using discrete lters: a simplied "regularized variational multiscale model" and an "enhanced eld model". Phys. Fluids, 19, 2007. [37] W. P. Jones and B. E. Launder. The prediction of laminarization with a two-equation model of turbulence. International Journal of Heat and Mass Transfer, 15:301{314, 1972. [38] A. N. Kolmogorov. Dissipation of energy in locally isotropic turbulence. Dokl. Akad. Nauk SSSR, 32:16{18, 1941. [39] A. N. Kolmogorov. Equations of turbulent motion of an incompressible uid. Izvestia Academy of Science, USSR, Physics, 6:56{58, 1942. [40] N. V. Kornev, I. V. Tkatchenko, and E. Hassel. A simple clipping procedure for the dynamic mixed model based on Taylor series approximation. Commun. Numer. Meth. Engng., 22:55{61, 2006. [41] P. K. Kundu. Fluid Mechanics. Academic Press, 1990. [42] B. E. Launder, G. J. Reece, and W. Rodi. Progress in the development of a Reynolds- stress turbulence closure. J. Fluid Mech., 68:537{566, 1975. [43] A. Leonard. Energy cascade in large-eddy simulations of turbulent uid ows. Adv. Gophys., 18:237, 1974. 166 [44] A. Leonard. Large-eddy simulation of chaotic convection and beyond. AIAA J., 97-0204, 1997. [45] A. Leonard and G. Winckelmans. A tensor-diusivity subgrid model for large-eddy simulation. CalTech Technical Report (TF-43). 1999. [46] M. Lesieur. Turbulence in uids. Kluwer Academic Pub., 1990. [47] D. K. Lilly. A proposed modication of the Germano subgrid-scale closure. Phys. Fluids, A 4 (3), 1992. [48] S. Liu, C. Meneveau, and J. Katz. On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet. J. Fluid Mech., 275:83{ 119, 1994. [49] K. Loh and J. A. Domaradzki. The subgrid-scale estimation model on nonuniform grids. Phys. Fluids, 11:3786{3792, 1999. [50] T. S. Lund. On the use of discrete lters for large-eddy simulation. Center for Turbulence Research, Annual Research Briefs, pages 83{95, 1997. [51] N. N. Mansour, J. H. Ferziger, and Reynolds W. C. Large-eddy simulation of a tur- bulent mixing layer. Thermosciences Div., Dept. of Mech. Eng., Stanford University, TF-11, 1978. [52] J. Mathieu and J. Scott. An Introduction to Turbulent Flow. Cambridge University Press, 2000. [53] O. McMillan and J. Ferziger. Direct testing of subgrid-scale models. AIAA J., 17(12):1340{1346, 1979. [54] C. Meneveau and J. Katz. Scale-invariance and turbulence models for large-eddy simulation. Annu. Rev. Fluid Mech., 32:1{32, 2000. [55] C. Meneveau, T. S. Lund, and W. H. Cabot. A Lagrangian dynamic subgrid-scale model of turbulence. J. Fluid Mech., 319:353{385, 1996. [56] C. Meneveau and T.S. Lund. Search for subgrid-scale parameterization by projection pursuit regression. Center for Turbulence Research, Proceedings of the Summer Program. 1992. [57] P. Moin and J. Kim. Numerical investigation of turbulent channel ow. J. Fluid Mech., 118:341{377, 1982. [58] J. O'Neil and C. Meneveau. Subgrid-scale stresses and their modeling in a turbulent plane wake. J. Fluid Mech., 349:253{93, 1997. [59] R. Panton. Incompressible Flow. Wiley, 1996. [60] N. Park, S. Lee, J. Lee, and H. Cho. A dynamic subgrid-scale eddy viscosity model with a global model coecient. Phys. Fluids, 18, 2006. 167 [61] N. Park and K. Mahesh. Reduction of the Germano-identity error in the dynamic Smagorinsky model. Phys. Fluids, 21, 065106, 2009. [62] U. Piomelli. Large-eddy and direct simulation of turbulent ows. Short course: 9e conference annuelle de la Societe canadienne de CFD, 2001. [63] U. Piomelli. Wall-layer models for large-eddy simulations. Prog. in Aerosp. Sci., 44,6:437{446, 2008. [64] U. Piomelli, P. Moin, and J. H. Ferziger. Model consistency in large eddy simulation of turbulent channel ow. Phys. Fluids, 31, 1988. [65] S. Pope. Turbulent Flows. Cambridge University Press, 2000. [66] L. Prandtl. Uber die ausgebildete turbulenz. Z. Angew. Math. Mech., 5, 136-139 1925. [67] L. Prandtl. Uber ein neues Formelsystem fur die ausgebildete turbulenz. Nacr. Akad. Wiss. Gottingen, Math-Phys., 6-19 1945. [68] L. F. Richardson. Weather Prediction by Numerical Process. Cambridge University Press, 1922. [69] P. Sagaut. Subgrid-scale modeling issues and approaches. in Implicit large-eddy simulation: computing turbulent uid dynamics. Cambridge University Press, 2007. [70] H. Schlichting. Boundary-Layer Theory. McGraw-Hill, 6th edition, 1968. [71] N. Schorghofer. Lessons in Computational Physics. http://www2.hawaii.edu/ nor- bert/compphysics.html, 2012. [72] U. Schumann. Subgrid scale model for nite dierence simulations of turbulent ows in plane channels and annuli. J. Comp. Phys., 18:376{404, 1975. [73] A. Scotti, C. Meneveau, and D. K. Lilly. Generalized Smagorinsky model for anisotropic grids. Phys. Fluids, 5:2306{2308, 1993. [74] J. Smagorinsky. General circulation experiments with the primitive equations. i. the basic experiment. Monthly Weather Review, 91, 3:99{164, 1963. [75] P. R. Spalart and S. R. Allmaras. A one-equation turbulence model for aerodynamic ows. AIAA J., 1992. [76] C. G. Speziale. Galilean invariance of subgrid-scale stress models in the large eddy simulation of turbulence. J. Fluid Mech., 156:55{62, 1985. [77] S. Stolz, N. A. Adams, and L. Kleiser. An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded ows. Phys. Fluids, 13, 4:997{1015, 2001. 168 [78] T. Tantikul and J. A. Domaradzki. Large eddy simulations using truncated Navier- Stokes equations with the automatic ltering criterion. Journal of Turbulence, 11, 2010. [79] T. Tantikul and J. A. Domaradzki. Truncated Navier-Stokes equations with the automatic ltering criterion: Reynolds stress and energy budgets. Journal of Tur- bulence, 12, 2011. [80] H. Tennekes and J. L. Lumley. A First Course in Turbulence. MIT Press, 1972. [81] E. R. Van Driest. On the turbulent ow near a wall. J. Aero. Sci., 23, 1956. [82] A. W. Vreman. The ltering analog of the variational multiscale method in large- eddy simulation. Phys. Fluids, 8:61{64, 2003. [83] A. W. Vreman. An eddy-viscosity subgrid-scale model for turbulent shear ow: algebraic theory and applications. Phys. Fluids, 16:3670{3681, 2004. [84] B. Vreman, B. Geurts, and H. Kuerten. On the formulation of the dynamic mixed subgrid-scale model. Phys. Fluids, 6:4057{4059, 1994. [85] B. Vreman, B. Geurts, and H. Kuerten. Large eddy simulation of the temporal mixing layer using the Clark model. Theor. Comput. Fluid Dyn., 8:309{324, 1996. [86] B. Vreman, B. Geurts, and H. Kuerten. Large-eddy simulation of the turbulent mixing layer. J. Fluid Mech., 339:357{390, 1997. [87] F. M. White. Viscous Fluid Flow. McGraw-Hill, 2nd edition, 1991. [88] D. C. Wilcox. Turbulence Modeling for CFD. La Canada: DCW Industries, 3rd edition, 2006. [89] G. S. Winckelmans. Some progress in large-eddy simulation using the 3-D vortex par- ticle method. Center for Turbulence Research, Proceedings of the Summer Program. 1995. [90] G. S. Winckelmans, T. S. Lund, D. Carati, and A. A. Wray. A priori testing of subgrid-scale models for the velocity-pressure and vorticity-velocity fomulations. Cen- ter for Turbulence Research, Proceedings of the Summer Program. 309-328, 1996. [91] G. S. Winckelmans, A. A. Wray, and O. V. Vasilyev. Testing of a new mixed model for LES: The Leonard model supplemented by a dynamic Smagorinsky term. Center for Turbulence Research, Proceedings of the Summer Program. 367-388, 1998. [92] G. S. Winckelmans, A. A. Wray, O. V. Vasilyev, and H. Jeanmart. Explicit-ltering large-eddy simulation using the tensor diusivity model supplemented by a dynamic Smagorinsky term. Phys. Fluids, 13, 5:1385{1403, 2001. [93] P. K. Yeung and J. G. Brasseur. The response of isotropic turbulence to isotropic and anisotropic forcing at the large scales. Phys. Fluids, 3:884{897, 1991. 169 [94] Y. Zang, R. L. Street, and J. R. Kose. A dynamic mixed subgrid-scale model and its application to turbulent recirculating ow. Phys. Fluids, 5:3186{3196, 1993. [95] Y. Zhou. Degrees of locality of energy transfer in the inertial range. Phys. Fluids, page 1092, 1993. 170

Abstract (if available)

Abstract In large-eddy simulation (LES) various subgrid-scale models have been proposed, all of which attempt to account for the unknown effects of the unresolved scales of turbulence on the resolved flow-field. The scale-similarity model is one such model, which is formulated using a secondary filter applied to components of the resolved velocity and its products. The scale-similarity model is based on the assumption that the velocities associated with neighboring scales in the flow produce turbulent stresses that are similar in character. The model uses an expression that weights the scales just below the LES cutoff in its approximation of the stresses of the unresolved scales. It is well-known, however, that the similarity model fails to accurately predict some of the most fundamental quantities in turbulent flows, perhaps the most important being the global energy transfer and the associated subgrid-scale dissipation. In previous research, additional dissipative terms have been added to the similarity model to improve the model's performance. ❧ In the present research, considerations of interscale energy transfer have been used to identify the deficiencies of the energy transfer role of the similarity model, specifically its inadequate removal of terms contributing energy to the smallest scales and its duplication of terms producing effects in the largest scales. These considerations are then used in the development of a new model, which shows more favorable characteristics of energy transfer. In this approach, partial nonlinear terms are used to decompose the nonlinear transfer present in LES and to formulate an expression capable of removing small-scale production terms depositing energy near the LES cutoff. The proposed model is formulated in the same vein as the scale-similarity model, consisting of test-filtered velocities and their products, but the new interscale transfer model offers improvements in its predictions of mean flow quantities and the global energy flux from the resolved to subgrid scales. The current research demonstrates that by implementing this model in a posteriori LES testing of wall-bounded flows, improved LES predictions are possible without the need for additional terms to augment subgrid-scale energy dissipation.

Conceptually similar

PDF

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

PDF

Large eddy simulations of laminar separation bubble flows

PDF

On the simulation of stratified turbulent flows

PDF

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

PDF

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

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Modeling and simulation of complex recovery processes

PDF

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

PDF

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

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

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

PDF

Nonlinear models for hippocampal synapses for multi-scale modeling

PDF

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

PDF

High temperature latent heat thermal energy storage to augment solar thermal propulsion for microsatellites

PDF

Point cloud data fusion of RGB and thermal information for advanced building envelope modeling in support of energy audits for large districts

PDF

A parametric study of the thermal performance of green roofs in different climates through energy modeling

PDF

Investigation of gas transport and sorption in shales

PDF

The effect of lattice structure and porosity on thermal conductivity of additively-manufactured porous materials

Asset Metadata

Creator Anderson, Brian Wayne (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 04/16/2012

Defense Date 03/30/2012

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

Tag differential filter,discrete filter,interscale energy,interscale energy transfer,interscale energy transfer model,interscale transfer,interscale transfer model,interscale-energy transfer,itm,large-eddy,large-eddy simulation,LES,OAI-PMH Harvest,scale-similarity model,similarity model,subfilter scale,subgrid,subgrid scale,subgrid scale dissipation,subgrid scale model,subgrid-scale,subgrid-scale model

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Domaradzki, Julian A. (committee chair), Blackwelder, Ron F. (committee member), Pottebaum, Tait S. (committee member), Shiflett, Geoffrey R. (committee member), Wellford, L. Carter (committee member), Wilcox, David C. (committee member)

Creator Email bwanders@usc.edu,bwapsu@hotmail.com

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

Unique identifier UC1113017

Identifier usctheses-c3-6676 (legacy record id)

Legacy Identifier etd-AndersonBr-610.pdf

Dmrecord 6676

Document Type Dissertation

Rights Anderson, Brian Wayne

Type texts

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

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

Repository Name University of Southern California Digital Library

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