ABSTRACT Title of dissertation: HYBRID RANS LES SIMULATION OF NON-EQUILIBRIUM BOUNDARY LAYERS Giuseppe De Prisco Doctor of Philosophy, 2007 Dissertation directed by: Professor Ugo Piomelli Department of Mechanical Engineering Hybrid simulations that couple the solution of the Reynolds-Averaged Navier- Stokes equations (RANS) to Large-Eddy Simulations (LES) have the ability to apply the high accuracy of LES only in regions of the ow that demand it, while using the less expensive RANS approach in regions of the ow where standard turbulence models are expected to be accurate, while LES is used in non-equilibrium ow regions. One issue that arises in these applications is the behavior of the ow in the transition zone between the RANS and LES regions. In the RANS zone the ow solution is either steady, or only contains information on the largest scales of motion; most or all of the Reynolds shear stress is provided by the turbulence model. In the LES region the resolved scales must supply most of the Reynolds shear stress. Typically, a transition zone exists in which the resolved eddies are gradually generated and grow. In this work, methodologies for the improvement of hybrid LES/RANS are studied, to shorten the transition from the smooth RANS eld to the LES, which requires energy- and momentum-supporting eddies. The method tested is based on the generation of synthetic turbulence with a realistic spectrum, and statistics obtained from the RANS. The eddies thus generated are then selectively forced by a control method that ampli es the bursts, and maintains a desired Reynolds shear stress downstream of the RANS/LES interface. This method allows to match the two techniques smoothly, and to minimize the extent of the region required to develop the realistic turbulent eddies. It was found to perform well in several non-equilibrium boundary layers achieved by imposing either a variable freestream velocity or a spanwise pressure gradient on a at-plate boundary layer. A nely resolved LES of the ow in a boundary layer subjected to strong acceleration was also performed. The ow in this con guration reverts to a laminar state and then retransitions to turbulence. Statistics and ow visualizations of this ow indicate the presence of two of the mechanisms that have been conjectured to cause the re-laminarization. HYBRID RANS LES SIMULATION OF NON-EQUILIBRIUM BOUNDARY LAYERS by Giuseppe De Prisco Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosophy 2007 Advisory Committee: Professor Ugo Piomelli, Chair/Advisor Professor Richard Calabrese, Professor James Duncan Professor Kenneth Kiger Professor Elias Balaras c Copyright by Giuseppe De Prisco 2007 Acknowledgements My special thanks to Professors Ugo Piomelli and Elias Balaras and Dr. Tony Keating for their guidance. Financial support of the Air Force O ce of Scienti c Research is gratefully acknowledged. E ort sponsored by the Air Force O ce of Scienti c Research, under grant number FA95500610116, monitored by Dr. John D. Schmisseur. The U.S. Govern- ment is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions con- tained herein are those of the authors and should not be interpreted as necessarily representing the o cial policies or endorsements, either expressed or implied, of the Air Force O ce of Scienti c Research or the U.S. Government. ii Table of Contents List of Figures v List of Abbreviations x 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Hybrid strategies: single grid and zonal approaches . . . . . . . . . . 4 1.3 Hybrid RANS/LES approach with single grid . . . . . . . . . . . . . 8 1.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2 Previous work . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Zonal Hybrid RANS/LES . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4.2 In ow generation methods . . . . . . . . . . . . . . . . . . . . 18 1.4.3 Previous work . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.5 Plan of the present dissertation . . . . . . . . . . . . . . . . . . . . . 25 2 Problem formulation 27 2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Numerical method and boundary conditions . . . . . . . . . . . . . . 29 2.3 Subgrid-scale modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3.1 Eddy-viscosity assumption . . . . . . . . . . . . . . . . . . . . 33 2.3.2 The dynamic procedure . . . . . . . . . . . . . . . . . . . . . 34 2.3.3 Averaging of the eddy-viscosity coe cient . . . . . . . . . . . 36 2.3.4 Filtering operations . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4 RANS equations and Reynolds stress modeling . . . . . . . . . . . . . 38 2.4.1 Spalart-Allmaras model . . . . . . . . . . . . . . . . . . . . . 39 2.4.2 k model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.5 Synthetic turbulence generation . . . . . . . . . . . . . . . . . . . . . 42 3 The controlled forcing method and its tuning 50 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Analysis of the controller parameters . . . . . . . . . . . . . . . . . . 53 3.2.1 Testing strategy: ZPG boundary layer . . . . . . . . . . . . . 53 3.2.2 The moving-average exponential window . . . . . . . . . . . . 55 3.2.3 The PI control . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 iii 4 Application of the controlled forcing 68 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.1.1 Zero-pressure-gradient boundary layer . . . . . . . . . . . . . 68 4.1.2 Adverse-pressure-gradient boundary layer . . . . . . . . . . . . 74 4.1.3 Favorable-pressure-gradient boundary layer . . . . . . . . . . . 77 4.1.4 3-D Boundary layer . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5 Favorable pressure gradient boundary layer 91 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6 Conclusion and future work 118 6.1 Hybrid RANS/LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.1 Control parameters . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.2 Resolved LES of the favorable pressure gradient boundary layer . . . 121 6.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Bibliography 124 iv List of Figures 1.1 Sketch of a hybrid method application with the boundary of the two domains parallel a) or orthogonal b) to the ow direction. . . . . . . . 5 1.2 Sketch of a hybrid method application. . . . . . . . . . . . . . . . . . 17 1.3 Sketch of the rescaling method (a) and recycling method (b). . . . . . 19 2.1 Distribution of the wavenumbers dni (top right) and modi ed wavenum- bers dnicn for di erent distances from the wall: top right gure is for y+ = 12, bottom left is y+ = 100, and bottom right is y = 5 . . . . . 48 2.2 Synthetic uctuation produced by the Batten?s decomposition for dif- ferent distance from the wall: bottom gure is for y+ = 12, central is y+ = 100, and top is y = 5 The distance are based on a LES calculation from which the Batten decomposition is evaluated. . . . . 49 3.1 Block diagram of the PI controller. . . . . . . . . . . . . . . . . . . . 51 3.2 Sketch of the geometric con guration. . . . . . . . . . . . . . . . . . . 54 3.3 (a) Domain-averaged turbulent kinetic energy for Tave = 1; KI = 5; KP = 30; Tave = 10; KI = 5; KP = 30; Tave = 100; KI = 1; KP = 1; Tave = 100; KI = 5; KP = 30. (b) Cf for the reference LES ( ) and the cases shown in (a). . . . . . . . . . . . . . . . . . . 56 3.4 Skin friction coe cient Cf with KI = 5 and KP = 30; rst control plane after one grid cell; rst control plane after one boundary layer thickness . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.5 Skin friction coe cient Cf with KI = 5 and KP = 30; hu0v0it is spanwise-averaged; hu0v0it is not spanwise-averaged. . . . . . . . 59 3.6 Bode plot of the frequency response of the PI system for di erent value of KI and xed KP = 30. KI = 5; KI = 20; KI = 30; It is evident that when KI increase, the gain at low frequency increase and the phase lag value starting at 90 is kept to higher frequency, i.e., the response of the PI system lags behind the input wave by 90 deg in the worst case . . . . . . . . . . . . . . . . . 62 v 3.7 Bode plot of the amplitude and phase frequency response of the PI system for di erent values of KP and xed KI = 5. KP = 5; KP = 20; KP = 30. When KP increases, the gain at high frequency increase and the phase lag value decreases at higher frequency 63 3.8 Domain-averaged turbulent kinetic energy. (a) KI = 5 and KP = 0; 15; 30; 500; (b) Kp = 30 and KI = 1, 5, 20, 30. . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.9 (a) Instantaneous and (b) integrated error for KP = 30 and KI = 5 (lines) and KI = 20 (lines with symbols). y+ = 13; y= o = 4:5. The curves for KI = 20 are shifted upwards by 0.005 in (a) and 0.1 in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.1 Zero-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; k RANS; Hybrid k /LES; SA RANS; Hybrid SA/LES. . . . . . . . . . . . . . . . . . . . . . . 70 4.2 Zero-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3 Skin friction coe cient Cf in the case of hybrid k /LES simulation. Reference LES; RANS k ; Hybrid k /LES with Tave = 10, KI = 5, KP = 30; Hybrid k /LES with Tave = 100, KI = 1, KP = 1 as in [30]; Hybrid k /LES with only Batten synthetic turbulence at the interface . . . . . . . . . . . . . . . . . . 72 4.4 Freestream velocity for the adverse-pressure-gradient calculation. . . . 74 4.5 Adverse-pressure-gradient boundary layer. (a) Skin friction coe - cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; k RANS, SA RANS; hybrid k /LES; Hybrid SA/LES. . . . . . . . . . . . . . . . . . . . . 76 4.6 Mean streamlines for the APG boundary layer. (a) Reference LES; (b) k RANS; (c) SA-RANS; (d) hybrid k /LES; (e) hybrid SA/LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 vi 4.7 Skin friction coe cient Cf in the case of LES/LES simulation (a) and hybrid RANS SA/LES (b). In (a) Tave = 38, KI = 5, KP = 30; Tave = 100, KI = 1, KP = 1 as in [30]. In (b) RANS SA model, hybrid SA/LES with Tave = 38, KI = 5, KP = 30; Tave = 100, KI = 1, KP = 1 as in [30]; Hybrid SA/LES with only Batten synthetic turbulence at the interface . . . . . . . . 79 4.8 Streamwise velocity uctuations in the APG boundary layer, in the y= o = 0:09 plane. (a) Reference LES; (b) hybrid k /LES; the red rectangle indicates the region where the control is active. . . . . . . . 80 4.9 Favorable-pressure-gradient boundary layer. (a) Skin friction coe - cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; hybrid SA/LES, Tave from LES data; hybrid SA/LES, Tave from RANS data; . . . . . . . . . 81 4.10 Skin friction parameter Cf for the favorable pressure gradient bound- ary layer: reference LES, RANS, hybrid RANS/LES with Tave = 2:4 and Tave = 1; Hybrid SA/LES with only Batten synthetic turbulence at the interface . . . . . . . . . . . . . . . . . . . 82 4.11 3D boundary layer. Flow angle = tan 1 U=W. Reference LES. SA-RANS; error based on hu0v0i; production-based error. 84 4.12 Grid resolution obtained in the 3D boundary layers +; + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.13 3D boundary layer. (a) Streamwise and spanwise skin friction coef- cients, Cf;x and Cf;z; (b) mean streamwise velocity pro les and (c) mean spanwise velocity pro les at the locations indicated by a ver- tical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. . . . . . . . . . . . . . . 86 4.14 3D boundary layer. (a) Streamwise and spanwise skin friction co- e cients, Cf;x and Cf;z; (b) secondary shear stress hv0w0i and (c) principal shear stress hu0v0i at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. . . . . . . . . . . . . . . . . . . 87 4.15 Contours of streamwise velocity fuctuations in a xz plane at y= o = 0:18. (a) Reference LES (only part of the domain is shown); (b) hybrid RANS/LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 vii 5.1 Sketch of the con guration. The computational domain is shown as a hatched area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Streamwise development of statistical quantities. (a) Freestream ve- locity, U1=Uo and friction velocity u =u ;o; (b) Acceleration param- eter K; (c) Momentum-thickness Reynolds number Re ; (d) Skin friction coe cient Cf. Experiments [97]; high-K case; low-K case; ZPG boundary-layer correlation. . . . . . . . . . . . 102 5.3 Wall-normal pro les of the mean streamwise velocity at the locations shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 0:5 units for clarity . . . . . . . 103 5.4 Wall-normal pro les of the mean streamwise velocity in inner coordi- nates at the locations shown in the top gure. Experiments [97]; high-K case; low-K case; U+ = 2:5 log y+ + 5.Each pro le is shifted by 18:0 units for clarity . . . . . . . . . . . . . . . . 104 5.5 Wall-normal pro les of the streamwise Reynolds stresses at the loca- tions shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 12:0 units for clarity . . . 106 5.6 Wall-normal pro les of the Reynolds shear stresses at the locations shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 0:8 units for clarity . . . . . . . 107 5.7 Contours of averaged quantities for the high-K case. The red line shows the local boundary layer thickness 99. (a) hu0u0i=U21; (b) hv0v0i=U21; (c) hu0v0i=U21; (d) Cuv = hu0v0i=urmsvrms. . . . . . . . . . . 108 5.8 Development of the Reynolds stresses along selected streamlines. (a) Streamline coordinates; (b) hu0u0i=U2o ; (c) hv0v0i=U2o ; (d) hu0v0i=U2o . The thick black line corresponds to the boundary-layer edge. Stream- line originating in: outer-layer; middle boundary layer; log region; viscous sublayer. . . . . . . . . . . . . . . . . . . . . 109 5.9 Development of a1 parameter and correlation coe cient along se- lected streamlines. (a) Streamline coordinates; (b) a1; (c) Cuv. The thick black line corresponds to the boundary-layer edge. Streamline originating in: outer-layer; middle of the boundary layer; logarithmic region; viscous sublayer. . . . . . . . . . . . . . 111 5.10 Instantaneous contours of u0 velocity uctuations in a plane parallel to the wall. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 viii 5.11 Instantaneous isosurfaces of Q = 0:002 in the outer layer colored by the streamwise vorticity. . . . . . . . . . . . . . . . . . . . . . . . 114 5.12 Contours of instantaneous u0v0 correlation in planes normal to the mean ow, and secondary (v w) velocity vectors. The two solid lines represent contours of U=U1 = 0:95 and 0.99. . . . . . . . . . . . 116 ix List of Abbreviations 3D Three-dimensional APG Adverse pressure gradient CPU Central processing unit DES Detached eddy simulation DNS Direct numerical simulation FPG Favorable pressure gradient LES Large eddy simulation LNS Limited Numerical Scales NLDE Non-linear disturbance equations NS Navier-Stokes RANS Reynolds averaged Navier-Stokes TKE Turbulent kinetic energy URANS Unsteady RANS WLES Wall-modeled LES ZPG Zero pressure gradient KP Proportional parameter Re Reynolds number based on the reference length and velocity scale Sij Strain rate tensor U A velocity scale Uo Free stream input velocity f generic function f ltered function k Turbulent kinetic energy p Pressure t Nondimensional time u (u1) Streamwise velocity component v (u2) Wall-normal velocity component w (u3) Spanwise velocity component x (x1) Streamwise coordinate y (x2) Wall-normal coordinate z (x3) Spanwise coordinate lter width Komlogorov eddy dissipation rate Kinematic viscosity Density ij The subgrid-scale stress w The wall stress Boundary-layer momentum thickness x Chapter 1 Introduction 1.1 Motivation Recent years have seen the spread of inexpensive parallel computer clusters that are able to reduce the computational time for the solution of numerical models; at the same time smart computational strategies that allow for faster simulations without compromising the accuracy of the results have been developed. These devel- opments have allowed the extension of the numerical solution of the Navier-Stokes equations to turbulent ows in more realistic con gurations than possible until now. The solution of turbulent ow problems is theoretically possible through the direct numerical simulation (DNS) of the Navier Stokes equation (with appropriate boundary conditions). This method constitutes the conceptually simplest approach to the problem of turbulence. Practically, however, the cost of DNS con nes this approach to simple application in terms of Reynolds number and geometry complex- ities. In fact, in DNS all the scales of motion must be resolved, from the integral (L) to the Kolmogorov ( ) scales. The grid size must be on the order of the smallest scale, and the computational domain must be on the order of the largest scale. This results in a number of grid points in each direction that is proportional to the ratio 1 between the largest and the smallest scale ( L ). De ning a Reynolds number based on the integral scale L, ReL, the number of grid points will be proportional to Re34L in each direction, so a total number of points proportional to Re94L is required by DNS. Moreover, the equations must be integrated for a time on the order of the integral time scale (T). Supposing that the time step of the calculation t is limited by the CFL condition, that is t < xV where V is the local velocity determined by the integral time and length scale, we have that the total number of points in time is T t again proportional to ( L ). In this way the total cost of a DNS calculation will be on the order of Re3L. This means that for high Re number, the DNS takes an unrealistic time to be completed: assuming that a computer power will increase by a factor of 5 every ve years, Spalart [1] estimated that DNS will not be applicable to the study of the ow over an airliner or a car until 2080. In large-eddy simulations (LES) the target is to simulate only the eddies con- taining the bulk of the energy of the ow. These scales are typically anisotropic and are dependent on the boundary conditions of the problem, therefore they are extremely di cult to be modeled analytically; the fact that the small, dissipative eddies are modeled reduces the cost of LES considerably compared with DNS. How- ever, when LES have to be applied to wall-bounded ows at high Reynolds numbers it is still very demanding because of the inner layer resolution. In fact, as found in Chapman [2], considering a at-plate boundary layer, the outer part has energy- carry structures on the order of the boundary layer thickness . Assuming that 2 the grid spacing is xed in the stream and spanwise direction and on the order of 0:1 , Chapman estimated a total number of grid points proportional to Re0:4. In the inner part of the boundary layer, on the other hand, the number of points must be based on the inner units; in fact, the dimension of quasi-streamwise vortices are constant in wall units (i.e., normalized by kinematic viscosity, wall stress and uid density). In this way, the total cost of the calculation in the inner layer is estimated be proportional to CfRe2L which, assuming Cf Re 0:2L , makes again this technique only suitable for moderate Reynolds number applications. A powerful way to use LES is either to apply it in non-bounded ow simula- tions, or to use LES calculations only on the outer part of the ow. Wall-modeled LES (WMLES), in which the inner layer is modeled either by the Reynolds-Averaged Navier-Stokes (RANS) equations or through approximate boundary conditions, can be a smart solution to the expensive full LES (see the review in Ref. [3]). Spalart [1], however, estimates that application of WMLES to external ow of aeronautical in- terest is still several decades away from being practical. His estimate assumes the presence, in these ows, of very thin boundary layers (near the wing leading-edge, for instance), where a large number of points is required to cover a small physical area. In environmental or oceanographic ows, in which the Reynolds numbers are comparably high, but the boundary layers are thicker, WMLES is already applicable to practical problems. Traditionally, high-Reynolds number ows have been predicted using RANS 3 models. However, RANS is designed to be accurate for problems such as thin shear layers. In many cases, especially in the presence of uid-dynamical non-equilibrium, RANS cannot give an accurate prediction. A possible avenue to solve uid- ow problems in practical engineering cases, then, is by using a combination of RANS and LES, in which the RANS equations are solved in quasi-equilibrium regions or where the models can be accurately tuned, while LES or WMLES are used in non-equilibrium regions or where the RANS models are expected to fail. Hybrid LES/RANS techniques have emerged as a tool that allows the simulation of complex uid ows within reasonable amounts of CPU time. In fact, in recent years the hybrid methods have received increasing attention; di erent terminologies have been used as LNS, DES, PANS, zonal approach, single grid approach, etc. The target of this chapter is to introduce and discuss the di erent methods and try to underline the drawbacks and key points of each. A review of previous work is, in some cases, necessary in order to understand the motivation of the new hybrid strategies. In the next section, an overview of hybrid strategies is presented. Detailed descriptions of each method can be found in Sections 1.3 and 1.4. The chapter is concluded with an analysis of the hybrid method used in the rest of the present dissertation. 1.2 Hybrid strategies: single grid and zonal approaches It is conceptually possible to divide the hybrid RANS/LES methods into two main categories: one in which a single grid is used and only the turbulence model 4 Figure 1.1: Sketch of a hybrid method application with the boundary of the two domains parallel a) or orthogonal b) to the ow direction. changes; another case occurs when the RANS/LES domains are separated, and two separate grids as well as sets of equations are used (zonal approach). In both cases, two di erent scenarios are possible, as shown in gure 1.1: the interface of the two domains can be parallel or orthogonal to the ow direction. This interface can be either a true discontinuous interface that separates the two domains as in the zonal methods, or can be representative of a region in which the model gradually switches from RANS to LES calculations as in the single grid methods. One issue that arises when hybrid RANS/LES methods are used is the behavior 5 of the ow in this transition zone between the RANS and LES regions, and how that region could a ect the downstream results. In the RANS zone the ow solution is either steady, or only contains information on the largest scales of motion if unsteadiness is present; most, or all, of the Reynolds shear stress is provided by the turbulence model. In the LES region, on the other hand, the resolved scales must supply most of the Reynolds shear stress and small scale structures must be present to provide it. In both the zonal approach where the interface between the two domains is explicit, and the single-grid method where the interface is an overlapping zone between RANS/LES, a transition zone exists in which the resolved, energy-containing eddies are gradually generated and grow. In this region the ow may be unphysical. Di erent hybrid methods, that are illustrated in this paragraph, have di erent ways to switch from RANS and LES to trigger instabilities and to develop small eddies to support LES calculations. In the next section each method will be addressed with emphasis on the literature available. In the single grid approach, the switch from RANS to LES is typically smooth and determined by the damping of the turbulent viscosity: in the LES region it must be damped to the levels implied by the sub-grid model. The damping can be a function of the grid spacing, as in DES where the level of the e ective sub-grid viscosity is a function of the local mesh spacing and/or distance from the wall; or the damping is a function of turbulence quantities as in the LNS of Batten et al. [8] where the turbulent viscosity is damped through a blending function that depends 6 on length and velocity scales; the damping can also be a combination of the two (Speziale?s latency parameter [7]). In the zonal approach, on the other hand, the switch between RANS and LES is, by de nition, explicit: two separate domains are de ned and two di erent sets of equations, RANS and LES, are solved. The problem, in this case, is how to transfer information between the two domains: appropriate boundary conditions are needed for each domain; these can give either a two-way coupling (RANS solution in uences the full LES and vice versa ) or a one-way coupling. As mentioned above, another issue is the generation of the energy-carry eddies in the transition zone. If a natural instability is present in the RANS/LES transition area the small eddies are generated by the natural ampli cation of noise or small uctuations. This instability is particularly e ective if adverse pressure gradients, separated shear layers or other very unstable ow features are present. In the case where such features are not present, other techniques must be implemented to generate turbulent uctuations quickly. In the next sections we will describe some of these methods: forcing methods that add energy in the LES domain through the addition of forcing terms in the momentum equations, the method that recycles and rescales uctuations from the downstream part of the calculation, and techniques that introduce coherent structures taken from a precursor calculation. Each of these methods will be illustrated in the next sections either applied to a single grid calculation or to zonal hybrid methods. 7 1.3 Hybrid RANS/LES approach with single grid 1.3.1 Introduction As pointed out in the rst section, the bottleneck of LES is the calculation of the ow dynamics close to the wall. To avoid this costly requirement, it is necessary to model the inner part of the ow in a less expensive way, from a computational perspective, than the full LES solution. Many hybrid methods try to achieve this end, by using RANS in the boundary layer (or in the inner layer) and LES away from solid walls. In the following, an overview of the literature of single grid approaches is presented. 1.3.2 Previous work One of the rst theoretical approaches to hybrid RANS/LES was due to Speziale [7] where the model Reynolds stress was damped by a latency factor that was a function of the grid size and the Kolmogorov scale. In this way, Speziale was able to merge RANS and LES depending on the grid size and/or Reynolds number. However, some parameters in the latency factor were never completely speci ed by Speziale. Moreover, the choice of the Kolmogorov length scale and the grid size as parameters to damp the eddy viscosity appear to be too drastic: for high Re numbers we can expect that no grid is ne enough to result in an LES calculation ( since the ratio between grid size and Kolmogorov scale remains large). The latency 8 factor of Speziale [7], appears to be more of an on/o switch from RANS to full DNS because of the Kolmogorov length scale used in the model. The method proposed by Batten et al., the Limited Numerical Scales (LNS) approach [8], is an attempt to improve Speziale?s approach: they de ned a parameter- free latency factor, based on some characteristic length and velocity scale of LES and RANS models. In the LNS approach the unresolved stress is a fraction of that predicted by a RANS model; the fraction was determined by a ratio of LES and RANS length- and velocity-scales. In this model the blending of the eddy viscosity is based entirely on RANS and LES quantities with no other empirical constants in- volved. One of the issues of the LNS approach regards the generation of small eddies in the LES region, when no natural instability, such as those present in massively separated ow, was present in the domain. In order to accelerate the generation of turbulent structures, Batten et al. [31] introduced a method to generate syn- thetic turbulence, based on Kraichnan?s [46] proposal, that takes into account the anisotropy of the ow. The method by Batten et al. [31] is based on the superpo- sition of sinusoidal modes with random frequencies and wave-numbers, with given moments and spectra. The wave-numbers are modulated to yield eddies that are more elongated in the direction of larger Reynolds stresses, thus introducing more realistic, anisotropic eddies into the ow. Batten et al. [31] performed calculations of the spatially developing ow in a channel in which the inlet section had a very coarse grid in the streamwise direction, capable of supporting a RANS solution but 9 not an LES. The grid was then suddenly clustered, and forcing terms were used to generate synthetic turbulence. The synthetic turbulence generation method proposed by Batten et al. [43, 31] was signi cantly more e ective than older, simpler, methods. For instance, Le et al. [44] performed calculations of a backward-facing step in which at the in ow they assigned a mean velocity pro le plus a superposition of random uctuations with given moments and spectra. The amplitude of the random uctuations was such that the bulk of the energy was contained in a range of well-resolved wave-lengths [45]. Since the uctuations lacked phase information, however, the turbulence levels de- cayed rapidly, and only some distance downstream the turbulent eddies regenerated. The most popular of the hybrid RANS/LES techniques that fall in the \single- grid category" is Detached Eddy Simulation [5]. In this model the switch between RANS and LES is determined by the grid size: when the mesh becomes small enough to resolve the energy-carrying eddies the eddy viscosity is reduced. DES was designed for the simulation of massively separated ows, in which the integral scales of motion in the separated ow region are determined by the geometry, and their size is not signi cantly a ected by the Reynolds number. Applying LES in this region, therefore, results in a relatively small increase of the cost of the calcu- lation compared, for instance, with an Unsteady RANS (URANS) model. Thus, in standard applications of this method, the thin attached shear layers are modeled by RANS, and only away from solid bodies the technique switches to LES, and the 10 RANS/LES interface occurs generally in the separated shear layers. The strong instability caused by the in ection points in the velocity pro les is extremely ben- e cial from the point of view of the generation of energy-carrying eddies, since any disturbance present in the ow is quickly ampli ed, resulting in a very short transi- tion region. The application of DES to ows with weaker instability mechanisms or with thick boundary layer, however, creates some problems: if the turbulent eddies do not grow su ciently fast, the simulations can become inaccurate as the mod- eled shear stress decreases while the resolved one does not increase su ciently fast. Such quantities as the mean velocity pro le, skin friction coe cient and separation location can be signi cantly in error. This case is typical [9] when the grid cells parallel to the wall are re ned and become of the same order of the boundary layer thickness within the boundary layer region, with the boundary layer attached to the wall. In this case, within the boundary layer, the grid spacing is ne enough to switch the calculation to LES solution by decreasing the amplitude of the eddy viscosity, but the turbulent eddies are not yet generated. Another challenging case is when the grid is re ned enough to switch the model to LES, but not ne enough to support resolved velocity uctuations internal to the boundary layer: in that region the solution is neither a pure RANS, nor a pure LES. In these kinds of scenarios, the grid plays a key role: it must be ne enough to support eddies and damp in an appropriate way the eddy viscosity; this double achievement is the cause of error of the DES model (gray region). 11 At this point it is important to stress that the problem mentioned in the DES calculation is not solely speci c to that kind of method. In fact, as mentioned in the previous section, a transition zone exists in all of the hybrid RANS/LES methods, in which energy-containing eddies are gradually generated and grow. This issue is of outstanding importance for the prediction of the correct behavior of the solution, especially in the case where separation of the ow is not determined by the geometry of the problem. One rst solution to the problem mentioned was already illustrated: the Batten et al. [31] decomposition. In the following paragraphs, we will present other solutions to overcome the limitations of the gray area. The problems of the gray area have motivated the use of zonal-DES approach [11], in which attached boundary layer regions are modeled by RANS independently of the grid spacing. In this way, the user can re ne the grid as much as desired in the LES region, and use RANS in all of the regions where an attached boundary layer is expected. This \manual" switching is, however, too elaborate to be applied in a complex 3D geometry and could be a source of error. Recently, to overcome the problem of the gray region, Spalart et al. [9] mod- i ed the DES approach and introduced the Delayed Detached Eddy Simulation (DDES), which maintains the RANS solution in the boundary layers independently of the grid spacing. They started from the idea of Menter et al. [75], who used appropriate functions in the SST RANS model to identify the boundary layers and prevent the switch to the LES model within that region; the DDES of Spalart et al. 12 [9] is applicable to any RANS model that involves an eddy viscosity. Spalart et al. [9] performed calculations on boundary layers, on a single and multi-element airfoil, a cylinder, and a backward-facing step and proved that the RANS solution was maintained in thick boundary layers, and switched in LES after massive separation. Another method that used one grid is the Partially Averaged Navier Stokes (PANS) solution by Girimaji et al. [10]. In this method the averaging of the Navier Stokes equation is performed only over a portion of the uctuating scales. The scale-resolution (cuto ) of the PANS closure is controlled by two model parameters: the fraction of unresolved kinetic energy and unresolved dissipation. The model production-to-dissipation ratio must be consistent with the turbulence physics at the cuto to give accurate results. Girimaji et al. [10] applied this method in a ow past a circular cylinder and a ow past a backward facing step with encouraging results. The trend in the most recent hybrid solutions to the one-grid hybrid RANS/LES is to apply a technique that resembles the wall-modeled LES. As described in [9], wall-modeled LES (WMLES) was introduced in the 1970s. An extensive review of wall-layers models can be found in Cabot et al. [13] and Piomelli et al. [3]. Here only a brief summary will be given. An early application of WMLES can be found in Balaras [6]: the Two-layer model (TLM); here a simpli ed (parabolized) version of the momentum equations was solved in the inner layer together with a turbulence-model equation; the bound- 13 ary conditions were taken from the external LES. The solution of the simple set of equations provided the wall-shear stress imposed in the LES calculation. In this case the coupling between the near-wall region and the outer region was supposed to be weak from inner to outer layer and strong in the opposite direction (two way coupling). Nikitin et al. [12], used DES as a wall-layer model in high Re numbers of turbulent channel ow, letting the model switch from RANS to LES behavior inside the boundary layer. In the inner layer the Reynolds shear stress was provided entirely by the turbulence model. In the outer part, the eddy viscosity was damped, which allowed the formation of turbulent eddies. The results that they obtained were promising: a clear logarithmic law was observed, and turbulence in the outer layer was sustained with a grid not particularly ne. The skin friction coe cient, however, was under-predicted. Improved results in an attached ow were obtained when backscatter forc- ing was used [15] to introduce small-scale uctuations in the transition region. This stochastic forcing generated small-scale uctuations that acted as \seeds" for the de- velopment of realistic, energy-carrying eddies in the LES region. Piomelli et al. [15] found that with the correct amount of backscatter they could successfully remove the shift in the log-law, and improve the prediction of the skin friction coe cient. Although this work showed a promising approach to the solution of the RANS/LES interface problem, the backscatter model proposed lacked physical justi cation, and 14 tuning the forcing magnitude was required as the grid or the Reynolds number were changed. To overcome this problem, Keating et al. [20] introduced a dynamic method to calculate the magnitude of the stochastic forcing, based on the fact that in the interface region, the resolved Reynolds shear stress should become larger than the modeled one. They applied the method to a turbulent channel ow. Davidson et al. [18] and Dahlstrom et al. [23], add turbulent uctuations, synthesized or from a DNS precursor simulation, to the momentum equation at the interface region RANS/LES to generate rapidly turbulent structure in the LES region. The source terms are scaled in order to match turbulent kinetic energy mod- eled with the RANS. They found that the turbulent length scale of the synthesized uctuation has a large impact on the results. In these cases, the turbulent length scale and the eddy viscosity are not the same in the LES and RANS domain. Tessicini et al. [21] instead, used an URANS close to the wall and an LES subgrid scale model in the outer part where they adjusted dynamically the constant C in the turbulent viscosity model of RANS to have a continuity between the two regions. The smooth transition between the constant value of C in the inner RANS domain to the value calculated dynamically, is obtained by an empirical smoothing exponential function. As can be seen from the last authors cited, better results were in general obtained by adding energy at the interface where the switch from RANS to LES was supposed to be. This is a problem very similar to the other approach zonal 15 RANS/LES, where an explicit interface between the two domains exist. This zonal approach will be presented in the next section. 1.4 Zonal Hybrid RANS/LES 1.4.1 Introduction We focus, in this work, on applications of hybrid RANS/LES methods of the type shown in Figure 1.2, in which RANS is used in equilibrium regions of the ow, while LES is performed in a small part of the domain, possibly including attached boundary layers as well as separated ow. As shown in Figure 1.2, the separation between RANS and LES calculations is pre-de ned, so a smooth transition of the eddy viscosity between the two domains does not have to be de ned. The interface RANS/LES de nes the two regions in a discontinuous way: the turbulent kinetic energy is zero in the RANS domain and not zero in the other. The main di culty is now that information must be exchanged between two solutions radically di erent in their energy content: RANS or unsteady RANS (URANS) in one domain and LES in the other. The interface problem is critical for applications in which the upstream ow is essential for predicting the correct behavior of the ow downstream; one case is when the ow undergoes separation, and errors in the separation prediction can a ect the ow for long distances downstream. In this type of con guration the generation of energy- and momentum-carrying eddies is expected to be much slower, 16 RANS LES Figure 1.2: Sketch of a hybrid method application. as the instability mechanism in the attached boundary layer is weaker than in the separated shear layer. Using the RANS solution alone to feed the LES domain is not su cient, but further modeling must be included in order to have an accurate LES simulation. One way to couple RANS and LES is by adding an external source of energy on the boundary of the two domains for the generation of small scales needed in the LES domain. The external source of energy can be either a stochastic reconstruction of the turbulent uctuations, or a precursor simulation to feed the required unsteady ow to the LES domain. These reconstructions, of course, must match the upstream eld given by the RANS model in a statistical sense. The problem of generating turbulent uctuations for the LES domain is anal- ogous to the problem of in ow conditions for LES and DNS simulations. In the following subsection we will describe the most common methods of in ow generation for spatially developing problems that have been used in LES. We will particularly 17 stress those methods that have a more immediate application to hybrid RANS/LES calculations. 1.4.2 In ow generation methods The rst applications of LES were to temporally developing ows, in which it was possible to use periodic boundary conditions in the direction of homogeneity. To extend LES to the more realistic spatially developing cases, initially modi cations of the periodic boundary conditions were developed that either added extra terms to the governing equations to account for the boundary-layer growth [34, 35, 36] or rescaled the uctuations before applying periodicity [102, 38] (see Figure 1.3 a). These methods have some shortcomings: rst, they require a longer computational domain, and can be applied only if some similarity law exists, at least in the initial part of the domain. Furthermore, it is not easy to control precisely the in ow parameters (boundary-layer thickness and wall stress, in particular) Schl uter et al. [40] adopted a method that improved the recycling by Lund et al. [102]: they introduced uctuations that already have some level of correlation at the in ow of the LES. To this end they speci ed the in ow by superposing to a mean ow obtained from the RANS, the uctuations obtained from a separate LES of a wall-bounded ow. These uctuations were rescaled to match the Reynolds stresses predicted by the RANS. A similar method, that uses an auxiliary calculation in which the mean pro le is assigned to match the RANS result was recently proposed 18 Figure 1.3: Sketch of the rescaling method (a) and recycling method (b). 19 by Medic et al. [41]. Another method to generate in ow conditions, (see Figure 1.3 b) consists of running a separate, precursor, calculation of an equilibrium ow in which either periodic boundary conditions or recycling arguments can be used. Then, a plane of velocity data orthogonal to the ow direction is stored at each time step. The sequence of planes is then read-in as in ow data for a separate calculation of the ow of interest, with a necessary space and time interpolation. This method has the advantage, over recycling methods, that the control of integral parameters can be achieved precisely but may require considerable additional computational resources. Recycling and precursor techniques are very e cient for simple geometries; however, the application to complex three-dimensional geometries, non-equilibrium ow and arbitrary meshes presents some challenges [43]. Because of the limitations of these approaches, several researchers have attempted to develop methods to generate realistic in ow conditions from synthetic turbulence. Le et al. [44] performed calculations of a backward-facing step in which at the in ow they assigned a mean velocity pro le plus a superposition of random uctuations with given moments and spectra. The amplitude of the random uctu- ations was such that the bulk of the energy was contained in a range of well-resolved wave-lengths [45]. Since the uctuations lacked phase information, however, the tur- bulence levels decayed rapidly, and only some distance downstream the turbulent eddies regenerated. 20 Synthetic turbulence generation methods like that one by Batten et al. [31] described in the previous section, have also been proposed by Smirnov et al. [47] and Klein et al. [48]. These methods achieve improved results (compared with the simple superposition of random uctuations) by a judicious assignment of the turbulence spectrum and by trying to match the ow anisotropy. An important feature of turbulence in wall-bounded ows is its structure, however; since none of these method contains realistic phase information between the modes, a fairly long adjustment region downstream of the in ow is unavoidable. In this region the initially random uctuations are selectively ampli ed by the ow, and realistic turbulent eddies are generated. A procedure to obtain ow eld with temporal correlation was applied by Davidson [32]; he generates new uctuating velocity with a prescribed time expo- nential correlation with a time scale proportional to the turbulent time scale. He applied the methods to an hybrid LES/RANS channel ow and found that the time scale was more important than the length scale and that both should not be based on physical values but related to the grid. Another method to generate synthetic turbulence was adopted by Sandham et. al [33]. It is based on the fact that outer and inner part of a boundary layer are dominated by large scale coherent structures and low speed streaks respectively; they, therefore, introduced streamwise and spanwise disturbances with deterministic phase information, and frequencies and wave numbers chosen to match typical size of 21 the coherent structures in the inner and outer part of the boundary layer. Moreover, the spanwise component was computed using the divergence-free condition. To break any remaining symmetries in the in ow condition, random noise was added. They tested his method for a at plate turbulent boundary layer. Recently, Mathey et. al [28] apply a vortex method to generate random uctu- ations representing a turbulent eld at the inlet of an LES domain. The generated velocity eld has temporal and spatial correlation. A perturbation is added on a speci ed mean velocity pro le via a uctuating two dimensional vorticity eld. Mathey et. al [28] validate the vortex method on fully developed turbulent chan- nel ow, pipe ow and separated hill obtaining results that compare well with the reference data. However, they only used methods based on simple random noise for comparision. These methods achieve improved results (compared with the simple superpo- sition of random uctuations) by a judicious assignment of the turbulence spectrum and by trying to match the ow anisotropy. As a result, their downstream develop- ment is more realistic than when random noise was used. In all cases, however, we continue to observe an initial decay of the turbulence levels (albeit shorter than with simpler speci cations) followed by an eventual establishment of realistic turbulence. This was shown by Keating et al. [50], who compared various types of in ow conditions for LES, among them the use of random noise, the adapted database method by Schl uter [40], the synthetic turbulence generation method proposed by 22 Batten et al. [31], and a method based on controlled forcing proposed by Spille- Koho and Kaltenbach [49]. This method is based on the selective ampli cation of strong bursts downstream of an in ow supplied by some synthetic turbulence generation method. A controller is used at several location downstream of the in- ow to determine the amplitude of a forcing term in the wall-normal momentum equation. This forcing term acts to reinforce the more realistic eddies, by requiring that a desired Reynolds shear-stress pro le (obtained from the RANS, or from ex- periments) be achieved. They found that the adapted database and the controlled forcing methods resulted in signi cantly shorter development lengths than any of the other techniques. 1.4.3 Previous work Applications of the type sketched in Figure 1.2 have been studied by zonal hybrid RANS/LES methods by several researchers. Labourasse and Sagaut [24], for instance, developed a coupling method to perform LES overlaid on a mean ow generated by RANS. The RANS parametrization yields the mean ow eld over the entire region, while the LES equations, which are formulated in perturbation form (Non-Linear Disturbance Equations | NLDE) are solved only in a small region of the ow. They applied this method to stationary and pulsed channel, as well as to the ow over a turbine blade (a con guration conceptually similar to that shown in Figure 1.2). The turbulent uctuations in their calculation were generated 23 by the natural ampli cation of acoustic perturbations and numerical errors in the adverse-pressure-gradient region on the suction side of the blade. The authors did not quantify how e ective this ampli cation was (i.e., what distance was required to develop a realistic turbulent ow. Terracol [25] performed calculations of the ow over a at-plate trailing edge with thickness comparable to the boundary layer thickness and a NACA0012 airfoil also using the NLDE method [24]. In the rst con guration he compared various treatments of the RANS/LES interface: allowing the ampli cation of existing in- stabilities to develop (the method used in [24]), a recycling method [26, 102], and a synthetic reconstruction of the uctuation eld [27]. He found that, in this ow, the ampli cation of existing disturbances was not su cient to establish a well-developed turbulent ow prior to the trailing edge, while the recycling resulted in correct turbu- lent statistics, but produced spurious peaks in the pressure spectrum; the synthetic reconstruction o ered the best result, and gave reasonable pressure spectra in the NACA0012 case. Schl uter et al. [39] presented results of the hybrid simulation of the ow in a gas turbine engine in which both the turbine and compressor were simulated by RANS while the combustor was computed using LES. The calculation was performed in separate stages: the out ow from the RANS solution of the compressor supplied in ow conditions for the LES, while the out ow of the LES was used to supply inlet conditions to the RANS calculation of the turbine. In order to generate turbulence 24 at the RANS/LES interface the authors used the results of a related investigation of in ow conditions for LES by Schl uter et al. [40]. Recently, Keating et al. [30] applied the controlled forcing method by Spille- Koho and Kaltenbach [49] coupled with the synthetic turbulent by Batten et al. [31] in a hybrid RANS/LES framework for the simulation of boundary layers in favor- able and adverse pressure gradients, and found that it produced physically realistic turbulence in short distances. They also observe that the quality of the RANS data used a ects the results signi cantly; this is true in particular in favorable and adverse-pressure-gradient boundary layers. 1.5 Plan of the present dissertation The controlled forcing methods has shown very good performance to generate realistic eddies in a short distance. This method, moreover, can be used both in a zonal calculation, as an in ow condition and in a sigle-grid method, through forcing terms added to the momentum equations. The purpose of this dissertation is to improve the controlled forcing method[49, 50, 30] by optimizing the controller parameters (which previously were assigned by trial-and-error) to shorten the transition region and the simulation transient, and also to explain their physical signi cance. The test cases used to evaluate the model will be a at-plate boundary layer subjected to zero, adverse and favorable pressure gradients, and a three-dimensional boundary layer. An additional contribution of 25 this work was the study of the favorable pressure gradient boundary layer. This ow was initially computed as a test case for the hybrid calculations. The richness of the physics we encountered prompted us to look at this problem in additional depth. In the following chapters we will rst describe the methodology used, then discuss the controller function, and the optimization of its parameters. The result of four applications will be presented next, followed by conclusions and recommen- dations for future work. The following chapter will describe some results relating to the relaminarization and retransition of boundary layers subjected to strong accel- eration. Finally, conclusions and recommendations for future work will be made. 26 Chapter 2 Problem formulation 2.1 Governing equations The numerical methods used in this work is the same as [51] The reference system is cartesian with x or x1 being the streamwise direction, y or x2 the wall- normal direction, and z or x3 the spanwise direction, and with velocity components u, v and w or u1, u2 and u3 respectively. The equations are in non-dimensional form, with the following de nition: x i = xiL u i = uiU o t = tUoL p = p LU2 o Re = LUo where L and Uo are the reference length and velocity. The non-dimensional equations can, therefore, be written as @ui @xi = 0 (2.1) @ui @t + @ujui @xj = 1 Re @2ui @xj@xj @p @xi (2.2) where the is omitted. The large-eddy simulation equations are derived from the previuous equation by appling a ltering operation, de ned as f (x) = Z f (x0) G x; x0; dx0 (2.3) 27 where G is a lter function with a characteristic length, . The ltering operation separates the turbulent motion into large and small scales. The lter functions used include the sharp Fourier cuto lter, the Gaussian lter and the top-hat (or box) lter. Both the Fourier cuto and Gaussian lter are de ned in spectral space, while the top-hat lter is de ned in physical space. As nite di erence methods are used in the present work, it is computationally more e cient to use the top-hat lter, which is de ned in physical space as G(x) = 8 >>> < >>> : 1= if jxj =2 0 otherwise (2.4) Appling the lter operator to the Navier-Stokes equation gives: @ui @xi = 0 (2.5) @ui @t + @ujui @xj = 1 Re @2ui @xj@xj @p @xi @ ij @xj (2.6) where the subgrid-scale stress term, ij is de ned as ij = uiuj uiuj (2.7) This term is unclosed and must be modelled. The closure of the subgrid-scale stress terms is known as subgrid-scale modelling. In this work we use an eddy- viscosity closure with dynamic determination of the coe cient. After describing the numerical technique used to solve the governing equations (2.6) and the SGS model employed, the nal subsection details the averaging procedures required to remove the numerical instability arising from sharp uctuations in the coe cients. 28 2.2 Numerical method and boundary conditions The Navier-Stokes equations are integrated in time using a fractional step tech- nique [62]. The wall-normal viscous and subgrid-scale di usion terms are advanced implicitly using the Crank-Nicolson method. The other terms in the equations are advanced using a low-storage third-order Runge-Kutta scheme [61]. The three-step (k = 1; 2; 3) fractional steps, time advancement scheme is, therefore, bui k tM (bui) = rk 1i (2.8) @2 ( p) @xi@xi = 1 2 k t @bui @xi (2.9) uki = bui 2 k t@ ( p)@x i (2.10) pk = pk 1 + p (2.11) where the explicit terms, rk 1i , are rk 1i = uk 1i + t kM uk 1i + kA uk 1i + kA uk 2i 2 k @p k 1 @xi (2.12) At the rst substep, when k = 1, uk 1 = un, and at the end of the third substep, the solution is un+1 = uk. The coe cients of the time advancement scheme are 1 = 4=15 2 = 1=15 3 = 1=6 (2.13) 1 = 8=15 2 = 5=12 3 = 3=4 (2.14) 1 = 0 2 = 17=60 3 = 5=12 (2.15) This time advancement scheme is stable under the following conditions: = t u i xi < p3 (2.16) 29 where i = 1; 2; 3 and = 4 Re t x2i < p3 (2.17) where i = 1; 3. A second-order nite di erence scheme is used for spatial discretisation of all terms [63]. While spectral methods are more accurate and less di usive than nite di erence schemes, they are more di cult to use in complex geometries. Moin and Mahesh [64] found that equivalent results to a spectral method can be obtained using a second-order nite di erence method if the mesh is doubled in each direction. Higher-order, conservative nite di erence schemes have been proposed [63]; however these have additional complications in the use of nonperiodic boundary conditions. A staggered grid [65] is used for its good conservation properties [63] , its accuracy and to prevent the \odd-even" decoupling of pressure. The following discrete di erencing and averaging operators [63] are de ned on the staggered grid: 1 1x i i+1=2 i 1=2x i+1=2 xi 1=2 (2.18) 1x i i+1=2 + i 1=22 (2.19) These operators are then used to discretise the Navier-Stokes equations yielding 1ui 1xi = 0 (2.20) 1ui 1t + 1uj1xiui1xj 1xj = 1 Re 1 1xj 1uj 1xj 1p 1xi 1 ij 1xj (2.21) 30 The subgrid-scale coe cient is evaluated using the dynamic procedure at each timestep. The subgrid-scale eddy-viscosity, and the scale-similar contribution to the subgrid-scale stress is evaluated at every Runge-Kutta substep. All the results are averaged in time for a long time period (at least 15 Large eddy turn over times) in order to have statistically steady quantities. The samples of turbulent quantities for each cases are recorded after the transient time is over. Periodic boundary conditions are used in the spanwise direction, with no-slip boundaries at the solid wall, and convective conditions [82] used at the out ow. The freestream condition for the zero pressure gradient boundary layer and the 3D boundary layer [102], @u @y = 0; v = U1 d dx ; @w @y = 0; (2.22) was used at the upper boundary, where is the boundary layer displacement thick- ness and U1 is the freestream velocity. We calculated d =dx using a linear re- gression on the calculated (x) distribution. In the case of favorable or adverse pressure-gradient boundary layer the U1 was assigned and the following expression were used: @u @y = @v @x; v = U1 d dx + ( h) d U1 dx ; @w @y = 0; (2.23) with h the higher of the domain. In the adverse pressure gradirnt boundary layer case, V1 was assigned. Other boundary condition are described with the speci c cases in Chapter 4. When a derivative of a ow variable is required at a wall (i.e., 31 for the wall-normal di usion term), a three-point one-sided approximation is used [66] The Poisson equation (2.9) resulting from the fractional step technique is solved using Fourier transforms in the spanwise directions followed by cyclic reduc- tion, and a direct tridiagonal matrix inversion in the wall-normal direction. Neu- mann boundary conditions are used for no-slip walls. The code is parallelized using Message Passing Interface (MPI). The compu- tation is distribuited between n processors (n was varied between 1 and 4 for all the simulations described in this dissertation). 2.3 Subgrid-scale modeling The term ij denotes the subgrid-scale (SGS) stress, which represent the cou- pling of small scales on the resolved turbulence. The closure of the subgrid-scale stress is known as subgrid-scale modeling. In the following subsection the eddy- viscosity closure is presented, and in the next subsection the dynamic modeling is described. Dynamic modeling removes the need for a priori speci cation of the model coe cient. In the last subsection the averaging procedure required to damp sharp uctuations in the model parameters is presented. 32 2.3.1 Eddy-viscosity assumption Eddy-viscosity subgrid-scale models relate the subgrid-scale stresses, ij, to the large-scale (i.e., resolved) strain rate tensor, Sij = 12 @ui @xj + @uj @xi , through a subgrid-scale eddy-viscosity, sgst : ij ij3 kk = 2 sgst Sij (2.24) The simplest eddy-viscosity subgrid-scale model is Smagorinsky?s (1963), where the eddy viscosity is modeled as the product of a length scale (proportional to the lter width, ) and a velocity di erence at that scale (proportional to jSj, where jSj = 2SijSij 1=2). A constant, Cs, is included, giving the following model for the subgrid-scale eddy viscosity: sgst = C2s 2jSj (2.25) The Smagorinsky constant must be evaluated in order to apply the model. In isotropic turbulence, if the lter width is in the inertial subrange, the constant takes values between Cs = 0:179 and Cs = 0:23 [52]. However in the presence of shear, near solid boundaries or near transition this value needs to be decreased. Methods such as van Driest near-wall damping [53] and intermittency factors [54] have also been used. Eddy-viscosity models can represent the dissipative e ects of the small scales accurately. However, they fail to reproduce the local stresses with adequate accuracy. A priori studies comparing the actual subgrid-scale stresses calculated using DNS 33 and those calculated using an eddy-viscosity subgrid-scale model show a correlation coe cient of between 0 and 0.25 [55, 56]. 2.3.2 The dynamic procedure Dynamic models evaluate the subgrid-scale model coe cient from information contained in the resolved- ow eld of a large eddy simulation. In this way the dynamic of the smallest resolved scales are implicitly assumed similar to the dynamic of the subgrid scale. The advantages of the dynamic model are that is not necessary to tune any constant for di erent ow eld problem, that the subgrid scale stress vanish in a laminar ow and at a solid boundary, and that they have the correct asymptotic behavior in the near wall region of a turbulent boundary layer. Thus the introduction of the dynamic procedure by Germano et al. [100] resolved one of the major di culties in using subgrid-scale models { the need for a priori input of model coe cients. Using the dynamic procedure, the coe cients in the model are determined during the calculation and are based on the energy content of the smallest resolved scales. A test lter is de ned with a width, b , larger than the grid lter width, . The calculation of the coe cients in the model is based on the Germano identity [100] : Lij = Tij b ij (2.26) 34 which relates the resolved turbulent stresses, Lij, Lij = duiuj buibuj (2.27) to the subgrid-scale stresses, ij, (equation (2.5)) and the subtest stresses, Tij, Tij = duiuj buibuj (2.28) If an eddy-viscosity model is used to parameterize the subgrid-scale and subtest stresses, ij = 2Cev 2jSjSij (2.29) Tij = 2Cev b 2jbSjbSij (2.30) then Germano?s identity ( 2.26) yields the following equation for Cev: CevMij = Lij (2.31) where Mij = a2 2jbSjbSij + 2 djSjSij (2.32) and a = b = is the ratio of the test lter width to the grid lter width. As equation( 2.31) has a single coe cient and ve independent equations, the system is overdetermined. Lilly [60] proposed a least-squares minimization of the error, giving the following equation for determining the coe cient: Cev = 12 hLijMijihM ijMiji (2.33) where h i is an appropriate average. 35 2.3.3 Averaging of the eddy-viscosity coe cient The averaging in the dynamic procedure is required to remove the very sharp uctuations in the coe cient which can lead to signi cant errors, not to mention numerical instability. These sharp uctuations are caused by a mathematical incon- sistency in the expression for Cs in the dynamic procedure. Meneveau et al. [80] proposed averaging along pathlines in order to reduce the noise. To have Galilean invariance, the time average have to be performed by following uid particles of the ow. This average method is known as the Lagrangian dynamic model, allowing for averaging in inhomogeneous ows in complex geometries. Here, the average is de ned as hfi = If = Z t 1 f(t0)W(t t0) dt0 (2.34) where the integral is carried out following a ow pathline, and W(t) is a weighting function. These are then used to replace the averages in the calculation for the coe cient. For example, in the dynamic eddy-viscosity model, Cs = 12 hLijMijihM ijMiji = 12 ILMI MM (2.35) where ILM = Z t 1 Lij(t0)Mij(t0)W(t t0)dt0 (2.36) IMM = Z t 1 Mij(t0)Mij(t0)W(t t0)dt0 (2.37) An exponential weighting function is usually chosen [80]: W(t) = T 1 exp( t=T) (2.38) 36 where the time scale, T, controls the memory length of the Lagrangian averaging and here is set to T = 1:5 ( 8ILMIMM) 1=8 (2.39) As an exponential weighting function is used, ILM and IMM are actually solutions to relaxation-transport equations. Meneveau et al. [80] suggests that the computa- tional cost of the solution of the two additional transport equations can be reduced by approximating them as In+1LM (x) = Hf Ln+1ij Mn+1ij + (1 )InLM(x un t)g (2.40) In+1MM(x) = Hf Mn+1ij Mn+1ij + (1 )InMM(x un t)g (2.41) where Hf g is the ramp function (Hfxg = x if x 0, and zero otherwise) and = t=T1 + t=T (2.42) Linear interpolation is used to evaluate the integrals at the position x un t. A thorough analysis of the Lagrangian averaging procedure is given in Meneveau et al. [80]. 2.3.4 Filtering operations In the present study, ltering of the large-scale governing equations is implicitly de ned by the second-order nite di erences on the numerical grid. Explicit ltering is performed at test lter level (for the dynamic method). Second-order top-hat 37 lters were used and were calculated using fi = fi + h 2 f 24 2f x2 = fi + h2f 24 f i 1 2fi + fi+1 h2 (2.43) where h is the grid width and hf is the lter width. This yields the well known Simpson?s rule and trapezoidal rule lters with lter widths of hf=h = 2 and p6 respectively. Trapezoidal rule lters were used in this work, comprising a test lter, bf, of width p6 , where is de ned as geometric mean in each direction?s lter width: bfi = 1 4 (fi 1 + 2fi + fi+1) (2.44) and a grid lter, f, of width p3 : fi = 18 (fi 1 + 6fi + fi+1) (2.45) Trapezoidal rule lters were selected over Simpson?s rule based lters due to slightly improved results in LES of turbulent channel ow using the dynamic mixed model. 2.4 RANS equations and Reynolds stress modeling The equation for the RANS model are formally identical to the (2.6) where the term ij ij = uiuj uiuj (2.46) are de ned as the Reynolds stresses and must be modeled. In this dissertation an eddy-viscosity model is used to close the equations (2.6): ij ij3 kk = 2 tSij (2.47) 38 with t evaluated through the one equation Spalart-Allmaras model [70], or through a two equations k model. In this section, the two basic models will be shown: the Spalart-Allmaras [70] model and the two-layer k model [71, 72]. 2.4.1 Spalart-Allmaras model The Spalart-Allmaras model relies directly on a transport equation for the turbulent viscosity. The form of this model leads to the following equations: @ @t +~u r( ) = cb1 (1 ft2)S + 1 (r ( + )r ) +cb2 (r )2 cw1fw cb 2ft2 d 2 + ft1Du2 (2.48) The turbulent viscosity t is computed from t = fv1 (2.49) fv1 = 3 3 + C3v1 (2.50) = (2.51) where the modi ed viscosity equals t away from the wall. The damping function fv1 is based on the well know logarithmic law of the wall and lets t go smoothly to zero near the wall. S is related to the modi ed magnitude of the vorticity: S = jSj + 2d2 fv2 (2.52) fv2 = 1 1 + f v1 (2.53) 39 where jSj is the magnitude ofm the vorticity and d the distance to the closest wall. The function fv2 is built as fv1 on the hypothesis of a classical logarithmic-layer behaviour. is the von Karman constant 0.41. The function fw recovers the decay of the destruction term in the outer layer and produces a realistic skin-friction coe cient: fw = g 1 + c6 w3 g6 + c6w3 1 6 (2.54) g = r + cw2 r6 r (2.55) r = S 2d2 (2.56) The function ft1 and ft2 makes it possible to prescribe transition to turbulence: ft1 = ct1gte ct2 ! 2t Du2 (d 2+g2t d2t) (2.57) ft2 = ct3e ct4 2 (2.58) where dt is the distance between the current point and the transition to turbulence point, !t is the wall vorticity at the transition point, Du is the di erence between the velocity at the current point and that at the transition point, and gt is gt = min 0:1; Du! t xt (2.59) where xt is the grid spacing along the wall at the transition point. The model constants have the following default values: cb1 = 0:1355, cb2 = 0:622, = 0:6667, cv1 = 7:1, cw1 = cb1 + 1+cb2 , cw2 = 0:3, cw3 = 2:0. 40 2.4.2 k model We tested the two-layer k model [71, 72] where in the outer layer, the standard k model is used; the equations for k and in this case are: @k @t + r (~uk) = r + t k rk + Pk (2.60) and for its dissipation rate is @ @t + r (~u ) = r + t r + C 1f 1 kPk C 2f 2 2 k ; (2.61) with the term Pk, the production of turbulent energy, de ned as Pk = huiuji@ui@x j : (2.62) and the turbulent viscosity: t = C f k 2 ; (2.63) with f = e 2:5=(1+Rt=50) (2.64) f 1 = 1:0 (2.65) f 2 = 1 0:3e R2t (2.66) as given function of the turbulent Reynolds number Ret = k 2 (2.67) and C = 0:09, C 1 = 1:55, C 2 = 2, k = 1, = 1:3 41 In the inner-layer, the eddy viscosity is computed using t;inner = C ? p k (2.68) where the length scale, ? is given by ? = yc? 1 e Rey=A ; Rey = y pk (2.69) where y is the distance to the nearest wall and A = 70, c? = C 3=4 . The eddy viscosity is obtained by combining the inner- and outer-layer eddy viscosities through a blending function, , which uses a hyperbolic tangent to vary smoothly from 0 to 1: t = t;outer + (1 ) t;inner; (2.70) e = 12 1 + tanh Re y 200 0:1Rey= tanh(0:98) : (2.71) 2.5 Synthetic turbulence generation The synthetic turbulence generation method of Batten et al. [31] is used to create a three-dimensional, unsteady velocity eld at the in ow plane of the LES region. The Batten method is easy to implement numerically and takes into account the anisotropy of the ow with a very simple solution that will be explained in the following. An intermediate velocity, vi is rst constructed, using a sum of sines and cosines with random phases and amplitudes: vi (xj; t) = r 2 N NX n=1 h pni cos b dnj bxnj + !nbt + qni sin b dnj bxnj + !nbt i ; (2.72) 42 where bxj = 2 xj=Lb; bt = 2 t= b; (2.73) are spatial coordinates normalized by the length- and time-scale of the turbulence. In the above, b = k=" and Lb = bVb are the turbulence time- and length-scales, and Vb = k1=2 is the velocity scale.The random frequencies !n = N(1; 1) are taken from a normal distribution N( ; 2) with mean = 1 and variance 2 = 1. The amplitudes are given by pni = ijk nj dnk; qni = ijk nj dnk (2.74) where ni ; ni = N(0; 1), and bdnj = dnj V cn : (2.75) are modi ed wavenumbers obtained by multiplying the wavenumbers, dni , by the ratio of the velocity scale Vb = Lb= b to cn, given by cn = s 3 2hu 0 lu0mi dnl dnm dnkdnk : (2.76) The wave-numbers dni = N(0; 1=2) are chosen from a normal distribution with vari- ance 1=2, resulting in a three-dimensional spectrum that behaves like d4 exp( d2). Although the wave-numbers dni are distributed isotropically in a sphere, dividing them by cn tends to elongate those wave-numbers that are most closely aligned with the largest component of the Reynolds-stress tensor, and contract those aligned with the smaller ones. This results in a more physically realistic spectrum of turbulence, 43 with eddies that (near the walls) are more elongated in x, and tend to be more spherical in the channel center. The synthetic turbulent uctuation eld is nally reconstructed by a tensor scaling: u0i = aikvk (2.77) where aik is the Cholesky decomposition of the Reynolds-stress tensor. To better understand the e ects of the decomposition (2.2) in the generation of synthetic turbulence, it can be written as: vi(xj; t) = 2 N 1 2 NX n=1 < h (pni jqni ) cos b dnj cxnj + !nbt + jsin b dnj cxnj + !nbt i (2.78) where < is the real part and j is the imaginary unit. From (2.78) we obtain vi(xj; t) = 2 N 1 2 NX n=1 < h (pni jqni ) ej(cdnj cxnj +!nbt) i ; (2.79) expanding the exponential in Taylor series gives vi(xj; t) = 2 N 1 2 NX n=1 < 2 4(pni jqni ) 0 @1 + 1X m=1 jm b dnj cxnj + !nbt m m 1 A 3 5 (2.80) where b dnj cxnj + !nbt = bdnj 2 xj Lb + c! n2 t b = bdnj 2 xj bpk + c!n2 t b (2.81) with b = k= . The wavenumber and frequency modulation depends on the time scale b that is large near the wall and small far from it. As we can see from (2.80), the in-phase and in-quadrature decomposition (2.2) has all the powers of 44 b dnj cxnj + !nbt . This results in an energetic non linear modulation (in wave-number and frequency) that spreads the energy of the signal on to a broad spectrum. Far from the wall, however, b = k= is large with respect to its value near the wall, so that in (2.80) we can neglect the non linear term, resulting in a smooth linear modulation (uniform distribution of the wave number) of the signal energy. In this way, the spectra module of the turbulence is mimicked by the decomposition (2.2). Figure 2.1 shows the distribution of the wavenumbers dni and modi ed wavenum- ber dnicn for di erent distances from the wall. It is clear that far from the wall the distribution of the modi ed wavenumber becomes more isotropic in space. The result of the synthetic uctuation produced by the Batten?s decomposition (2.2) for di erent distance from the wall is shown in Figure 2.2. It is evident that close to the wall the structure are strongly anisotropy, stretched in the direction of the larger Reynolds stress. Far from the wall, on the other hand, the distribution becomes more isotropic. In the present work either the Spalart-Allmaras (S-A) or k model [71, 72] is used in the RANS section of the domain. From one of the two models the time scale and Reynolds stresses for the synthetic turbulence generation are computed. The Spalart-Allmaras model gives the averaged ow eld, and the turbulent viscosity t, while the k model gives the turbulent kinetic energy k, and the dissipation rate other than the averaged eld, of course. From the k model we have: t = C k2= (2.82) 45 where C = 0:09; therefore, the Reynolds shear stress hu0v0i are: hu0v0i = t @u@y (2.83) and the normal Reynold stresses are assumed to be equal: hu0u0i = hv0v0i = hw0w0i = 23k: (2.84) Therefore, the wavenumbers used in the present hybrid RANS/LES formulation are also isotropic. The time-scale for the generation of synthetic turbulence is the ratio of k over : b = k (2.85) If the Spalart-Allmaras model is used, to relate the TKE, k, to the Reynolds shear stress, hu0v0i, we use the experimental result [73, 74] j hu0v0ij = t @u @y = a1k (2.86) where a1 = pC and C is typically 0.09 as in the previous case. The normal Reynolds stresses are evaluated as in the S-A model, and for the time-scale, we use the de nition of eddy-viscosity from the k turbulence model to express in terms of k and t [75]: = c k2= t : (2.87) so the time scale b = k = 1p c @U@y (2.88) 46 In the simulations presented here, the number of modes, N, used to generate the synthetic eld was 200. This number was required to ensure that the resulting statistics were independent of the number of modes used. Further details of the method can be found in [31]. 47 ?2 ?1 0 1 2 ?2 ?1 0 1 2 . ?40 ?20 0 20 40?40 ?20 0 20 40 . . ?20 0 20?30 ?20 ?10 0 10 20 30 . ?100 ?50 0 50 100 ?100 ?50 0 50 100 d1 bd1 bd1 bd 1 d2 bd2 bd2 bd2 Figure 2.1: Distribution of the wavenumbers dni (top right) and modi ed wavenum- bers dnicn for di erent distances from the wall: top right gure is for y+ = 12, bottom left is y+ = 100, and bottom right is y = 5 48 x z 20 40 60 80 100 120 10 20 x z 20 40 60 80 100 120 10 20 x z 20 40 60 80 100 120 10 20 Figure 2.2: Synthetic uctuation produced by the Batten?s decomposition for dif- ferent distance from the wall: bottom gure is for y+ = 12, central is y+ = 100, and top is y = 5 The distance are based on a LES calculation from which the Batten decomposition is evaluated. 49 Chapter 3 The controlled forcing method and its tuning 3.1 Introduction The controlled forcing method for the turbulence generation consists of two parts: a method to generate synthtic turbulence at hte in ow (the method by Batten et al. [31] is used here) and the addition of a forcing term to the wall-normal momentum equation that ampli es the velocity uctuations in that direction, thus enhancing the production term in the shear-stress budget. The forcing amplitude is determined by a PI controller, which has two components acting in concert: a proportional part and an integral one; the two outputs are added together to form the actuating signal of the system to control, the Navier-Stokes equation system. PI controllers give a robust performance over a wide range of operating con- ditions and are widely used due to their easy implementation. They are frequently used when the mathematical model of the system to control is complex or unavail- able. Typically, the PI tuning for complex systems is based on semiempirical rules which have been proven in practice (for example the Cohen Coon tuning method [76] when a rst-order approximation of the system to control is available). In other cases, for unknow non-linear system for instance, the primary target in the PI tun- 50 Figure 3.1: Block diagram of the PI controller. ing is the stabilization of the controlled system and, marginally, the steady state performances (e.g., minimize the steady state error). In the controlled-forcing method, the system to control is the ltered NS equa- tion system with a dynamic subgrid model. The desired output is the time-averaged Reynolds shear stress, generally obtained from the RANS solution. The output of the system is the instantaneous u0v0 correlation; an exponentially weighted moving- average lter is added to obtain smooth Reynolds stresses. The block diagram of the resulting system is shown in Figure 3.1. First of all, we describe how the PI controller works in the closed-loop system using the schematic in Figure 3.1. The variable e(t), the di erence between the desired output hu0v0ides and the actual output hu0v0it, is sent to the PI that computes the integral of the error. The output of the PI (the control signal) is equal to KP times the magnitude of the error, plus KI times the integral of the error. The control signal is sent to the system to control, and a new output u0v0 is obtained. The new output is sent back (through the moving average lter) to obtain the new 51 error signal e(t). This process goes on and on. The error is de ned as e(y; z; t) = hu0v0ides(xo; y) hu0v0it(xo; y; z; t) (3.1) where hu0v0ides(xo; y) is the target Reynolds shear stress at the control plane x = xo, which is obtained from the RANS solution, and hu0v0it(xo; y; t) is the Reynolds shear stress or production, averaged over some time interval. The magnitude of the forcing is set to f(xo; y; z; t) = r(y; z; t) u(xo; y; z; t) huit(xo; y; z; t) (3.2) where r(t) = KP e(t) + KI Z e(t0)dt0 (3.3) is the output of the PI controller. The signi cance and optimal values of the two constants, KP and KI, will be discussed later. The forcing so de ned is added to the y-momentum equation. Enhancing the v0 uctuations through events with large u0 [the term in square brackets in (3.2)] has the e ect of accelerating the production of either Reynolds shear stress or the production of turbulent kinetic energy (TKE). The error is evaluated and forcing is performed at several planes within a few integral scales of the in ow (see below). In the next section we will analyze and the time-averaging window lter to relate it to local time scales of the ow. Then, the PI controller will be tuned in order to have stable results and minimize the error in a short distance downstream the controlled region. 52 3.2 Analysis of the controller parameters The controller described above has three parameters: the constants KP and KI and the width Tave of the averaging window used to obtain hu0v0it. In past applica- tions [49, 50, 30] no investigation was carried out in which they were systematically varied. In this section we examine the transient and steady-state responses of the ow to these parameters, and relate them to physical properties of the ow. The target is to choose them to obtain statistically steady-state Reynolds stresses that match the desired ones in a short distance without destabilizing the Navier-Stokes equation system. While we cannot claim that to have found true optimal values of these parameters, we will call as \optimal" the value that gives realistic turbulence in the shortest distance, and with the shortest transient. In the next subsection it is presented the test case where the controller will be tuned. 3.2.1 Testing strategy: ZPG boundary layer We test the PI controller in a at-plate, zero-pressure-gradient (ZPG) bound- ary layer. In general the strategy to implement a RANS/LES calculation is the following (see Figure 5.1): rst we perform a reference LES of the entire domain of interest; then we perform a RANS calculation of the equilibrium region, and use the RANS statistics to assign the in ow of the LES. Comparisons between the hybrid RANS/LES and the reference calculation allow us to evaluate the e ectiveness of the method. However, to test the PI control we will perform an LES calculation 53 Reference LES U (x), W (x) 88 LES U (x), W (x) 88 RANS RANS/LES interface Overlap region Control planes Figure 3.2: Sketch of the geometric con guration. instead of a RANS calculation: rst, we carried out the reference LES on a domain 240 o 25 o 25 o in the streamwise, wall-normal and spanwise directions, respec- tively; here o is the displacement thickness at the in ow. From this calculation, the time average Reynolds stress and the turbulence scales were extracted at x= o = 80. At this section a LES of dimension 120 o 25 o 25 o with the synthetic turbulence generation and the controlled forcing at the in ow began; in this way it was possi- ble to analyze the e ect of the KI, KP and Tave parameters without the in uence introduced by a RANS model. The control planes were distributed over a length of 27:5 o downstream of the in ow of the LES domain, with a spacing between planes of one boundary-layer thickness to allow the ow to re-adjust after the forcing (test calculations in which the forcing was distributed continuously were also carried out, with no signi cant di erences). The in ow Reynolds number (based on freestream velocity Uo and displacement thickness) was Re = 1000. The grid spacings in the 54 streamwise and spanwise directions were x= o = 1:25 and z= o = 0:385, while 64 points were used in the wall-normal direction (with y+min 1:0). The recycling and rescaling method of Lund [102] was used as in ow boundary condition for the reference LES. 3.2.2 The moving-average exponential window We begin by discussing the moving-average (MA) lter. In the classical termi- nology used in the feedback control theory, this block is called \Measuring Device"; the output of the NS equations are the three instantaneous velocity uctuations u0i, which are the inputs to the MA block. Its output is an appropriate time-averaged Reynolds stress that can be compared to the result of the RANS. We use an expo- nentially weighted moving-average lter, which places more emphasis on the most recent data available. At time-step k we use the following expression: hxki = 1 tT ave hxk 1i + tT ave xk (3.4) t is the time-step, and x is the data to average. As new data is accumulated, the contribution of the oldest data decreases. The value of t=Tave determines the memory of the lter. We de ne an attenuation time Td as the time after which the contribution of the signal at some time to the averaged data is negligible (say, 5% of the original contribution). It is easy to show 55 0 400 800 1200 1600 20000 1 2 3 4 x 10 ?3 Domain?averaged TKE t U o/d * o 0 50 100 150 200 250 6 8 10 x/d*o C f ? 10 3 (a) (b) Figure 3.3: (a) Domain-averaged turbulent kinetic energy for Tave = 1; KI = 5; KP = 30; Tave = 10; KI = 5; KP = 30; Tave = 100; KI = 1; KP = 1; Tave = 100; KI = 5; KP = 30. (b) Cf for the reference LES ( ) and the cases shown in (a). 56 that, in this case, Td Tave = t Tave log 0:05 log 1 tTave (3.5) since 1 >> t=Tave, (3.5) gives Tave ? Td=3. We expect that Td should be of the order of a large-eddy turnover time (LETOT) =u , which results in Tave ? 10 in the present calculation. Values much larger than this result in long transients, as the error a ects the input to the forcing for a long time, and the forcing does not adjust rapidly enough to the present ow conditions, but is strongly a ected by past events. This is shown in Figure 3.3(a) where the domain-averaged turbulent kinetic energy (TKE) is plotted as a function of time; in all the calculations we used the optimal values of KP and KI, except for one case in which we employed the values used in [50, 30], that is KP = KI = 1. A longer transient can be observed in the case Tave = 100 compared with Tave = 10. With the lower values of KP and KI used in early applications of this method the transient is still quite long, the magnitude of the oscillations is signi cant, and reasonable ow statistics require very long averaging times. An optimal choice of KP and KI, on the other hand, reduces the amplitude of the uctuations, so that even an incorrect value of Tave is less damaging. A low value of Tave results in incorrect prediction of the ow: since the short averaging gives values of hu0v0it that are too far from a stationary sample, and thus may result in excessive forcing even if the ow has reached a realistic state. In Figure 3.3(b) the skin-friction coe cient Cf = 2 w= U21 (where w is the wall stress) is shown for the same two cases, and it appears that the steady state 57 0 50 100 150 200 2503 4 5 6 x/d*o C f ? 10 3 Figure 3.4: Skin friction coe cient Cf with KI = 5 and KP = 30; rst control plane after one grid cell; rst control plane after one boundary layer thickness results are also strongly in uenced by the time window of the MA lter. Excessively small values of Tave result in incorrect steady state results, while the other values of Tave are in reasonable agreement with each other. The case with Tave = 10 (i.e., Td of the order of one LETOT) seems, altogether optimal, since it results in somewhat more accurate statistics and a much shorter transient (and reduced CPU costs). Unless explicitly speci ed, Tave = 10 for all following cases. On a related note, it should be mentioned that previous applications of the controlled forcing method [49, 50, 30] used spanwise averaging (in addition to the time-averaging in the MA block) to supply a smoother signal to the forcing; also the rst control plane was located one boundary-layer thickness downstream of the in ow plane. We now try to better address these choices. Regarding the control plane, the result in term of Cf coe cient in Figure 3.4 shows that the best choice is to have the rst control plane immediately after the in ow plane. This was not evident in [50, 30] because the parameters were not in the optimal range, expecially 58 0 50 100 150 200 2503 4 5 6 x/d*o C f ? 10 3 Figure 3.5: Skin friction coe cient Cf with KI = 5 and KP = 30; hu0v0it is spanwise-averaged; hu0v0it is not spanwise-averaged. the time window Tave. If this is not tuned well, the control has a weak sensitivity on the position of the control planes. Regarding the spanwise average, Figure 3.5 shows calculations in which no spanwise averaging was performed, compared to a simulation in which the u0v0 cor- relation is averaged in the spanwise direction as well as in time. Better result are obtained when the error re ects more closely the local and instantaneous conditions of the ow. We observe more rapid adjustment of the ow towards the reference LES, as the forcing is more responsive to the local state of the ow. In the next subsection we will turn the attention to the PI controller. This will allow us to evaluate the in uence of its parameters on the transient time and steady-state results. 59 3.2.3 The PI control In general a proportional controller a ects the rise-time and the transient time of the dynamic system: it modi es the bandwidth of the frequency response of the closed loop system and the gain of the zero frequency. If KP is large, the output of the integral block is more sensitive to the high-frequency variation of the error, and the stability of the system can be a ected. The integral control, on the other hand, gives a large gain at low frequency and reduces the break frequency (the point where the amplitude spectrum of the transfer function is zero); its e ect is to eliminate the steady-state error and attenuate the high-frequency disturbances. In fact, the output of the integral control will change over time as long as an error exists, but the phase lag between the input and output of the integral block increases for all of the frequencies (the phase of the integrator block starts at 90o). If KI is large, the phase lag is extended to higher frequency, so the system can became oscillatory and potentially unstable. To better understand these considerations, we evaluate the transfer function of the PI control, i.e., the Laplace transform of (3.3) divided by the Laplace trasform of the error e(t): G(s) = KP + KIs (3.6) where s is the complex variable s = +j!. It is simple to see the frequency response of the PI control: setting s = j! (Fourier domain) we have that the magnitude and 60 the phase of G(!) are jG(!)j = r K2P + K 2 I !2 (3.7) = tan 1 KI!K P : (3.8) To see the behaviour of these function in the frequency domain, the Bode plot is used in gures 3.6 and 3.7 (i.e. jG(!)jdB = 20 log10(jG(!)j) vs log10(!) and vs Log10(!) ). From these gures is evident that the integral control introduce a phase lag that increases with KI Figure 3.8 shows how the transient is a ected by KP and KI. If KP = 0 (only the integral control is activated) the calculation diverges: in the initial transient the forcing is proportional to e(t0)dt0 and the initial time-step is small because of the initial unphysical ow in the computational domain. Because of the small time-step the integral action is weakened, so the integral control needs a longer time before it can give a correct control signal. As the simulation advances, however, the time-step is further decreased because no fast correction is present in the domain, so the correct integral action is even more delayed until the simulation goes unstable. The case of high KP is the opposite: as soon as the simulation starts there is the fast correction due to the proportional control. Apart from the length of the transient, however, the steady-state results did not di er much in all the converged calculations, showing little sensitivity of the ow to the proportional controller (at least for KP > 1). Similarly, we found instability for KI 30, as can be expected by the theory of the PI feedback control. If KI was close to zero, only the proportional control 61 10?2 10?1 100 101 102 20 40 60 80 100 Magnitude (dB) w (rad/sec) 10?2 10?1 100 101 102 ?80 ?60 ?40 ?20 0 w (rad/sec) Phase (deg) KI KI Figure 3.6: Bode plot of the frequency response of the PI system for di erent value of KI and xed KP = 30. KI = 5; KI = 20; KI = 30; It is evident that when KI increase, the gain at low frequency increase and the phase lag value starting at 90 is kept to higher frequency, i.e., the response of the PI system lags behind the input wave by 90 deg in the worst case 62 10?2 10?1 100 101 102 0 20 40 60 80 Magnitude (dB) w (rad/sec) 10?2 10?1 100 101 102 ?80 ?60 ?40 ?20 0 w (rad/sec) Phase (deg) KP KP Figure 3.7: Bode plot of the amplitude and phase frequency response of the PI system for di erent values of KP and xed KI = 5. KP = 5; KP = 20; KP = 30. When KP increases, the gain at high frequency increase and the phase lag value decreases at higher frequency 63 0 500 1000 15000 .001 .002 .003 .004 Domain?averaged TKE t Uo/d*o 0 500 1000 15000 .001 .002 .003 .004 Domain?averaged TKE (a) (b) Figure 3.8: Domain-averaged turbulent kinetic energy. (a) KI = 5 and KP = 0; 15; 30; 500; (b) Kp = 30 and KI = 1, 5, 20, 30. 64 0 200 400 600 800 1000?5 0 5 10 x 10 ?3 e(t) t Uo/d*o 0 200 400 600 800 1000?0.2 ?0.1 0 .0 0.1 0.2 t Uo/d*o Integral e(t?)dt? (a) (b) Figure 3.9: (a) Instantaneous and (b) integrated error for KP = 30 and KI = 5 (lines) and KI = 20 (lines with symbols). y+ = 13; y= o = 4:5. The curves for KI = 20 are shifted upwards by 0.005 in (a) and 0.1 in (b). 65 was activated and the steady state error (the di erence between the desired Reynold stress and the one obtained from the calculation) was large. Figure 3.9 shows the error and its integral at two locations, one in the wall layer, the other in the outer region of the boundary layer, for KP = 30 and two values of KI. Increasing KI reduces the amplitude of the dominant frequency of the error, as well as its ampli- tude, at least in the outer layer. Note, however, that the forcing signal (which has the integral error multiplied by KI) does not change its magnitude. 3.2.4 Conclusions As a result of the tests described in this section, we conclude that some local- ization of the error, both in time and space, is desirable. Compared with the work of Keating et al. [30], we obtained improved results using a shorter time-averaging window, and removing the spanwise averaging. Averaging over a time interval of the order of the local integral time-scale resulted in shorter transient and more rapid development of physically realistic turbulent eddies. Removal of the spanwise aver- aging used in previous investigations also gave improved results, as the error (and hence the forcing) re ected more accurately the local state of the ow. The con- troller parameters, the constants KP and KI, play a less signi cant role. Giving excessive weight to the integral error, or insu cient one to the proportional part, resulted in instabilities of the ow. For a wide range of values of KP and KI, on the other hand, the ow statistics were found to be insensitive to the constants. We 66 observed that the controller coe cients mostly a ect the length of the transient; for a wide range of both KI and KP , however, the ow statistics are fairly insen- sitive to the parameter values. In the next section we will apply the method with the improved coe cients to four ows, to determine its performance (with the new coe cients) in actual cases. In the next Chapter, we will report results obtained applying the PI control to di erent boundary layer cases, including boundary layers undergoing adverse and favorable pressure gradients, and a three-dimensional boundary layer. 67 Chapter 4 Application of the controlled forcing 4.1 Introduction In this Chapter the application of the controlled forcing is presented for di er- ent boundary layer con gurations: zero pressure gradient (ZPG), adverse pressure gradient (APG), favorable pressure gradient (FPG), and three dimensional bound- ary layer (3D). It is important to test the controlled forcing for these di erent ow, to prove that the tuning discussed setup in the previous Chapter, especially the optimal range found for the PI control, is not a \single-case" optimal tuning, but can be extended to di erent con guration. The only parameter that need to be matched with the local characteristic of the ow is Tave, the width of the moving average lter that must match with the intrinsic scales of the turbulent structures. 4.1.1 Zero-pressure-gradient boundary layer In this Section we report results obtained from calculations of the zero-pressure- gradient (ZPG) boundary layer. In the previous Chapter the reference LES provided the in ow statistics; here, they are calculated from the RANS; the modeling errors that would a ect real hybrid RANS/LES calculations are going to play an important 68 role. Figure 4.1 compares a case in which the k model is used for the RANS with one that employs the SA model. One can observe the e ect of modeling errors: First, the k model gives a prediction of the wall stress in better agreement with the reference LES than the SA model; the k model also gives k and (which are required by the synthetic turbulence generation method) directly, while with the SA model they must be estimated using (2.86). The skin-friction coe cient matches the reference LES value after only 20 o from the end of the control region (approximately 2 boundary layer thicknesses) with the k model, while it requires a longer distance with the SA model. Both the mean velocity pro les and the Reynolds stresses are in good agreement with the reference calculation, except perhaps in the outer layer where the growth of realistic structures is somewhat slower (as was discussed in [50]). Figure 4.3 shows the Cf obtained with the coe cients Tave = 100, KI = KP = 1 as in [50], and the case when only sysnthetic turbulence is applied at the in ow without any control. One limitation of the current approach lies in the fact that the error is based on the Reynolds shear stress, which is not a coordinate-invariant quantity; in com- plex geometries, the error might be ill-de ned. The proposed RANS/LES merging approach, however, is predicated on the RANS calculation being accurate in the interface region, which in practice restricts the interface to lie in a thin, attached shear layer, where the de nition of directions along and normal to the shear layer is 69 10?2 10?1 100 101 0 2 4 6 8 x 10 ?3 y / d*o ? / U 2 o 100 101 102 103 0 10 20 30 40 y+ U + 0 50 100 150 200 2503 4 5 6 7 x / d*o C f ? 10 3 (a) (b) (c) Figure 4.1: Zero-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the loca- tions indicated by a vertical line in part (a). Reference LES; k RANS; Hybrid k /LES; SA RANS; Hybrid SA/LES. 70 10?2 10?1 100 101 0 2 4 6 8 x 10 ?3 y / d*o ? / U 2 o 100 101 102 103 0 10 20 30 40 y+ U + 0 50 100 150 200 2503 4 5 6 7 x / d*o C f ? 10 3 (a) (b) (c) Figure 4.2: Zero-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the loca- tions indicated by a vertical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. 71 0 50 100 150 200 250 2 2.5 3 3.5 4 4.5 5 5.5 6 x / d*o C f ? 10 3 Figure 4.3: Skin friction coe cient Cf in the case of hybrid k /LES simulation. Reference LES; RANS k ; Hybrid k /LES with Tave = 10, KI = 5, KP = 30; Hybrid k /LES with Tave = 100, KI = 1, KP = 1 as in [30]; Hybrid k /LES with only Batten synthetic turbulence at the interface 72 unique. Despite this practical observation, it would be desirable to develop a more universal error de nition, and in particular one that is invariant to the choice of the frame of reference. One possible choice to satisfy this requirement is to base the error on the production of turbulent kinetic energy, which satis es the desired invariance properties. Thus, we de ne a production-based error eP (y; z; t) = hu0iu0jideshSijides + hu0iu0jithSijit (4.1) (where the dependence on xo; y; z and t has been omitted, and Sij is the strain- rate tensor). For the ZPG boundary layer examined here this de nition of the error should give the same result as (3.1), since the shear stress is the leading term in the de nition of the production. We found results comparable to those obtained using the shear stress in the de nition of the error, with two notable di erences: rst, eP decays more rapidly than e away from the wall (since the velocity gradient goes to zero); this results in less forcing in the outer layer, so that the correct Reynolds- stress pro le away from the wall is established further downstream, compared with the case in which the error is de ned by (3.1). Secondly, we observed a much longer transient before the error stabilized; may be due to the decreased forcing in the outer layer again, which slows the production of hu0v0i, and hence the decrease of the error there. 73 0 50 150 250 3500 0.2 0.4 0.6 0.8 1 x / d*o U ? / U o Figure 4.4: Freestream velocity for the adverse-pressure-gradient calculation. 4.1.2 Adverse-pressure-gradient boundary layer Hybrid simulations were next carried out in an adverse-pressure-gradient (APG) boundary-layer that undergoes separation. The con guration of these simulations is similar to that of [30], with an in ow Reynolds number, Re = 1260, and a pro le of V1 imposed at the top boundary that results in a signi cant deceleration of the ow (see Figure 4.4). Due to the strong adverse pressure-gradient, the ow separates; a favorable pressure gradient then closes the recirculation bubble. The computational domain used in the full-domain LES was 380 o 64 o 20 o. This LES, which used the rescaling/recycling method at the in ow, had a grid of 384 192 64. The RANS calculations for the hybrid cases were extended to cover the entire ow domain to remove di culties of placing out ow boundary 74 conditions close to reversed- ow regions; both SA and k models were used for the RANS solution. The LES region started at x= o 80, and used the RANS data at that location for the synthetic turbulence generation, as well as target for the controlled forcing algorithm. The =u obtained from the k model at the interface location was approximately 110, which agrees well with the value from the full LES. To match this value, Tave was xed to 38, so that the memory Td was approximately 113. KI was xed at 5 and KP to 30. Figure 4.5 shows the Cf, mean velocity pro les and Reynolds shear stresses at three locations. As the boundary layer is subjected to the adverse pressure- gradient, Cf decreases until the ow separates at x= o 170. A weak separation bubble, which has a length of approximately 80 o, can be observed in the full-domain LES. The results from the RANS simulation in terms of Cf are close to those of the full LES. Switching to the LES is, however, bene cial (the prediction of Cf and the dimensions of the separation bubble are both more accurate). Shortly after separation (the second pro le in Figure 4.5) some di erences between the hybrid cases can be observed, whereas inside the separation bubble the LES give similar results. In this ow, the ampli cation of turbulence in the separated shear layer acts as a powerful mechanism to accelerate the generation of realistic eddies. However, the mean streamline shown in Figure 4.6 show that the separation bubble in the hybrid calculation is still larger then the full LES, especially when the SA model is applied in the RANS region. As the zero pressure gradient boundary layer, the 75 10?2 10?1 100 101 0 2 4 6 x 10 ?3 y/d*o ? / U 2 o 10?2 10?1 100 101 102 0 1 2 3 y/d*o U/U o 0 50 100 150 200 250 300 350 400 0 2 4 6 x/d*o C f ? 10 3 (a) (b) (c) Figure 4.5: Adverse-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; k RANS, SA RANS; hybrid k /LES; Hybrid SA/LES. 76 k model gives better performances in the hybrid simulation Keating et al. [30] observed a much stronger dependence of the hybrid results on the accuracy of the RANS model (see Figure 11 in [30]), as is shown in gure 4.7. This e ect is due to the use of suboptimal PI controller parameters and averaging window. With the more localized averaging used here, the initial decay of the wall stress that was observed by Keating et al. [30] does not occur. Figure 4.8 shows contours of streamwise velocity uctuations in a plane near the wall for the reference LES and a hybrid case. We observe a very rapid devel- opment of a physical streaky structure. The onset of separation (which is strongly a ected by the ow state immediately before the separation occurs) is predicted well, and the length-scales of the ow in the separated- ow region are remarkably similar, con rming again the robustness and e ectiveness of the method. 4.1.3 Favorable-pressure-gradient boundary layer Simulations of a boundary layer subjected to a favorable pressure-gradient (FPG) were performed. In this ow, if the acceleration is rapid enough, re-laminarization and re-transition of the ow may occur. The case studied here matches the ex- periments of [97]. The freestream acceleration (from U1 = 1 to U1 3) be- gins at approximately 100 o and is completed by 450 o. The ow acceleration is achieved by imposing a streamwise velocity pro le U1(x) at the top boundary of the domain [102]. Its magnitude was calculated from the acceleration parameter, 77 Figure 4.6: Mean streamlines for the APG boundary layer. (a) Reference LES; (b) k RANS; (c) SA-RANS; (d) hybrid k /LES; (e) hybrid SA/LES. 78 0 50 150 250 350 400 0 2 4 C f ? 10 3 x/d*o 0 50 150 250 350 400 0 2 4 x/d*o C f ? 10 3 (a) (b) Figure 4.7: Skin friction coe cient Cf in the case of LES/LES simulation (a) and hybrid RANS SA/LES (b). In (a) Tave = 38, KI = 5, KP = 30; Tave = 100, KI = 1, KP = 1 as in [30]. In (b) RANS SA model, hybrid SA/LES with Tave = 38, KI = 5, KP = 30; Tave = 100, KI = 1, KP = 1 as in [30]; Hybrid SA/LES with only Batten synthetic turbulence at the interface 79 z / d * o 100 150 200 250 3000 10 20 (a) x/d*o z / d * o 100 150 200 250 3000 10 20 -0.1 0 0.1 (b) Figure 4.8: Streamwise velocity uctuations in the APG boundary layer, in the y= o = 0:09 plane. (a) Reference LES; (b) hybrid k /LES; the red rectangle indicates the region where the control is active. K = ( =U21)(dU=dx), experimentally obtained from [97]. The acceleration parame- ter, K, and the free-stream velocity distribution U1(x) are the same as [30]. Since K exceeds the critical value for re-laminarization (Kcrit = 3:5 10 6) for an ex- tended region of the ow, turbulence is expected to be damped in the region of high acceleration. The ow then re-transitions once the acceleration is removed. A full-domain LES was rst performed using the rescaling/recycling method at an in ow Reynolds number, Re = 1260. For this simulation the domain length was 476 o, the width was 20 o and the height was 20 o. A somewhat coarse grid (512 64 64) was used, since the comparisons were primarily being made between this full LES and a hybrid RANS/LES simulations. One hybrid calculations was performed that included a RANS domain 350 o long, and an LES region that started at x= o = 225 and extended to 476 o (i.e., had a length 251 o). Synthetic turbulence with controlled forcing was introduced at the RANS/LES interface. 80 10?2 10?1 100 101 0 0.01 0.02 y/d*o ? / U 2 o 10?2 10?1 100 101 0 0.5 1 1.5 2 y/d*o U/U ? 0 100 200 300 400 500 2 4 6 8 x/d*o C f ? 10 3 (a) (b) (c) Figure 4.9: Favorable-pressure-gradient boundary layer. (a) Skin friction coe cient Cf (b) mean velocity pro les and (c) pro les of the Reynolds shear stress at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; hybrid SA/LES, Tave from LES data; hybrid SA/LES, Tave from RANS data; . 81 0 100 200 300 400 500 2 3 4 5 6 7 8 x / d*o C f ? 10 3 Figure 4.10: Skin friction parameter Cf for the favorable pressure gradient boundary layer: reference LES, RANS, hybrid RANS/LES with Tave = 2:4 and Tave = 1; Hybrid SA/LES with only Batten synthetic turbulence at the interface All the di erent RANS model tested give signi cant modeling error; in par- ticular, no model was able to predict the relaminarization and the re-transition of the accelerating boundary layer correctly. Furthermore, with the SA model, the timescale =u at the interface location was found to be approximately 2, while the value obtained from the full LES calculation is twice as large in the region where the control was applied. Two hybrid calculation were performed, one that used a value of Td close to the time-scale predicted by the RANS (giving Tave = 1), and another that matched the LES timescale (Tave ? 2:4). Figure 4.9 compares the results of the various simulations. The mean data at the interface, which is supplied by the RANS, is not accurate; the controller tries to drive the LES towards the RANS 82 region, but when Tave is small a smooth hu0v0i cannot be obtained, and excessive forcing is applied, resulting in an initial overshoot of the skin-friction coe cient. Once the control is released, the ow adjusts fairly rapidly towards the full-domain LES. When a longer average is used, the controller is more successful in driving the LES shear stress towards the desired distribution (Figure 4.9(c)) and the LES ad- justs more rapidly and with no overshoot. Note nally that this is a very challenging test case: RANS models cannot be expected to predict relaminarization accurately; furthermore, since the energy of the uctuations is decreasing because of the ac- celeration, adding energy through the forcing term is moving the ow away from its equilibrium. This is evident in Figure 4.10 where are shown the two cases with control presented before, and the case when only synthetic turbulence is used at the in ow: in this last case the ow is drown through the full LES con guration soon after the rst grid cell; however, the retransition to turbulent ow is leaded. De- spite the fact that this can be considered an o -design application for the controlled forcing, the results are altogether acceptable. 4.1.4 3-D Boundary layer Simulations were then performed on a 3D boundary layer obtained by applying a spanwise pressure gradient to a at-plate boundary layer. The pressure gradient increased from zero, at x= o=100, up to 1:5 10 3 at x= o=150, and remained constant thereafter. The magnitude of the pressure gradient was set in such a way 83 0 50 100 150 200 250 300?50 ?40 ?30 ?20 ?10 0 10 x/d*o a Wall Freestream Figure 4.11: 3D boundary layer. Flow angle = tan 1 U=W. Reference LES. SA-RANS; error based on hu0v0i; production-based error. as to turn the ow by 45o near the wall and 24o at the free stream by the end of the domain. Figure 4.11 shows the turning angle near the wall and in the freestream. Once again, two simulations were carried out: a full-domain LES, and one hybrid RANS/LES calculations. The computational domain of the reference LES was 291 o 20 o 20 o with a grid of 720 100 292. An in ow plane from a precursor simulation was used with an in ow Reynolds number of Re = 1260. Because of the turning of the ow, the grid had to be re ned in the streamwise direction in the downstream area. We de ne a local coordinate frame with in the ow direction and in the lateral direction in the wall plane: + = x+ cos + z+ sin ; + = x+ sin + z+ cos : (4.2) We maintained streamwise and spanwise grid resolutions of + < 50 and + < 13 as is shown in gure 4.12. Here, is the ow angle at the wall and wall units are 84 0 50 100 150 200 250 3000 10 20 30 40 50 60 x/d*o Grid resolution + Dz+ Dx+ Figure 4.12: Grid resolution obtained in the 3D boundary layers +; + de ned using the resultant wall stress u = ( 2xy;w + 2xz;w)1=4; xy;w = @U@y y=0 ; xz;w = @W@y y=0 : (4.3) In the hybrid cases, the RANS equations was solved using the Spalart-Allmaras model. The LES region started at x= o 200 where the ow angle was 32o near the wall and 12o in the free-stream. The SA model predicts a value of =u = 40 at the RANS/LES interface; consequently, we set Tave = 14 (corresponding to Td ? 42). We also performed simulations with the production-based error, in addition to those with the standard de nition of e. Figure 4.13 shows the skin friction coe cient and mean velocity pro les in the streamwise and spanwise directions. The development of Cf;x is similar to that of the ZPG case (with an initial decrease followed by a quick recover). For the lateral ow we observe remarkably good agreement with the reference simulation beginning from the end of the controlled region, approximately two boundary layer thicknesses 85 100 101 102 103 0 5 10 15 20 y+ W + 100 101 102 103 0 10 20 30 40 y+ U + 0 50 100 150 200 250 300?0.005 0 0.005 x/d*o C f Cf;x Cf;z (a) (b) (c) Figure 4.13: 3D boundary layer. (a) Streamwise and spanwise skin friction coe - cients, Cf;x and Cf;z; (b) mean streamwise velocity pro les and (c) mean spanwise velocity pro les at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. 86 100 101 102 103 0 0.5 1 1.5 2 2.5 y+ ? + 100 101 102 103 ?0.2 0 0.2 0.4 0.6 0.8 y+ ? + 0 50 100 150 200 250 300 ?0.005 0 0.005 x/d*o C f Cf;x Cf;z (a) (b) (c) Figure 4.14: 3D boundary layer. (a) Streamwise and spanwise skin friction coe - cients, Cf;x and Cf;z; (b) secondary shear stress hv0w0i and (c) principal shear stress hu0v0i at the locations indicated by a vertical line in part (a). Reference LES; SA RANS; error based on hu0v0i; production-based error. 87 (b) (a) Figure 4.15: Contours of streamwise velocity fuctuations in a xz plane at y= o = 0:18. (a) Reference LES (only part of the domain is shown); (b) hybrid RANS/LES. downstream of the in ow. The principal shear stress, hu0v0i and the secondary ones hv0w0i develop quite rapidly as is shown in 4.14. Notice that by adding the forcing into the v equation we amplify the production of the secondary stress, hv0v0i@W=@y, as well as that of hu0v0i. Figure 4.15 shows streamwise velocity uctuations in a plane parallel to the wall. We observe again a very rapid build-up of the streaky structure starting from the synthetic turbulence. 4.2 Conclusions In this Chapter we show further validation of the controlled-forcing method for turbulence generation at RANS/LES interfaces. The method has shown itself to be robust and e cient. By studying the zero-pressure-gradient boundary layer (which, despite its simplicity is one of the most challenging cases examined) we developed 88 guidelines for the determination of the parameters of this method. Perhaps the most important parameter is the length of the averaging used to determine the computed Reynolds stresses in the de nition of the error (3.1) or (4.1). Of the ows examined, we found the zero-pressure-gradient and favorable- pressure-gradient boundary layers to be the most challenging: in the zero-pressure- gradient case no external mechanism amplify the instability of the ow, so that the generation of turbulence is entirely due to the synthetic turbulence and controlled forcing; since the latter is less e ective in the outer layer (where the turbulence production mechanisms are weaker), we observe the slow development of the outer layer. The favorable-pressure-gradient boundary layer presents a di erent set of problems: The acceleration results in relaminarization of the ow. Therefore, the injection of energy due to the forcing drives the ow in the wrong direction. In addition, RANS models are not very accurate in the strong acceleration region, which results in three sources of error: rst, the mean velocity at the interface is far from the result obtained with the full-domain LES. Second, the target stresses are not close to the \correct" ones; thus, the controller drives the solution towards an incorrect one. Finally, the time-scale obtained from the RANS is substantially di erent from that obtained from the reference calculation, so that the moving average is suboptimal. The method worked well in the adverse-pressure-gradient case, in which the deceleration (and, to an even greater extent, the separation) ampli es the distur- 89 bances, thus aiding the generation of turbulent eddies. This result is consistent with the ndings of Terracol [25], who had some success using synthetic turbulence gener- ation in the adverse-pressure-gradient region on the suction side of a turbine blade. In the 3D boundary layer the turning of the streamline imposed by the spanwise pressure gradient drives the lateral ow, with secondary stresses that develop from nearly zero at the RANS/LES interface. The controller proves to be able to seed the turbulence enough that the growth of the secondary stresses, and hence the lateral mean velocity pro le, can be predicted accurately. In summary, we have shown that the controlled forcing method is robust and e cient. The least favorable results are obtained in mildly unstable ows like the ZPG pressure gradient, or in re-laminarizing ones. Even in these situations, within 5 boundary layer thicknesses of the end of the control region a realistic ow was established. In realistic applications, especially in curved ows and in adverse pres- sure gradients (an important class of applications) the controlled forcing gives very rapid development of realistic eddies. 90 Chapter 5 Favorable pressure gradient boundary layer 5.1 Introduction Turbulent boundary layers subjected to a favorable pressure gradient (FPG) (i.e., one that results in freestream acceleration) are common in many engineering applications, such as airfoils and curved ducts. While the canonical zero-pressure- gradient (ZPG) boundary layer is relatively well understood, FPG boundary-layers are less well known. The simplest case of an accelerating boundary layer is the \sink ow", in which the acceleration parameter K = ( =U21)dU1=dx is constant with the streamwise distance x. This ow has been studied experimentally and numerically; it is known that, for strong acceleration (K > 3 10 6) turbulence cannot be maintained, and the ow re-laminarizes. The sink ow is, of course, an idealization: in real cases, a large acceleration parameter cannot be sustained for long distances, and, in practical applications, a region of FPG and streamwise acceleration is followed by one with constant or adverse pressure gradient (such is the case for the ow on an airfoil downstream of maximum thickness). In these conditions, full re-laminarization may not occur, and the pressure gradient may leave the ow in a \laminarescent" [83] state. As the pressure gradient is removed, 91 the ow may then return to a turbulent state; the re-transitioning process may be strongly a ected by the state of the turbulence at the end of the acceleration region. For these reasons it is important to understand the mechanisms of relaminar- ization. Reviews of current knowledge can be found in several articles by Narasimha and Sreenivasan [84, 85] and by Sreenivasan [83]; here, only the main ndings are summarized. Re-laminarization can be caused by three mechanisms: (1) a decrease in the Reynolds number accompanied by an increase in the viscous dissipation; (2) strati cation or ow curvature, or (3) ow acceleration. Experimental investigations of re-laminarization due to ow acceleration started in the early 1960s. Wilson [86] found the rst evidence that the ow did not follow the semi-empirical theories for a fully turbulent boundary layer: the measured heat- ux rates on a convex surface of a blade were considerably lower than the predicted values. Wilson [86] conjectured the possibility of a \reverse transition" of the ow due to the low local value of the Reynolds number. However, Patel and Head[87] later observed that there is no explicit correlation between a low Re and the relaminarization process, as long as the initial Re is high enough to allow the turbulence to be self-sustained. Senoo [88] studied the boundary layer on the end-wall of a turbine nozzle cas- cade; he found that it was laminar in the region of the throat, while the upstream layer was turbulent. His ndings did not explain the phenomenon of relaminariza- tion. Moreover, the e ects of secondary ow were not clear in his study. Studying a similar con guration, Launder [89, 90] observed that relaminariza- 92 tion of the ow, at least in terms of macro-scale properties and integral parameters, starts near the end of the acceleration region. In the upstream region, the boundary layer was found be turbulent but with a departure from the universal law of the wall. Only further downstream the boundary layer became close to the laminar state. Moreover, he documented an increase of the energy at low wave-numbers in the near-wall region and an even more pronounced one in the outer part of the boundary layer. He found an explanation of that energy shift in the time-lag be- tween inner and outer regions: the latter is relatively distant from the turbulence production peak, so the eddies in that region are not fast enough to adjust their frequency as the free stream is accelerated. He, therefore, proposes a picture in which the ow dynamics are completely dominated by the near wall region. Laun- der [89, 90] also observed that the turbulence does not vanish, but an increasing fraction of it plays a passive role in the boundary-layer development. Because of this inability of turbulent structures to adjust to the ow acceleration and the rapid increase of momentum, the viscous stresses grow larger than the turbulent ones; con- sequently, the dissipation exceeds production, leading to a decay of turbulence and to laminarization. However, Badri-Narayan et al. [91] found that dissipation never exceeds production in an accelerated boundary layer, and that both production and dissipation decrease. Kline et al. [92] found a correlation between the drastic drop of eruptions in the bu er layer and the acceleration parameter K: for K that approach the value of 93 Kmax 3:5 10 6 the burst rate tended to zero. Narasimha and Sreenivasan[85] sug- gested that the bursting frequency decreases exponentially in an accelerated bound- ary layer before the ow relaminarizes, but never goes to zero completely. Patel and Head [87] showed a strong correlation between the distribution of the shear-stress gradient near the wall and the relaminarization. They de ned a non- dimensional shear-stress-gradient parameter, and they inferred the onset of relami- narization from the position where this parameter has a minimum value of 0:009. However, Narasimha and Sreenivasan [85] observed that this non-dimensional pa- rameter reaches that value before the turbulent boundary-layer actually reverts to a laminar-like state; in this sense, the Patel and Head [87] parameter predicts deviation from the universal law for ZPG boundary layers, but not necessarily relaminariza- tion. Blackwelder and Kovasznay [93] noted an increase of the Reynolds stress and turbulent kinetic energy (TKE) along streamlines in the outer part of the boundary layer as the ow was accelerated. In the inner layer, on the other hand, they found a decrease of those quantities along streamlines. Moreover, they studied the space-time correlations of the large structure in the outer layer. They found no change compared to the ZPG case for the uu correlations; this suggested that the acceleration had no strong e ect on the large streamwise structure (those close to the boundary-layer edge). On the other hand, they noted a drastic change of the normal velocity component of the outer-layer structures, as the vv correlations lost 94 the anti -symmetrical part found in the ZPG case. Narasimha and Sreenivasan [84] point out that relaminarization is the result of the domination of pressure forces over the slowly responding Reynolds stresses in the outer layer, and the generation of a new laminar sub-layer that is maintained stable by the acceleration. In their model the turbulent structures in the outer layer are only distorted and not destroyed by the rapid acceleration. The new sub-layer and the distorted outer layer do not interact, but they only provide the appropriate boundary conditions to each other. Rapid distortion theory could be applied in the outer part to predict the turbulence intensities there, but the interior part of the layer was not well understood and the wall-normal components were not well predicted. In the same period, Falco [94] performed a smoke-visualization study, and suggested that in the relaminarization process large scale structures existed upstream of the contraction that accelerates the ow: the boundary layer in the later stages of the acceleration is dominated by an array of large scale streamwise vortices; thus, an inner-outer layer interaction could exist and the relaminarization seems to begin from the outer region with a strong coupling between inner and outer parts. This visual observation was recon rmed by Ichimiya et al. [95], who conjecture that non-turbulent uid on the outside of the boundary layer penetrates near the wall, so that the beginning of relaminarization is due to the outer region together with the change in ejection-sweep phenomena. Recent experimental studies of FPG boundary layers have been performed by 95 Fernholz and Warnack [96] and by Warnack and Fernholz [97]. Their measurements showed that the Reynolds number had little e ect on the relaminarization, and the pressure-gradient e ects were dominant. They found a strong increase in the anisotropy of the normal Reynolds stresses (which decreased) in the outer region of the boundary layer, but observed the opposite e ect in the inner layer. By calculating the integral length-scales, they found that the near-wall vortices are stretched in the streamwise direction, but are only slightly smaller in the wall-normal direction. They measured high levels of atness of the instantaneous skin friction coe cient, which indicated the presence of intermittent high-frequency bursts in the re-laminarizing ow. Other recent experiments by Escudier et al.[98] found that the streamwise RMS velocity in the inner layer scales with the local freestream velocity, while in the outer layer it is \frozen" (i.e., remains relatively constant) in the accelerated region. They inferred that the frequency content of the turbulence was only changed at the highest frequencies, and that the bulk of the turbulence generated in the thick upstream boundary layer prior to acceleration was at frequencies too low to cause signi cant Reynolds stresses within the accelerated boundary-layer. Numerical calculations of a boundary layer with variable acceleration param- eters (as opposed to the sink ow studied by Spalart [34]) were carried out by Piomelli et al. [99], who examined the e ect of the acceleration on the near-wall vortical structures. They observed that the near-wall streaks became more elon- 96 gated and showed fewer undulations. Since they found that the vorticity levels in the acceleration region were similar to those in the zero-pressure gradient (ZPG) boundary-layer, and that the vortex scales in the cross-plane were unchanged, they suggested that the additional vortex-stretching due to the streamwise velocity gra- dient must be counterbalanced by other mechanisms. However, from the results presented, it was unclear which physical phenomena provide this balance. The mechanisms involved in the relaminarization of the turbulent boundary layer in an FPG are still poorly understood. What is clear, however, is that the turbulence in the outer layer remains frozen through the acceleration, and is, there- fore, strongly dependent on the conditions of the upstream boundary layer (i.e., neither on the local near-wall behavior nor on the local freestream velocity). Closer to the wall, the ow undergoes a process of laminarization, in which the skin fric- tion coe cient drops sharply. Finally, after the acceleration is completed, the ow quickly re-transitions to an equilibrium boundary layer. Several questions are still open; among them are: (1) Does the outer layer turbulence play a part in the re- laminarization phenomenon? (2) How do the inner and outer layers interact during and after the acceleration? (3) How does the re-transition to turbulence takes place (and why does it take place so abruptly)? In an attempt to clarify at least the rst of these points, large-eddy simula- tions (LES) of boundary layers in FPG are performed with di erent acceleration parameters. The rst simulation, which matches the high-acceleration experiment 97 of Warnack and Fernholz [97] (Kmax 4 10 6) (where K = ( =U21)(dU1=dx) is the acceleration parameter), shows a substantial reduction in turbulent kinetic energy production, and the ow becomes laminar-like in the acceleration region. A second simulation is performed at a lower acceleration (Kmax 3 10 6) (obtained by scaling down the high-K case) resulting in less of a \laminar-like" behavior in the acceleration (for comparison, on a 0.3m-long NACA0012 airfoil at Re = 1:5 106, K has values of the order 3 10 6 in the region between 2% of the chord and the point of maximum thickness). In the following we will present the numerical approach employed. Then will show results of the simulations, including statistics and ow visualizations; we will describe the results of some numerical experiments in which the ow was arti cially altered, to isolate the structures responsible for the re-transition. Finally, we shall draw some conclusions. 5.2 Problem formulation In this study we perform large-eddy simulations (LES) of a at-plate boundary layer in the presence of an accelerating freestream. Most experimental measurements in FPG boundary layer are performed on a at surface in which the pressure gra- dient is imposed through contouring of the opposite wall of the wind-tunnel, or by including a contoured body above the at wall to produce the desired acceleration 98 FPG ZPG APG Figure 5.1: Sketch of the con guration. The computational domain is shown as a hatched area. (Figure 5.1). In our case we imposed directly a variable freestream velocity1 U1(x) on the top boundary of the domain [102]. The other two velocity components were obtained by requiring that the vorticity in the freestream is zero. An unsteady in ow boundary condition was obtained from a separate simulation that used the recycling/rescaling method [102], while a convective out ow boundary condition was used at the downstream boundary.[82] Periodic conditions were used in the spanwise direction. The simulations were performed on a domain of size 476 o 20 o 20 o (in the streamwise, wall-normal and spanwise directions, respectively), using 1136 104 192 grid points for the high-acceleration case, and 1024 64 128 grid points for the low-acceleration one. The results obtained using this resolution compare well with coarser calculations. The Reynolds number, based on freestream velocity at the in ow, Uo, and on the displacement thickness at the in ow, o, is 1,260. 1In the following, U1 denotes the (variable) freestream velocity, whereas Uo = U1(0) = 1 is the reference velocity for the at-plate region upstream of the FPG region. 99 The equations of motion were integrated for 3171 o=Uo time units. Statistical data were obtained by averaging over the last 2415 time units, and over the spanwise direction. In the following, time-averaged quantities are denoted by angle brackets, and uctuating ones by a prime. 5.3 Results In the present investigation we studied two cases: one with high-acceleration, another with lower acceleration. Figure 5.2(a) shows the freestream velocity, U1=Uo, and the resulting friction velocity u =u ;o for the two cases of high and low acceler- ation. Figure 5.2(b) shows the acceleration parameter, K. In the case of high K, the freestream velocity at the out ow is almost three times that at the in ow. The momentum-thickness Reynolds number Re = U1 = , with = Z 1 0 U U1 1 UU 1 dy (5.1) (where U = hui) is shown in Figure 5.2(c); Re decreases as the ow begins to ac- celerate (x= o > 140). This re ects profound changes in the velocity pro le, which result in signi cant decrease of and will be discussed further later. The skin-friction coe cient is shown in Figure 5.2(d). Again, it can be observed that, although the mean freestream velocity increases in both cases, Cf begins to decrease near the loca- tion of maximum K (this decrease begins to occur downstream of the corresponding decrease in Re ). After the pressure gradient is relaxed, rapid re-transition towards an equilibrium turbulent value occurs (x= o ? 310). The chain-dotted line shows the 100 equilibrium ZPG value of Cf obtained from the correlation [104] Cf = 0:0576Re 0:2x . At the in ow, while the computations show good agreement with the experimen- tal correlation, the experiments have lower skin friction, suggesting that pressure gradient e ects are already signi cant on the upstream boundary layer. The region of relaminarization is predicted well by the LES. The re-transitioning occurs more abruptly than in the experiment, and the downstream Cf is higher in the LES, even taking into account the initial shift. Several factors can account for this di erence. First, the di erent geometry: in the experiment the measurements are made on the inside wall of a cylinder, whereas the calculation uses a at plate [94]; the ratio between cylinder radius and boundary-layer thickness varies between 12 and 9, so curvature e ects may play a role. Furthermore, despite the fact that the grid in the streamwise direction was re ned in the re-transition region, the three-fold increase in the friction velocity results in marginal resolution of the boundary layer in the fully turbulent region downstream of the acceleration: in the upstream region we have x+ ? 28, z+ ? 6:2, while in the recovery region x+ ? 56 and z+ ? 19. The lower acceleration case shows similar behavior to the high-acceleration one; however the re-laminarization is less severe. Re and Cf are reduced by a much smaller amount in the acceleration region, and recovery takes place earlier than for the high-K case. Figure 5.3 shows the mean velocity pro les in outer coordinates at several loca- tions in the ow. If the velocity pro les are plotted in wall coordinates (Figure 5.4), 101 0 100 200 300 400 5001 1.5 2 2.5 3 3.5 x / d*0 U ? /U o ; u t / u t o 0 100 200 300 400 5000 1 2 3 4 x 10?6 x / d*0 K 0 100 200 300 400 5002 3 4 5 6 x 10 ?3 x / d*0 C f 0 100 200 300 400 5000 500 1000 1500 2000 x / d*0 Re q U? / Uo ut / uto (a) (b) (c) (d) Figure 5.2: Streamwise development of statistical quantities. (a) Freestream ve- locity, U1=Uo and friction velocity u =u ;o; (b) Acceleration parameter K; (c) Momentum-thickness Reynolds number Re ; (d) Skin friction coe cient Cf. Ex- periments [97]; high-K case; low-K case; ZPG boundary-layer correlation. one can observe the existence of a logarithmic layer (following the standard law, U+ = 2:5 log y+ + 5) at the in ow and in the mild acceleration region (x= o < 150). As the FPG becomes signi cant, the slope of the logarithmic region decreases (a well-known e ect of acceleration [34]). The two cases are in good agreement up to the point of maximum K. Thereafter, the high-K case departs signi cantly from the equilibrium boundary layer pro le, becoming more laminar-like. The recovery 102 10?3 10?2 10?1 100 101 102 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 y/q U/U ? 0 50 100 150 200 250 300 350 400 4500 2 4 x/do* K ? 10 3 Figure 5.3: Wall-normal pro les of the mean streamwise velocity at the locations shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 0:5 units for clarity 103 10?1 100 101 102 103 104 0 15 30 45 60 75 90 105 120 135 150 U + y+ 0 50 100 150 200 250 300 350 4000 5 K x 10 6 Figure 5.4: Wall-normal pro les of the mean streamwise velocity in inner coordinates at the locations shown in the top gure. Experiments [97]; high-K case; low-K case; U+ = 2:5 log y+ + 5.Each pro le is shifted by 18:0 units for clarity 104 of the inner layer to an equilibrium logarithmic law occurs quite rapidly, between x= o ? 330 and 370. The agreement with the experimental data is very good. One should observe that in the region of high acceleration there is a signi cant region of well-mixed uid, in which the normal velocity gradient is nearly zero (from y= > 5 at x= o = 320, for instance). Very good agreement is also observed in the prediction of the normal Reynolds stresses (Figure 5.5). The decrease of the magnitude of the stresses following the maximum of the acceleration is evident, and is due to the increase of u2 , which is used to normalize the stresses, rather than to a decrease of the stresses themselves, which remain approximately equal to their upstream value (see the discussion below). Figure 5.6 shows the Reynolds stresses in inner units; again, there is very good agreement with the experimental results. The Reynolds stresses decrease sig- ni cantly in the region where K has the maximum value. In many cases ow re- laminarization is due to decorrelation between the wall-normal and the streamwise uctuation, which remain signi cant but do not contribute to the Reynolds shear stress. In this ow, the cause of the decrease of hu0v0i appears to be di erent. Figure 5.7 shows contours of the normal Reynolds stresses hu0u0i and hv0v0i, of the shear stresses hu0v0i, and of the correlation coe cient Cuv = hu 0v0i (hu0u0ihv0v0i)1=2 = hu0v0i urmsvrms : (5.2) From this gure it appear that, while the correlation coe cient decreases signi - cantly around x= o ? 300, the decrease of the v0 uctuations is much more dramatic. 105 10?1 100 101 102 103 104 0 10 20 30 40 50 60 70 80 90 100 y+ + 0 50 100 150 200 250 300 350 400 4500 5 K x 10 6 Figure 5.5: Wall-normal pro les of the streamwise Reynolds stresses at the locations shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 12:0 units for clarity 106 10?1 100 101 102 103 104 0 1 2 3 4 5 6 7 < u?v? > + 0 50 100 150 200 250 300 350 400 4500 5 y+ K x 10 6 Figure 5.6: Wall-normal pro les of the Reynolds shear stresses at the locations shown in the top gure. Experiments [97]; high-K case; low-K case. Each pro le is shifted by 0:8 units for clarity 107 y / d o* 0 2 4 6 8 0.02 0.01 0 (a) y / d o* 0 2 4 6 8 0.002 0.001 0 (b) x/do* y / d o* 0 100 200 300 4000 2 4 6 8 0 -0.001 -0.002 (c) x/do* y / d o* 0 100 200 300 4000 2 4 6 8 0.4 0.2 0 (d) x/d*o K ? 1 0 6 0 100 200 300 400.0 1.0 2.0 3.0 4.0 Figure 5.7: Contours of averaged quantities for the high-K case. The red line shows the local boundary layer thickness 99. (a) hu0u0i=U21; (b) hv0v0i=U21; (c) hu0v0i=U21; (d) Cuv = hu0v0i=urmsvrms. The same phenomenon is better illustrated in Figure 5.8. Here, we identi ed four streamlines (one in the outer layer, two in the logarithmic region and one in the bu er layer), and plot the streamwise development of the Reynolds stresses along each streamline. This method allows one to account for the signi cant thinning of the boundary layer in the high-acceleration region, and for its subsequent thickening in the re-transition region. Focusing our attention rst on the streamwise stresses, hu0u0i=U2o (Figure 5.8(b)), we observe that in the outer layer their level is essentially frozen (i.e., they remain constant and equal to their upstream value) until a location well after the maxi- mum acceleration point (x= o ? 300). They then begin to increase, an increase that occurs later for the streamline located farther away from the wall. Along the 108 0 50 100 150 200 250 300 350 400 4500 2 4 6 8 y / d o* x/d*o (a) 0 50 100 150 200 250 300 350 400 45010 ?1 100 101 10 3 ? /U o2 x/d*o 0 50 100 150 200 250 300 350 400 450 10?2 10?1 100 101 10 3 ? /U o2 0 50 100 150 200 250 300 350 400 45010 0 101 102 10 3 ? ?/U o2 (b) (c) (d) Figure 5.8: Development of the Reynolds stresses along selected streamlines. (a) Streamline coordinates; (b) hu0u0i=U2o ; (c) hv0v0i=U2o ; (d) hu0v0i=U2o . The thick black line corresponds to the boundary-layer edge. Streamline originating in: outer-layer; middle boundary layer; log region; viscous sublayer. 109 streamline in the bu er layer the increase in hu0u0i=U2o begins earlier, as soon as the freestream velocity begins to increase (x= o ? 120). A di erent behavior can be ob- served for the wall-normal stresses hv0v0i=U2o (Figure 5.8(c)): along the outer-layer streamlines the stresses also appear to be frozen to their upstream value. They begin to increase well after the peak acceleration, and also after the rise of the streamwise ones. Near the wall, on the other hand, we observe a signi cant decrease (by over one order of magnitude) of the wall-normal Reynolds stresses (consistent with the observations of Blackwelder and Kovasznay [93]) which is re ected in a similar de- crease of the Reynolds shear stress hu0v0i=U2o (Figure 5.8(d)) in the bu er layer and in the viscous sublayer. In the logarithmic region, on the other hand, the increase in the streamwise uctuation level balances the decrease of the wall-normal ones, resulting in constant shear stress until the recovery region (we will show later that the correlation coe cient does not vary very much in this region). Note that the data showed in this gure were normalized using the in ow freestream velocity, Uo, to emphasize advection e ects. If one used either the local freestream velocity, U1, or the friction velocity u , which increase through the acceleration region, all these quantities would be observed to decrease through the region in which the pressure gradient is applied. Figure 5.9 shows the structure parameter a1 = jhu0v0ij=hu0iu0ii and the corre- lation coe cient Cuv. The structure parameter is a measure of the e ciency of turbulence in extracting Reynolds shear stress from the available turbulent kinetic 110 0 100 200 300 400 5000 0.2 0.4 0.6 0.8 C uv x/d*0 0 50 100 150 200 250 300 350 400 4500 0.1 0.2 a 1 0 100 200 300 400 5000 2 4 6 8 y / d * 0 (b) (a) (c) Figure 5.9: Development of a1 parameter and correlation coe cient along selected streamlines. (a) Streamline coordinates; (b) a1; (c) Cuv. The thick black line corresponds to the boundary-layer edge. Streamline originating in: outer- layer; middle of the boundary layer; logarithmic region; viscous sublayer. 111 energy and, in equilibrium ows, it has a value close to 0.15. The acceleration causes a signi cant departure from its equilibrium value: we observe a signi cant decrease of a1 over the lower half of the boundary layer. The decrease is particularly strong in the logarithmic layer (the value of a1 is lower in this region even in equilibrium ZPG boundary layers). The correlation coe cient Cuv, on the other hand, remains close to the canonical value for at-plate boundary layers (Cuv ? 0:4 0:5) everywhere, with variations of less than 10%. Thus, the relaminarization does not seem to be due to so much to a decorrelation between the frozen uctuations, but rather to a re-organization of the ow that results in much lower wall-normal uctuations that, although reasonably well-correlated with the streamwise ones, can only produce a much reduced shear stress. The decreased mixing due to the turbulent transport, in turn, causes the decrease of the skin-friction coe cient that is considered one of the symptoms of relaminarization. The signi cant changes to the turbulent statistics observed above must be accompanied by similar alterations of the turbulent structure, which we will now describe. The contours of streamwise velocity uctuations in an xz plane near the wall are shown in Figure 5.10. We can observe the regular streaky structure of the boundary layer near the in ow. In the region of high K we observe uctuations of magnitude comparable to those in the equilibrium region; they form, however, very long streamwise streaks, without the kinks characteristic of the burst event. This indicates the change towards a very stable and more orderly inner layer, as pointed 112 Figure 5.10: Instantaneous contours of u0 velocity uctuations in a plane parallel to the wall. out by Narasimha and Sreenivasan [84]. Figure 5.10 also shows a fast re-transition to turbulence when the pressure gradient is turned o (x= o ? 310). In Figure 5.11 the coherent structures in the outer layer are visualized though isosurfaces of the second invariant of the velocity gradient tensor Q = 12 @u0 i @xj @u0j @xi (5.3) (see Dubief and Delcayre[105]). We note that the outer layer vortices become aligned in the streamwise direction in the acceleration region. This is most likely a kinematic e ect, as the dominant component of the velocity gradient in this region, @U=@x, has the e ect of stretching and re-orienting the coherent eddies into the streamwise direction. We note, however, that the more orderly structure of the ow observed in the inner layer (Figure 5.10) is re ected in a more orderly outer layer structure. 113 x/d * o 180 200 220 240 260 280 300 320 340 360 y / d * o 0 10 20 z/d * o 0 10 20 30 40 0.01 0.005 0 -0.005 -0.01 Outer layer Q=0.0002 Figure 5.11: Instantaneous isosurfaces of Q = 0:002 in the outer layer colored by the streamwise vorticity. 114 An obvious question that needs to be addressed is whether the inner layer, with its reduced burst frequency, is unable to \scramble" the outer layer structures, or if the outer layer forces a more orderly inner-layer structure. The outer layer structures certainly a ect the inner layer: Figure 5.12 shows contours of the instantaneous u0v0 correlation in planes normal to the mean ow, and secondary (v w) velocity vectors. One can observe that these vortices occur mostly in the well-mixed region mentioned above (the two thick lines that denote contours of u = 0:95U1 and 0:99U1 show the extent of this well-mixed region); they may, in fact, be responsible for it, as they induce vigorous motions of high-speed uid towards the wall, and low-speed uid outwards. Moreover, some of the motions induced by these large vortices result in signi cant values of the u0v0 correlation (at x= o = 321 and z= o = 11, or at x= o = 260 and z= o = 4, for instance). 5.4 Conclusions We performed LES of the ow in an accelerating boundary layer, using two dif- ferent levels of acceleration. The low-acceleration case remains in quasi-equilibrium, with a logarithmic law observed through most of the ow (albeit with decreased slope). The high-acceleration case results in relaminarization and re-transition of the ow. The computed statistics are in good agreement with the experimental data [96], which gives us con dence that the LES can be used to study the physics of this complex ow. 115 z/d*o y / d * o 0 5 10 150 2 4 6 8 -0.002 0 0.002 x/d*o=190 z/d*o y / d * o 0 2 4 6 80 1 2 3 4 -0.002 0 0.002 x/d*o=260 z/d*o y / d * o 0 5 10 150 2 4 6 8 -0.002 0 0.002 x/d*o=321 z/d*o y / d * o 0 5 10 150 2 4 6 8 -0.002 0 0.002 x/d*o=380 Figure 5.12: Contours of instantaneous u0v0 correlation in planes normal to the mean ow, and secondary (v w) velocity vectors. The two solid lines represent contours of U=U1 = 0:95 and 0.99. Examining the ow development along streamlines we observe that in the outer layer the turbulent uctuations appear to be largely frozen to their initial state, and the ow is dominated by advection. A notable feature of the ow is that the correlation coe cient Cuv does not decrease very signi cantly. The decrease of the Reynolds shear stresses that is observed is mostly due to the damping of the wall-normal uctuations. We observed changes in the turbulent structures both in the inner and in the outer layers. The acceleration a ects the outer-layer eddies by changing their structure and shape; in particular, large coherent structures are formed that are oriented in the streamwise direction. This results in the formation of a well-mixed outer layer, in which the turbulence production is decreased, and the turbulence 116 advected from upstream remains frozen. The inner layer is also a ected: because of the strong acceleration, the ow becomes more orderly, with longer, more two- dimensional streaky structures and decreased frequency of bursts. However, a fast re-transition to turbulence is observed as soon as the applied pressure gradient is negligible. This may be due to the coherent structures in the outer part of the ow that trigger the re-transition to turbulence. Two possible scenarios can explain the ow behavior: one, which matches the results of Blackwelder [93], Launder [89, 90], Narasimha and Sreenivasan [84], and Sreenivasan[83], is that the inner layer is made stable by the pressure gradient, and the turbulence in the outer layer remains relatively high and \frozen": once the stabilizing in uence of the pressure gradient is removed, transition occurs very rapidly, following a process that resembles bypass transition due to high freestream turbulence. A di erent picture was conjectured by Falco [94] and later by Ichimiya et al. [95]: the relaminarization seems to begin from the outer region with a strong coupling between inner and outer parts. The outer layer structures could induce strong incursions of more quiescent, outer-layer uid into the wall region, and strong ejections of inner-layer uid into the outer ow. Our data show that both of these mechanisms are present. Additional work is required to determine which of them is dominating. 117 Chapter 6 Conclusion and future work 6.1 Hybrid RANS/LES This work shows that a PI feedback control applied to the Navier-Stokes equa- tions accelerates the adjustment of synthetic turbulence generated at the in ow plane towards well-resolved LES turbulence. It has been proved that the PI con- trol is stable and robust with di erent non-equilibrium ows. In this Chapter we summarize the major ndings. 6.1.1 Control parameters The averaging window used to determine the error, Tave, has a strong e ect on the control, and it seems be the most critical parameter: if the averaging time is too large, the control acts on quantities that carry information not related with the actual state of the ow. On the other hand, if the averaging time is too short the input of the PI control does not have smooth stationary information, so it can input excessive energy into the system. If this parameter is not well tuned, two main e ects can be expected: the transient time of the calculation increases and the simulation requires large amounts of CPU time; also, the steady state error can be 118 large and a ect the results in the region downstream the last control plane. Another issue that has an in uence on the steady state results is the spanwise average that was used to de ne the error by previous researchers. The error, if no spanwise average is performed, re ects more accurately the local state of the ow, so that the PI control can locally adjust the control signal. Both the tuning of the time averaging window and the absence of spanwise average re ect the fact that in the feedback control is important to feed the PI control with the correct instantaneous information that allow for a correct response. Any delay can be dangerous. To summarize the e ect of KI and KP in the control it is conceptually useful to divide the computational domain in two regions: the controlled one, where the PI control is activated and the region downstream the control. The controlled region does not have a uid-dynamical interest because it is not a region where realistic turbulence exists, but it is seen as a transition from RANS to LES, in which a mismatch from the desired quantities is tolerated. However, from the point of view of the control, in this region the PI control tries to match the target: the two constants KI and KP determinate the steady-state match of the desired quantities. If these two constants are not well-tuned, we can have a large steady state error within that region. A residual error, however, has a weak e ect on the result in the region downstream the controlled one, as can be seen in the work of Keating et al. [30]. In fact, downstream of the last control plane the ow is biased to the exact solution of the NS equation. This is true especially when an external forcing is 119 applied to the ow (3D Boundary layer, Adverse Pressure gradient boundary layer): the bias in these cases is even stronger than in a boundary layer where no external force is applied (as the zero pressure gradient case). However, if the external forcing has the opposite trend of the control, the natural bias could not work. This is the case of the favorable pressure gradient: the acceleration damps the turbulent kinetic energy in the inner layer, instead the control tries to add energy towards a RANS solution with strong modeling errors. As was shown in the results, this kind of ow was extremely tedious to control. The same is true for the ZPG boundary layer: there is no external force and the solution in the region downstream the control is more sensible to the choice of the parameters. The coe cient KI and KP have a strong in uence in the stability of the system. Destabilization can be obtained if the PI system (or any other systems close in a feedback loop) introduce a lag in the loop. If KI is too large or KP too small the output of the PI control can be a ected by a time lag that can destabilize the Naveir-Stokes equations. However, for stable values of KP and KI, the results outside the control region are insensitive to these parameters. 6.1.2 Applications In the hybrid RANS/LES results, we found that an important part is played by the modeling error in the RANS solution. In the zero pressure gradient and adverse pressure gradient boundary layer, we found that a rapid adjustment to the 120 LES results is obtained when the k model is used in the RANS region. The fact that this model gives independent values of k and , while with the SA model must be estimated empirically, is a particularity desirable feature. In the favorable pressure gradient, on the other hand, any of the RANS model was unable to predict correctly the behavior of the ow. In the 3D boundary layer, we have seen that a recovery of the main shear stress hu0v0i and the secondary shear stress hv0w0i is obtained, with the forcing applied only to the hu0v0i stress. However, if the wrong time-averaging is used, the secondary stress is not matched. In the case of the zero pressure gradient and 3D boundary layer, we tried to force the turbulent kinetic energy production in order to have a control that was coordinate-invariant. We tested a production-based error in these two ows obtaining the same result as in the Reynolds stress-based error; the only di erence was in the controlled region: here the steady state error was large in the outer part of the boundary layer. 6.2 Resolved LES of the favorable pressure gradient boundary layer Large-eddy simulation of at-plate boundary layers in favorable pressure gra- dients (FPG) are performed for two di erent acceleration parameters. The high- acceleration case is in good agreement with the experimental data by Fernholz and Warnack [96, 97]. Substantial reduction in turbulent kinetic energy and shear stress production, and a reduction of the bursting frequency indicate that the inner part of 121 the accelerated boundary layer is in a laminar-like state when the pressure-gradient parameter K exceeds a threshold value. Downstream of this region, the boundary layer has a fast re-transition to turbulence. In the low K case, the boundary layer does not depart signi cantly from equilibrium. In the outer part of the boundary layer, the turbulent kinetic energy and the shear stress correlation maintained the same level as before the acceleration began. Two main hypotheses have been made: the inner and outer layers do not interact; or, the relaminarization is the conse- quence of inner-outer layer coupling. The results of our well-resolved calculations of the FPG boundary layer highlighted some signi cant di erences in the inner-outer layer interactions compared to a canonical ZPG case. They did not, however, allow us to determine whether the relaminarization process begins in the inner layer and propagates outwards, or if the reorganization of the outer layer forces the inner layer to become more stable. 6.3 Future work An important extension of the hybrid RANS/LES calculation is to compress- ible ow. The use of controlled forcing in this case can be challenging, because the control that we examined adds energy into the domain. One possible consequence is that not all the energy goes in turbulent uctuations but a fraction of that can be found as internal energy; this would make the method less e ective and could require di erent control strategies. One rst step to allow the extension to com- 122 pressible ow, could be the nding of more accurate control coe cients. This can be possible if a simpli ed model of the Navier-Stokes equation is found, so that well-known tested models can be applied to the feedback control. Regarding the favorable pressure gradient boundary layer, in order to give an answer to the issue of inner-outer layer coupling, a possible strategy is perform falsifying experiments [106]. In a rst set of simulations the ow can be manipulated to suppress the outer layer eddies and observes the behavior of the inner layer. In a second set of calculations the turbulence in the inner layer can be removed, and study the development of the outer layer in isolation. Another possibility is to study a single horseshoe vortex subjected to acceleration to determine whether the outer- layer re-organization is due to the freestream acceleration only, or if it is a ected by the decreased frequency of the bursts also. 123 Bibliography [1] Spalart, P.R. 2000 \Strategies for turbulence modelling and simulations" Int. J. Heat Fluid Flow, 21, 252-263. [2] Chapman, D.R. Computational aerodynamics, development and outlook. AIAA J 17 1293-1313 [3] Piomelli, U. and Balaras, E. 2002 \Wall-layer models for large-eddy simula- tions," Annu. Rev. Fluid Mech. 34, 349 (2002). [4] Sagaut,P. Deck,S. Terracol,M. 2006 Multiscale and multiresolution approaches in turbulence Imperial College Press . [5] Spalart, P.R. Jou, W.H. Strelets, M. and Allmaras, S.R. 1997 \Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach," In Advances in DNS/LES, edited by C. Liu and Z. Liu, (Greyden Press, Columbus, OH). [6] Balaras, E Benocci, C 1994 Subgrid-scale models in nite-di erence simulations of complex wall bounded ows. AGARD CP 551 pp 2.1-2.5 [7] Speziale,C. 1997 Turbulence modelling for time-dependent RANS and VLES:A review. AIAA Paper 97-2051 [8] Batten,P. Goldberg,U. and Chakravarthy,S. 2000 Sub Grid Turbulence Model- ing for Unsteady Flow with Acoustic Resonance AIAA Paper 2000-0473 . [9] Spalart, P.R., Deck, S., Shur, M.L., Squires, K.D., Strelets, M. Kh. and Travin, A., 2006, A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theo. Comput. Fluid Dyn., 20, 181{195. [10] Girimaji, S.S. Abdol-Hamid, K.S. 2005 Partially-averaged Navier Stokes Model for Turbulence: Implementation and Validation AIAA 2005-502 pp 1-14 [11] S. Deck, 2005 Zonal-detached-eddy simulation of the ow around a high-lift con guration AIAA J. 43:11,2372{2384. 124 [12] Nikitin, N.V., Nicoud, F., Wasistho, B., Squires, K.D. and Spalart, P.R. 2000 An approach to wall modeling in large-eddy simulations. Physics of Fluids, 12, 1629-1632. [13] Cabot, W. Moin, P. 2000 Approximate wall boundary conditions in the large- eddy simulation of high Reynolds number ows. Flow Turbulent Comb 63, 269 [14] Hamba, F. 2003 A hybrid RANS/LES simulation of turbulent channel ow. Theoret. Comput. Fluid Dyn., 16, 387-403. [15] Piomelli, U. Balaras, E. Pasinato, H. Squires, K.D. and Spalart,P.R. 2003 The inner-outer layer interface in large-eddy simulations with wall-layer models Int. J. Heat Fluid Flow, 24, 538. [16] Davidson, L. and Dahlstrom, S. 2005 Hybrid LES-RANS: An approach to make LES applicable at high Reynolds number. 19, 415-427. [17] Temmerman, L., Hadziabdic, M., Leschziner, M. A. and Hanjalic, K. 2005 A hybrid two-layer URANS-LES approach for large eddy simulation at high Reynolds numbers. International journal of heat and uid ow, 26, 173-190. [18] Davidson, L. and Billson, M. 2006 Hybrid LES-RANS using synthesized tur- bulent uctuations for forcing in the interface region. International journal of heat and uid ow, 27, 1028-1042 [19] Hamba, F. 2006 A hybrid RANS/LES simulation of high-Reynolds-number channel ow using additional ltering at the interface. Theoret. Comput. Fluid Dyn., 20, 89-101. [20] Keating, A. and Piomelli, U. 2006 A dynamic stochastic forcing method as a wall-layer model for large-eddy simulation. J. of Turbulence, 7(12), 1-24. [21] Tessicini, F. Temmerman,L. Leschziner, M.A. 2006 Approximate near-wall treatments based on zonal and hybrid RANS-LES methods for LES at high Reynolds numbers.Int. J. Heat and uid ow, 27, 789-799. [22] Benarafa, Y. Cioni, F. Ducros, F. Sagaut, P. 2006 RANS/LES coupling for unsteady turbulent ow simulation at high Reynolds number on coarse meshes Comput. Methods Appl. Mech. Engrg. 195, 2939-2960 125 [23] Dahlstrom S., Davidson L. 2003 Hybrid LES-RANS with additional conditions at the matching region Turbulence heat and mass transfer 4 Ney York, Walling- ford (UK) [24] Labourasse,E. and Sagaut,P. 2002 Reconstruction of turbulent uctuations us- ing a hybrid RANS/LES Approach Journal of Comp. Physics 182, 301-336 [25] Terracol, M. 2005 A zonal RANS/LES approach for noise sources prediction. In Engineering Turbulence Modelling and Experiments 6 edited by W. Rodi and M. Mulas, Elsevier, Amsterdam, 699{708. [26] Sagaut, P., Garnier, E., Tromeur, E., Larchev^eque, L., and Labourasse, E, 2004, Turbulent in ow conditions for LES of compressible wall-bounded ows, AIAA J. 42, 469{477. [27] Sandham, N.D., Yao, Y.F., Lawal, A.A., 2003, Large-eddy simulation of tran- sonic turbulent ow over a bump Int. J. Heat FLuid Flow, 24, 584{595. [28] Mathey, F. Cokljat, D. Bertoglio, J.P. and Sergent, E., 2006, Assessment of the vortex method for Large Eddy Simulation inlet conditions Progress in Comp Fluid Dynamics 6, 58{67 [29] Sanchez-Rocha, M. Kirtas, M. Menon, S.2006 Zonal Hybrid RANS-LES method for static and oscillating airfoils and wings AIAA Paper 2006-1256 [30] Keating, A., De Prisco, G., and Piomelli, U., 2006 Interface conditions for hybrid RANS/LES calculations Int. J. Heat FLuid Flow, 27, 777{788. [31] Batten P., Goldberg, U. and Chakravarthy, S.2004 Interfacing statistical tur- bulence closures with large-eddy simulation, AIAA J. 42, 485{492. [32] Davidson,L. 2005 Hybrid LES-RANS: Inlet boundary conditions Proc. 3nd Na- tional Conf. Comput. Mech., 7{22. [33] Sandham, N.D., Yao, Y.F., Lawal,A.A. 2003 Large-eddy simulation of transonic turbulent ow over a bump Int. J. Heat FLuid Flow, 24, 584{595. [34] Spalart,P.R. 1986 Numerical study of sink- ow boundary layers J. Fluid Mech., 172, 307 (1986). 126 [35] Spalart,P.R. 1986 Direct simulation of a turbulent boundary layer up to Re =1410 J. Fluid Mech., 187, 61. [36] Spalart, P.R. and Watmu ,J.H. 1993 Experimental and numerical study of a turbulent boundary layer with pressure gradients J. Fluid Mech., 249, 337. [37] Lund,T.S. Wu,X. and Squires,K.D. 1998 Generation of in ow data for spatially- developing boundary layer simulationsJ. Comput. Phys., 140, 233. [38] Ferrante A. and Elghobashi,S.E. 2004 A robust method for generating in ow conditions for direct simulations of spatially-developing turbulent boundary layers J. Comput. Phy. 198, 372. [39] Schl uter, J.U., Wu, X., Kim, S. Alonso, J.J. and Pitsch, H., 2004, Coupled RANS-LES computation of a compressor and combustor in a gas turbine engine. AIAA Paper 2004-3417. [40] Schl uter,J.U. Pitsch,H. and Moin,P. 2004 LES in ow conditions for coupling with Reynolds-averaged ow solvers AIAA J., 42, 478. [41] Medic,G. You,D. and Kalitzin,G. 2006 An approach for coupling RANS and LES in integrated computations of jet engines CTR Annual Research Briefs [42] Schluter J.U., Pitsch H., Moin P. 2002 Consistent boundary conditions for integrated LES/RANS conditions AIAA paper [43] Batten P., Goldberg, U. and Chakravarthy, S.2002 LNS- An approach Towards Embedded LES, AIAA 2002-0427, 1{10. [44] Le,H. Moin,P. and Kim,J. 1997 Direct numerical simulation of turbulent ow over a backward-facing step J. Fluid Mech., 330, 349. [45] Lee,S. Lele,S. and Moin,P. 1992 Simulation of spatially evolving compressible turbulence and application of Taylors hypothesis," Phys. Fluids A 4, 1521. [46] Kraichnan,R.H. 1969 Di usion by a random velocity eld Phys. Fluids 13, 22-31 [47] Smirnov,A. Shi,S. and Celik,I.2001 Random ow generation technique for large eddy simulations and particle-dynamics modeling J. Fluids Eng., 123, 359. 127 [48] Klein,M. Sadkiki,A. and Janicka,J. 2003 A digital lter based generation of in ow data for spatially developing direct numerical or large eddy simulations J. Comput. Phys., 186, 652. [49] Spille-Koho ,A. and Kaltenbach,H.-J. 2001 Generation of Turbulent In ow Data with a Prescribed Shear-Stress Pro le. In DNS/LES Progress and chal- lenges, edited by C. Liu, L. Sakell, T. Beutner (Greyden, Columbus, OH, 2001), 319{326. [50] Keating, A., Piomelli, U., Balaras, E. and Kaltenbach, H.-J., 2004, A priori and a posteriori tests of in ow conditions for large-eddy simulation, Phys. Fluids 16, 4696{4712. [51] Keating, A. 2003 Large-eddy simulation of Heat transfer in turbulent channel ow and in the turbulent ow downstream of a backward-facing step, PhD dissertation, The University of Queensland, Brisbone, Australia [52] Lilly, D.K. 1967, The representation of small-scale turbulence in numerical sim- ulation experiments Proceedings of the IBM Scienti c Computing Symposium on Environmental Sciences, Yorktown Heights, NY [53] Moin, P. and Kim, J. 1982 Numerical Investigation of Turbulent Channel Flow, Journal of Fluid Mechanics 118, 341-377 [54] Piomelli, U. and Zang, T.A. and Speziale, C.G. and Hussaini, M.Y. 1990, On the large-eddy simulation of transitional wall bounded ows, Physics of Fluids A 1 257-265 [55] Clark, R.A. and Ferziger, J.H. and Reynolds, W.C. 1979 Evaluation of subgrid- scale turbulence models using an accurately simulated turbulent ow Journal of Fluid Mechanics 91, 1-16 [56] Liu, S. and Menveau, C. and Katz, J. 1994, On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet Journal of Fluid Mechanics 275, 83-119 [57] Vreman, B. and Geurts, B. and Kuerten, H. 1997 Large-eddy simulation of the turbulent mixing layer Journal of Fluid Mechanics, 339, 357-390 [58] Bardina, J. and Ferziger, J.H. and Reynolds, W.C. 1980, Improved subgrid- scale models for large-eddy simulations, AIAA Paper 80-1357 128 [59] Germano,M. Piomelli,U. Moin,P. and Cabot,W. 1991 A dynamic subgrid-scale eddy viscosity model Phys. Physics of Fluids 3, 1760. [60] Lilly, D.K. 1992 A proposed modi cation of the Germano subgrid-scale closure model, Physics of Fluids A, 4, 633-635 [61] Spalart, P.R. and Moser, R.D. and Rogers, M.M. 1991 Spectral methods for the Navier-Stokes equations with one in nite and two periodic directions Journal of Computational Physics 96, 297-324 [62] Akselvoll, K. and Moin, P. 1996 An e cient method for temporal integration of the Navier-Stokes equation in con ned axisymmetric geometries Journal of Computational Physics 125 454-463 [63] Morinishi, Y. and Lund, T.S. and Vasilyev, O.V. and Moin, P. 1998 Fully con- servative higher order nite di erence schemes for incompressible ow Journal of Computational Physics 143, 90-124 [64] Moin, P. and Mahesh, K. 1998 Direct numerical simulation: A Tool for Turbu- lence Research, Annual Review of Fluid Mechanics 30, 539-578 [65] Harlow, F.H. and Welch, J.E. 1965 Numerical calculation of time-dependent viscous incompressible ow of uid with free surface. The Physics of Fluids 8, 2182-2189 [66] Akselvoll, K. 1995 Large eddy simulation of turbulent con ned coannular jets and turbulent ow over a backward facing step, PhD thesis, Stanford University [67] Piomelli, U. 1993 High Reynolds number calculations using the dynamic subgrid-scale stress model, Physics of Fluids A 5, 1484-1490 [68] Ghosal, S. and Moin, P. 1995 The basic equations for the large eddy simulation of turbulent ows in complex geometry, Journal of Computational Physics 118, 24-37 [69] Cabot, W. and Moin, P. 1993 Large eddy simulation of scalar transport with the dynamic subgrid-scale model, Large eddy simulation of complex engineering and geophysical ows, Cambridge University Press [70] Spalart, P.R. and Allmaras, 1992 S.R. A one-equation model for aerodynamic ow AIAA paper 92-0439 129 [71] Launder, B.E., and Spaulding, D.B., 1972, Lectures in mathematical models of turbulence, Academic Press, London. [72] Chen, H.C., and Patel, V.C., 1988, Near-wall turbulence models for complex ows including separation. AIAA J. 26, 641{648. [73] Bradshaw, P. Ferriss, D.H. Atwell, N.P. 1967 Calculation of boundary layer development using the turbulent energy equation Journal Fluid Mechanics 23, 31-64 [74] Townsend, A.A. Equilibrium layers and wall turbulence 1962 Journal Fluid Mechanics 11 97-120 [75] Menter, F.R. Eddy viscosity transport equations and their relation to k model Journal Fluids Engineering 119, 876-884 [76] Cohen, G.H. and Coon, G.A. 1953 Theoretical consideration of retarded control, Trans. ASME 75, pp 827-834 [77] Balaras,E. Benocci,C. Piomelli,U. Two-layer approximate boundary conditions for large-eddy simulations 1996 AIAA J., 34, 1111. [78] Benarafa Y., Cioni F., Ducros F., Sagaut P. 2006 RANS/LES coupling for unsteady turbulent ow simulation at high Reynolds number on coarse meshes" Comput. Methods Appl. Mech. Engrg, 195, 2939-2960 [79] Spalart, P.R. and Allmaras, S. R. 1994 La Recherche A erospatiale, 1, 5{. [80] Meneveau,C. Lund,T.S. and Cabot,W.H. 1996 A Lagrangian dynamic subgrid- scale model of turbulence J. Fluid Mech. 319, 353{. [81] Kim,J. and Moin,P. 1985 Application of a fractional step method to incom- pressible Navier-Stokes equations J. Comput. Phys., 59, 308. [82] Orlanski,I. 1976 A simple boundary condition for unbounded hyperbolic ows J. Comput. Phys. 21, 251. [83] Sreenivasan, K.R. (1982) Laminariscent, Relaminarizing and Retransitional Flows. Acta Mech. 44, pp.1{48. 130 [84] Narasimha,R. and Sreenivasan, K.R. (1973) Relaminarization in highly accel- erated turbulent boundary layers. J. Fluid Mech. 61, pp.417{447. [85] Narasimha,R. and Sreenivasan, K.R. (1979) Relaminarization of Fluid Flows. Adv. Appl. Mech. 19, pp.221{309. [86] Wilson, D. (1954) Convective heat transfer to gas turbine blades. Proc. Inst. Mech. Eng. [87] Patel, V. and Head, M. (1968) Reversion of turbulent to laminar ow. J. Fluid Mech. 34, pp 371-392 [88] Senoo, Y. (1972) The boundary layer on the end wall of a turbine nozzle cascade. Transaction of the ASME. MIT Gas Turbine Laboratory. 35, pp 1711-1720 [89] Launder, B.E., (1963) The Turbulent Boundary Layer in a Strongly Negative Pressure Gradient Gas. MIT Gas Turbine Laboratory Report 71, pp. 1{28 [90] Launder, B.E., (1964) Laminarisation of the Turbulent Boundary Layer by Acceleration. MIT Gas Turbine Laboratory Report 77, pp. 1{50 [91] Badri Narayanan, M., Rajagopalan, S. and Narasimha, R. (1977) Experiment on the ne structure of turbulence. J. Fluid Mech. 80, pp. 237{257. [92] Kline, S., Reynolds, W., Schraub, F. and Runstadler, P. (1967) The structure of turbulent boundary layers. J. Fluid Mech. 30, pp. 741{773. [93] Blackwelder, R.F. and Kovasznay, L.S.G. (1972) Large-scale motion of a tur- bulent boundary layer during relaminarization J. Fluid Mech. 53, pp. 61{83 [94] Falco, R.E. (1980) Combined simultaneous ow visualization/hot-wire anemometry for the study of turbulent ows J. Fluids Eng. 102, pp. 174{182 [95] Ichimiya, M., Nakamura, I., and Yamashita, S. (1998) Properties of a relam- inarizing turbulent boundary layer under a favorable pressure gradient Exp. Thermal Fluid Sci. 17, pp. 37{48 [96] Fernholz, H.H. and Warnack, D. (1998) The e ects of a favourable pressure gradient and of the Reynolds number on an incompressible axisymmetric turbu- lent boundary layer. Part 1. The turbulent boundary layer.J. Fluid Mech. 359, pp. 329{356. 131 [97] Warnack, D. and Fernholz, H.H. (1998) The e ects of a favourable pressure gradient and of the Reynolds number on an incompressible axisymmetric tur- bulent boundary layer. Part 2. The boundary layer with relaminarization.J. Fluid Mech. 359, pp. 357{381. [98] Escudier, M.P., Abdel-Hameed, A., Johnson, M.W. and Sutcli e, C.J. (1998) Laminarisation and re-transition of a turbulent boundary layer subjected to favourable pressure gradient. Exp. Fluids 25, pp. 491{502. [99] Piomelli, U., Balaras, E., & Pascarelli, A. (2000) Turbulent structures in accel- erating boundary layers. J. Turbulence 1 (001) pp. 1{16. [100] Germano, M., Piomelli, U., Moin, P. and Cabot W. H. (1991) A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3, pp. 1760{1765. [101] Lilly, D. K. (1992) A proposed modi cation of the Germano subgrid-scale closure method. Phys. Fluids A 4, pp. 633{635. [102] Lund, T. S., Wu, X., and Squires, K. D. (1998) Generation of turbulent in ow data for spatially-developing boundary layer simulations, J. Comput. Phys. 140, pp. 233{258. [103] Chorin, A. J. (1968) Numerical solution of the Navier-Stokes equations. Math. Comput. 22, pp. 745. [104] Schlichting, H. (1979) Boundary layer theory, McGraw Hill, New York. [105] Dubief, Y., and Delcayre, F. (2000) On coherent-vortex identi cation in tur- bulence. Journal of Turbulence 1, 1-22. [106] Jimenez, J. and Pinelli, A. 1999, The autonomous cycle of near-wall turbu- lence. J. Fluid Mech., 389, 335 { 359 132