Large Eddy Simulation with the unstructured collocated arrangement Sofiane Benhamadouche A thesis submitted to the University of Manchester for the degree of Doctor of Philosophy School of Mechanical, Aerospace and Civil Engineering January 20, 2006 ii Declaration No portion of the work referred to in this thesis has been submitted in support of an application for another degree or qualification of this or any other university, or other institution of learning. iii iv Acknowledgments I would like to express my sincere gratitude to Professor Dominique Lau- rence, my supervisor, who allowed me to carry out this Phd, in parallel to my research work at Electricit´e De France in Paris, with the financial support of an EPSRC grant. I take also this opportunity to thank my wife Rym for her patience during my travels to Manchester and for her support during the period of this thesis. I also gratefully acknowledge Dr. Y. Addad, I. Afgan, V. Boyer, S. Gant, V. Garibian, V. Guimet, N. Jarrin, Dr. C. Moulinec, A. Revell, C. Robinson, P. Shepherd, J. Uribe and others who I forgot for discussions, support, coffee breaks and other social events. v vi Abstract Several industrial applications of CFD, such as thermal fatigue, aero-acoustic noise or fluid-structure interaction need to capture the unsteady behaviour of an incompressible flow. Large Eddy Simulation is the only reasonable approach in CFD toobtainatime-dependentsolution,aswelldemonstratedonacademiccases. However,performingLES onfullyunstructuredfinite-volumegrids,asrequiredfor complexindustrialgeometries, aswellastolocallyadaptthegridresolutiontothe largespatialvariationsofturbulentscales, remainschallenging. Themainpurpose of the present work is to deal with the collocated arrangement implemented in the in-house E´lectricit´e De France (EDF) code Code Saturne and to show whether reasonable LES computations are feasible with this discretization. Several numerical issues such as numerical dissipation are first analyzed in a discrete sense. Itis foundthattheRhie andChowinterpolation, widelyusedinall the collocated approaches, introduces a numerical dissipation which drains energy equallyatallthewave-numbers. RemovingthisinterpolationonregularCartesian grids allows, with a second order Crank-Nicolson scheme and sub-iterations on the predictor-corrector algorithm for pressure/velocity coupling, to strictly conserve kinetic energy. Nevertheless, the interpolation is maintained as it has a stabilizing effect on skewed grids and as it avoids odd-even decoupling phenomenon. It is also foundthata’symmetrical’formulationallowstostrictlyconservekineticenergyfor the convection term. The use of this formulation alters the second order precision of the code. Thus, a second order approach for the convection term is kept with implicit gradient reconstruction for non-orthogonal grids. Finally, sub-iterations on the predictor-corrector algorithm for pressure/velocity coupling are retained as this allows to reduce time splitting errors. A particular attention is then paid to the Decaying Isotropic Turbulence in order to identify the implicit filter induced by the discretization and to find the appropriate Smagorinsky constant. It is found that the implicit filter induced by the numerical discretization with a Smagorinsky model resembles to a Gaussian filter with a filter-width equal to 2h where h is the grid-spacing in a regular mesh. The value of the Smagorinsky constant is between 0.16 and 0.18 what corresponds to consensual values found in the literature for highly accurate discretizations. The explicit filter employed for the dynamic model is a smoothing operator that involves all the neighbours of a cell which share a node with it. The filter-width of this explicit filter is equal, on uniform meshes, to 3h. The averaging of the Smagorinsky constant is done for its numerator and denominator separately. It is also shown that local averaging gives reasonable results with the DIT which is encouraging for industrial simulations in which no homogeneous direction is available. The last two academic cases to be tested are the standard and the oscillating channelflows. AfterusualtestsonstructuredCartesianmeshes, theabilityofnon- conformingmeshestoadapttheresolutionlargeeddyscalevariationsatthewallis pointed out. The global over-estimation of the mean velocity often observed with vii Cartesian meshes is cured by introducing local refinements in the near-wall region using the capabilities of a collocated approach to naturally handle non-conforming meshes. Thissolutionisverypromisingtorunindustrialcasesatmoderatetohigh Reynolds numbers without artificial wall functions. To conclude the channel flow simulations, a distributed non-conforming mesh is created following the behaviour oftheKolmogorovscalesineachdirectioninordertoobtainaveryfinemeshcom- patible with a Direct Numerical Simulation resolution but with a lower number of computationalcells thaninthe conformingcases. The results are verysatisfactory compared to previous DNS calculations on structured grids. This capability can be exploited to run DNS calculations at higher Reynolds numbers with much less nodes than in the structured case. Finally, The potential of LES to predict complex vortical separation is re- ported using a 3D axi-symmetric bump considered very challenging by the CFD community as all Reynolds Averaged Navier-Stokes attempts totally failed. The calculations give satisfactory quantitative results. The main physical structures observed in the experiment are captured. The quantitative results are globally satisfactory except for the wall-normal componentofthe velocity whichis farfrom the experimental results and for the Reynolds stresses which are overestimated. viii Nomenclature Lower-case Roman d The dimension of the problem i Component in a Cartesian coordinate system (i ∈ {1,2,3}) k Inner iteration k Residual turbulent kinetic energy r k Resolved turbulent kinetic energy or the turbulent kinetic energy in a e continuous sense l Integral length-scale 0 m Flux of the velocity through an inner face between the cells Ω and Ω IJ I J n External normal to a control volume (k n k= 1) n Normal vector to the common face to the cells Ω and Ω (goes from I to J, IJ I J k n k= 1) IJ n Time iteration p Pressure t Time variable tn Time at the nth iteration u Velocity magnitude or u velocity component 1 u ith velocity component i u Velocity component parallel to the wall τ u Velocity vector ∂u u∗ Friction velocity (ρu∗2 = µ | τ |w) ∂n v u velocity component 2 w u velocity component 3 x ith Cartesian coordinate i x Cartesian coordinate x, y Position vectors y Distance from the wall or Cartesian coordinate z Cartesian coordinate u∗y y+ Non dimensional distance to the wall (y+ = ) ν ix Upper-case Roman C Smagorinsky constant S C Square of the Smagorinsky constant (C = C2) D D S C Universal Kolmogorov constant (C = 0.15) Kol Kol Z 1 E Total kinetic energy ( u.udΩ) 2 Ω or the density of turbulent kinetic energy (k = R E(K)dK) e F Centre of gravity of inner or boundary faces G Kernel of a spatial filter (R G = 1) I Centre of gravity of the control volume Ω I I0 Projection of I on the normal to the face which crosses F I Identity tensor J The centre of gravity of the neighbour of a control volume Ω I K, K Wave-number and Wave-number vector K Cut-off wave-number c O Intersection of IJ with the face common to the control volumes Ω and Ω I J Re Reynolds number 1 S Rate of strain tensor (S = (∇u+t ∇u)) 2 q |S| or S Rate of strain invariant (|S| = S = 2S S ) ij ij S Surface of the common face to the control volumes Ω and Ω IJ I J dS Elementary surface T Temperature x