Models vs. Physical Laws, Vienna, 2nd Feb 2011 Single-point second-moment turbulence models – why, where and where not M A Leschziner Shock-induced separation The holy grail We are promised a ‘model-free’ CFD world A Boeing 747 is not a homogeneous square box! Hybrid LES-RANS Courtesy: ANSYS, Germany Some scales and estimates Aircraft: Re 108 Nη 1019 Time steps: Nτ 106 − 107 Nodes: Current estimate of time of realisation: 2080 Current estimate for LES: 2045 (based on resolution at Taylor scale) Current capability: RANS and RANS-LES hybrids 95%+ of all engineering CFD is based on RANS The cost Mesh: 1019 Cost: 5000 CPU years per 1 second of flying at 1Tflop throughput 50 CPU years 5000 CPU years GPUs? BPUs? DNS - Status Model-free DNS used to ¾ investigate fundamental physics; ¾ examine subgrid-scale models (a-priori testing) ¾ Develop, calibrate and validate RANS models Largest channel-flow DNS: Reτ = 964 , 2.7x109 nodes (Del Alamo et al, 2004) Example: insight into origin of drag reduction by spanwise wall oscillation (Touber & Leschziner, 2010) Travelling surface wave of spanwise motion ¾ Reτ = 500 (→ 1000) ¾ Drag reduction up to 40% kx ¾ 0.5x109 nodes, 1M CPU hours ω Streamwise drag Fundamental mechanism of streak response Streak formation and re-orientation mechanisms Conditional sampling and averaging Decomposition of small streaks/super-streaks Modulation mechanisms Reduction of wall-normal Linear analysis (GOP) fluctuations around streaks in % due to actuation Streak decay, regeneration, reorientation and modulation The “RANS” equations Time-averaged framework: ∂ρU i U j ∂ P ∂ ⎛ ∂U i ∂U j μ⎜ =− + + ∂x j ∂xi ∂x j ⎜⎝ ∂x j ∂xi ⎞ ∂ ρ ui u j + Bi ⎟⎟ − x ∂ j ⎠ Unsteady – URANS framework ¾ Triple decomposition % U= U { + u{ + u{' Mean Periodic 144244 3 Stochastic Phase-average U u% ∂ρU jU i ∂x j ( ∂ 2U i ∂P ∂ = +μ + ρ −u%i u% j − ui'u 'j ∂xi ∂x j ∂x j ∂x j “coherent” ) ∂ρ u%i ∂ρU j u%i ∂p% ∂ 2u%i ∂ + = +μ + ρ u%i u% j − u%i u% j ∂t ∂x j ∂xi ∂x j ∂x j ∂x j ( ( ) ) ∂ ∂ρU i u% j ' ' ' ' + ρ ui u j − < ui u j > − ∂x j ∂x j ¾ Requires closure equations for the periodic and stochastic terms: too complex in practice – URANS use RANS models + ∂ / ∂t Nature of Modelling Reynolds stresses related to known or determinable quantities: ⎛ Sij , Ωij , Skl Skl , Ω kl Ω kl , ⎞ ⎜ ⎟ ui u j = fij ⎜ length-scale surrogates ⎟ ⎜ turbulence invariates ⎟ ⎝ ⎠ Sij Strain tensor Ωij Vorticity tensor Ultimately, need to relate to stresses and mean velocity. Modelling principles – not only “ad-hoc curve fitting” ¾ strong fundamental foundation; ¾ resolution of anisotropy; ¾ correct response to shear and normal straining; ¾ correct response to curvature and body forces; ¾ frame-invariance (“objectivity”); ¾ realisability; ¾ correct approach to 2-component turbulence at wall and fluid-fluid interfaces; ¾ satisfactory numerical stability; ¾ economy. Model types – basic classification About 150 models & major variations, many meant for restricted flow classes −ui u j 2 Linear eddy viscosity −ui u j = 2ν t Sij − kδ ij models 3 Algebraic Second-moment closure (Reynoldsstress models) Differential (turbulence transport) −uiθ “Algebraic” secondmoment closure Non-linear eddy viscosity models 1-equation 2-equation + v2f Eddy-viscosity transport Turbulence-energy transport Turbulence energy + various length-scale surrogates Wall-normal energy component + various length-scale surrogates Defects of linear eddy-viscosity models Linear EVM: ¾ Well suited to thin shear flow ¾ Much less well suited to separated and highly 3d flow ¾ No resolution of anisotropy ¾ Wrong sensitivity to flow curvature, rotation, normal straining and body forces ¾ Reliant on ad-hoc corrections Defects are rooted in ¾ Inapplicability of linear stress-strain relations ¾ Isotropic nature of viscosity, relating to scalar turbulence properties ¾ Calibration by reference to simple, near-equilibrium flows ¾ Excessive extrapolation to complex condition. Only fundamentally credible alternative ¾ Modelling based on exact equations for the Reynolds stresses ¾ Strong resistance from engineering community - complexity Reynolds-Stress-Transport Modelling Introduce the Reynolds decomposition U i = U i + ui etc. into the NS equations. Subtract from this the corresponding RANS equation. Repeating the above, but with the indices i and j interchanged. Add the two equations. Time-averaging the result: ∂U j ⎧ ∂U i ⎫ ∂ui ∂u j = − ⎨ ui uk + u j uk + f u + f u 2 ν − ⎬ i j j i Dt ∂ ∂ ∂xk ∂xk 14 4 244 3 x x k k ⎭ ⎩14444 123 14 243 24444 3 F ij C ( Dui u j ij ) ε ij Pij pu j ∂ui u j ⎫⎪ pu i ∂ ⎧⎪ p ⎛ ∂ui ∂u j ⎞ + + ⎜ ⎟⎟ − δ ik + δ jk -ν ⎨ui u j uk + ⎬ ⎜ ∂xk ⎪⎩ ∂xk ⎪⎭ ρ ⎝ ∂x j ∂xi ⎠ ρ ρ 1442443 14444444 4244444444 3 Φ ij dij Pressure-velocity ¾ Cij , Pij , Fij , Φ ij , ε ij and dij represent, respectively, stress convection, production by strain, production by body forces (e.g. buoyancy ), dissipation, pressure-strain redistribution and diffusion The Argument for resolving anisotropy Production is a key process: it drives the stresses. It requires no approximations if stresses and velocity are known It is reasonable to assume, tentatively: Stress = Production x Time (capital = interest rate x time) Exact equations imply complex stress-strain linkage ∂U j ⎧ ∂U i ⎫ + u j uk ρ uiu j ←⎯ → − τ ⎨ uiuk ⎬ + τ × Body -force production ∂ ∂ x4 x k3⎭ k ⎩1444 24444 Pij Hence, simple EVM stress-strain linkage is inapplicable Analogous linkage between scalar fluxes and production ⎧ ∂Φ ∂U i ⎫ + uiϕ ⎬ + τ ϕ × Body -force production ∂ xk ∂ xk ⎭ ⎩1444 424444 3 → − τ ϕ ⎨ uiuk ρ uiϕ ←⎯ Pu ϕ i Hence, Fourier-Fick law (eddy-diffusivity approximation) ρ uiϕ = − not valid μt ∂Φ σ ϕ ∂xi The equations for thin shear flow Only one shear strain, only one shear stress Duv p ⎛ ∂u ∂v ⎞ ∂ ⎛ 2 pu ⎞ μ ∂uv 2 ∂U = −v + ⎜ + ⎟ − ⎜ uv + − ε12 ⎟+ ρ ⎠ ρ ∂y ∂y ρ ⎝ ∂y ∂x ⎠ ∂y ⎝ Dt ⎫ ⎪ ⎪ ⎪ μ ∂v2 p ∂v ∂ ⎛ 3 2 pv ⎞ ⎪ +2 − ⎜v + + − 2 ε ⎟ 22 ⎬ ρ ∂y ∂y ⎝ ρ ⎠ ρ ∂y ⎪ ⎪ 2 μ ∂w p ∂w ∂ ⎪ − ε 33 +2 − vw2 + 2 ⎪⎭ ρ ∂z ∂x ρ ∂y ( ) μ ∂u 2 ∂U Du 2 p ∂u ∂ 2 = −2uv +2 − − ε11 u v +2 ρ ∂x ∂y ρ ∂y ∂y Dt Dv 2 = Dt 0 Dw2 = Dt 0 ( ) Anisotropy ∑= 0 ∑ = 2ε ∑ = k − equation Anisotropy in simple shear Homogeneous shear Channel flow ¾ Development in time of stresses normalized by k Strain rate x time ¾ Normal and shear stresses The importance of anisotropy: expansion (deceleration) Positive generation Du 2 2 ∂U = −2u + .... Dt ∂x Dv ∂U = v2 + .... Dt ∂x Negative generation 2 Dw2 ∂U = w2 + .... Dt ∂x 1 Dk 2 2 2 ∂U = − (2u − v − w ) + ..... 2 Dt ∂x Low or negative k-production, relative to very high EVM production Eddy viscosity k ⎛ ∂U ⎞ ε ⎜⎝ ∂x ⎟⎠ 2 2Cμ 2 Anisotropy in expansion and contraction aij ≡ ui u j 2 − δ ij k 3 S * = k / ε Sij Sij / 2 Round impinging jet Wall-normal stress and mean velocity Anisotropy in plain strain Other sensitivities Strong effect of curvature on anisotropy and shear stress. ∂V ∂x y x Strong effects of rotation on anisotropy and shear stress y x S P Inapplicability of Fourier-Fick law in scalar transport ¾ Production of flux vector: ∂Φ ∂U i Puiφ = −uiuk − uiϕ ∂ xk ∂ xk Reynolds-Stress-Transport Modelling Closure of exact stress-transport equations ∂U j ⎧ ∂U i ⎫ = − ⎨ uiuk + u j uk ⎬ + ( Pressure − velocity ) Dt3 ∂ xk ∂ xk ⎭ ⎩1444 12 4 24444 3 = AdvectiveTransport Duiu j Cij Pij = Production + Diffusion − Dissipation Pressure-velocity, dissipation and diffusion require approximation About 10-15 major closures forms Modern closure aims at realisability, 2-component limit, coping with strong inhomogeneity and compressibility Additional equations for dissipation tensor ε ij At least 7 pde’s in 3D (up to 17 in heat/scalar transport) Numerically difficult in complex geometries and flow Can be costly Dissipation and pressure-velocity are major sources of error The exact dissipation-rate equation Modelled dissipation-rate equation ⎡⎛ ∂U i k k εε% ⎞ ∂ε ⎤ − Cε 2 ⎢⎜νδ kl + Cε uk ul ⎟ ⎥ − Cε 1 ui u j ∂x j k ε ε ⎠ ∂xl ⎦ ⎣⎝ + "special" model fragments ∂ Dε = Dt ∂xk Generalised gradient-diffusion approximation ε {Cε 1 (Production of k ) − Cε 2 fε (Dissipation of k )} k In energy equilibrium, Pk = ε , and the imbalance is absorbs by diffusion Transport equations for ε ij are too complex as basis for modelling Anisotropy in dissipation – algebraic approximations of the form: ui u j 2 ε ij = f e εδ ij + (1 − f e ) ε 3 k In most models, f e = 1 reflecting assumption of small-scale isotropy Closure – stress diffusion Regarded as least influential (suggested by DNS/LES). Represented as gradient-diffusion with tensorial diffusivity. Simplest model: ∂ Diffij = − ∂xk ⎧⎪ k ∂u j u j ⎫⎪ ⎨cd uk u m ⎬ x ∂ ε m ⎪ ⎪⎩ ⎭ Based on observation that the most important fragment in the exact diffusion term is uk ui u j . It can be shown, via transport equations for triple correlation, uk ui u j , that the production of these triple correlations is by gradients of stresses of the form Pijk = −uk u m ∂u j u j ∂xm + ..... Suggests (also on dimensional grounds) Diffij = − ∂ c(time scale) × (production ijk )} { ∂xk Closure – pressure-strain / velocity Extremely important: responsible for redistribution among normal stresses. Regarded as the hardest term to model Pressure-velocity dictates energy transfer and hence v But 2 v 2dictates uv Duv p ⎛ ∂u ∂v ⎞ ∂ ⎛ 2 pu ⎞ μ ∂uv 2 ∂U = −v + ⎜ + ⎟ − ⎜ uv + − ε12 ⎟+ ρ ⎠ ρ ∂y ∂y ρ ⎝ ∂y ∂x ⎠ ∂y ⎝ Dt ⎫ μ ∂u 2 ∂U Du 2 p ∂u ∂ 2 = −2uv +2 − − ε11 u v +2 ⎪ ρ ∂x ∂y ρ ∂y ∂y Dt ⎪ ⎪ μ ∂v 2 Dv 2 p ∂v ∂ ⎛ 3 2 pv ⎞ ⎪ = +2 − ⎜v + + − 0 2 ε ⎟ 22 ⎬ ρ ∂y ∂y ⎝ ρ ⎠ ρ ∂y Dt ⎪ ⎪ 2 μ ∂ Dw2 p ∂w ∂ w ⎪ − ε 33 = +2 − 0 vw2 + 2 ⎪⎭ ρ ∂z ∂x ρ ∂y Dt ( ) ( ) Closure – pressure-strain Subject to constrains: ¾ Isotropisation: transfer of energy from largest stress to lower ones ¾ Inhibition of isotropisation at walls/interfaces (splatting, reflection) ¾ shear stresses have to decline as isotropisation progresses Guidance provided by ‘exact’ integration for pressure-fluctuations and substitution in pressure-velocity correlation p ⎛ ∂ui ∂u j ⎞ Aij Φ ij = ⎜ + ⎟⎟ ⎜ ρ ⎝ ∂x j ∂xi ⎠ * ⎧⎛ 2 1 ⎪ ∂ ul um ⎞ ⎛ ∂ui ∂u j + = ⎨⎜ ⎟ ⎜⎜ ∫ 4π V ⎪⎝ ∂xl ∂xm ⎠ ⎝ ∂x j ∂xi ⎩ V x ⎞ ⎫⎪ dV (x* ) ⎟⎟ ⎬ * − x x ⎠ ⎪⎭ * *⎫ ⎧ 1 ⎪ ⎛ ∂um ⎞ ⎛ ∂ui ∂u j ⎞ ⎛ ∂U l ⎞ ⎪ dV (x* ) + + ⎟⎜ ⎨2⎜ ⎟ ⎬ ⎟ ⎜ 4π V∫ ⎪ ⎝ ∂xl ⎠ ⎜⎝ ∂x j ∂xi ⎟⎠ ⎝ ∂xm ⎠ ⎪ x − x* ⎭ ⎩ Bijkl x* + body-force and surface terms Closure – pressure-strain Suggests the general Ansatz: ∂U k = A a + k B a ε { } { } Φ ij ij ij ijkl ij ∂xl ⎧⎪ ui u j 2 ⎫⎪ − δ ij ⎬ ⎨aij ≡ k 3 ⎪⎭ ⎪⎩ (+ body-force and wall terms) Most complex model is cubic Much more popular is the quasi-linear form ε⎛ 2 2 ⎞ ⎛ ⎞ Φ ij = −C1 ⎜ ui u j − δ ij k ⎟ − C2 ⎜ Pij − δ ij Pk ⎟ 3 3 k⎝ ⎠ ⎝ ⎠ (+ body-force and wall terms) This is a sink term in the second-moment equations, depressing anisotropy in proportion to anisotropy of stresses and productions Ensures that anisotropy in stresses and productions drives energy from above-average normal stresses to below-average ones Coefficients sensitized to anisotropy invariants, turbulence Reynolds number…..in lieu of non-linear expansions Closure – pressure-strain Model Performance Construction and calibration rely heavily on highly-resolved experimental & simulation data Done mostly by reference to thin-shear-flow data Models work well for many flows Notable exception: flow separating from curved surfaces (2d & 3d) Associated with dynamics of highly unsteady separation (& preseparation) U-velocity close to wall Separation from curved surface Shear-stress profiles with different models relative to LES 3 (x/h=2.0) (x/h=6.0) LS- ε AL- ε WJ- ω CLS- ε AJL-ω LES LS- ε AL- ε WJ- ω CLS- ε AJL-ω LES 2.5 2 1.5 1 0.5 0 -0.04 -0.02 uv/U2b 0 0.02 -0.04 -0.02 uv/U2b 0 0.02 Model developments Model defects are difficult to cure, but efforts are ongoing Example: re-examination of dissipation and pressure-velocity interaction terms in separation from curved ramp Foundation: highly-resolved simulation – near DNS, 25M nodes ReH=13700; ReΘ = 1150 Second moments, invariants, budgets of all second moments…. Part of larger study on separation control with synthetic jets Experimental data Starting point Choice of basic model, based on full computation LES: (x/h)s = 0.87 & (x/h)r=4.21 JH: (x/h)s = 0.79 & (x/h)r=4.15 Used to illustrate path to improvement 1.09 separation Shima: (x/h)s = 0.57 & (x/h)r=5.60 SSG+C: (x/h)s = 0.84/1.50 & (x/h)r=1.11/4.10 Secondary recirculation Defect identification Focus on shear stress in separated shear layer Defect identification Budgets for uv and uu uv - LES uu - LES uv - JH uu - JH Model fragmentation - dissipation A-priori study of dissipation-rate equation Isolated solution of equation LES strains and stresses input into equation Only output is dissipation Examination of a range of corrections in efforts to procure agreement with LES data for dissipation rate Model fragmentation - dissipation Ongoing efforts to sensitize dissipation to mean-flow/turbulence length scales Model fragmentation – dissipation components A-priori study of dissipation anisotropy – stresses and ε from LES into 2 Various proposals 3 f s = 1 − AE ε ui u j + (ui u j n j nk + u j u k ni nk + uk u l n j nk ni n j ) f d * ε ij = f d = (1 + 0.1Ret ) −1 3 u p uq k n p nq f d 1+ 2 k ε ij = f sε ij* + (1 − f s ) δ ijε Weighting function sensitized to anisotropy invariant A = 1 − 9 ( A2 − A3 ) 8 A2 = aij aij ; A3 = aij a jk aki ; aij = Use of anisotropy invariant ui u j 2 − δ ij k 3 Model fragmentation – dissipation components Component ε22 Model fragmentation – pressure-velocity Quasi-linear approximation ε⎛ 2 2 ⎞ ⎛ ⎞ Φ ij = −C1 ⎜ ui u j − δ ij k ⎟ − C2 ⎜ Pij − δ ij Pk ⎟ (+ wall-reflection terms) k⎝ 3 3 ⎠ ⎝ ⎠ Coefficients sensitized to anisotropy invariants, in compensation to the omission of high-order fragments C1 = C + AE C = 2.5 A[min{0.6, A2 }]1/ 4 f ⎧⎪⎛ Ret ⎞3/ 2 ⎫⎪ f = min ⎨⎜ ⎟ ,1⎬ 150 ⎠ ⎪⎩⎝ ⎪⎭ C2 = 0.8 A1/ 2 Model fragmentation – pressure-velocity Sensitivity of coefficients to pressure-velocity interaction of uu Shock-induced Separation on 3D Jet-Afterbody - RSTM General view and surface-pressure coefficient ReL =2x107 0.5 M nodes CFL=O(1000) Desktop workstation, a few CPU hours Shock-induced Separation on flat plate – LES Touber and Sandham, 2010 Reτ =3000, M=2.3 20 M nodes, 240,000 CPU hours Concluding remarks Fundamentally, Second-moment closure is far superior to eddyviscosity modelling. In reality, closure is extremely challenging, because the anisotropy is an extremely influential model element and is difficult to approximate. Redistribution and dissipation are especially influential. Many ways of construction models, but all involve calibration. Does involve “curve-fitting”, but is based on rational principles and physically tenable assumptions. Little used, because of “the-simpler-the-better” attitude. Second-moment closure is inappropriately complex in (most) thin shear flows, but the only fundamentally solid approach in complex strain.