Exactly Divergence-free Hybrid Discontinuous Galerkin Method for Incompressible Turbulent Flows
TTechnische Universit¨at M¨unchen
Department of InformaticsMaster’s Thesis in Computational Science and Engineering
Exactly Divergence-free HybridDiscontinuous Galerkin Method forIncompressible Turbulent Flows
Xaver Mooslechner BSc a r X i v : . [ c s . C E ] S e p echnische Universit¨at M¨unchen Department of InformaticsMaster’s Thesis in Computational Science and Engineering
Exakt divergenzfreie hybridediskontinuirliche Galerkin Verfahren f¨ur inkompressible turbulente Str¨omungenExactly Divergence-free HybridDiscontinuous Galerkin Method forIncompressible Turbulent Flows
Author: Xaver Mooslechner BSc1 st Examiner: Univ.Prof. Dr. Folkmar Bornemann2 nd Examiner: Univ.Prof. Dr. Joachim Sch¨oberl, Technische Universit¨at WienAssistant advisor: Dr. Philip Lederer, Technische Universit¨at WienSubmission date: 22.04.2020 isclaimer
I hereby declare that this thesis is entirely the result of my own work except where other-wise indicated. I have only used the resources given in the list of references.Vienna, 22.04.2020 cknowledgement
I wish to express my gratitude to my supervisor, Prof. Dr. Joachim Sch¨oberl, for providingme the opportunity to write this master thesis within his work group at the technicaluniversity of Vienna. I want to thank you for the incredible amount he has taught meduring this work. At the same time, I want to express my appreciation to Prof. Dr.Folkmar Bornemann for taking the official supervision of my thesis on him.I also want to thank Dr. Philip Lederer for introducing and enhancing my knowledge inthis topic and leading me to find solutions for arising problems.Furthermore, I want to thank my girlfriend Caroline for all the support during my studies. urzfassung
Diese Arbeit besch¨aftigt sich mit der Untersuchung eines H (div)-konformen hybriden diskon-tinuirlichen Galerkin Verfahrens f¨ur inkompressible turbulente Str¨omungen.Die Diskretisierungsmethode liefert einige g¨unstige physikalische und l¨osungs-orientierendeEigenschaften, welche von Vorteil sein k¨onnen f¨ur das Aufl¨osen von rechenintensiven tur-bulenten Strukturen. Eine herk¨ommliche Methode zur Diskretisierung der Navier-StokesGleichungen mit den bekannten Taylor-Hood Elementen ist auch gegeben, um einen Ver-gleich wiedergeben zu k¨onnen. Die vier Hauptmethodiken zur Simulation von turbulen-ten Str¨omungen sind erl¨autert: die Reynolds-gemittelte Navier-Stokes Simulation, dieGrobstruktursimulation, die Methode zur Mehrskalenvariationsrechnung und die direktenumerische Simulation. Die Grobstruktursimulation und Mehrskalenvariationsrechnungzeigen gute Ergebnisse in der Berechnung von traditionell schwierigen turbulenten Str¨omungs-f¨allen. Diese Genauigkeit kann nur durch direktes berechnen der Navier-Stokes Gleichungen¨ubertroffen werden, was jedoch mit sehr hohen Kosten verbunden ist. Eine sehr verbre-itete Herangehensweise ist die Reynolds-Mittelung, da diese sehr kosteng¨unstig ist. DiesePrinzipien wurden an beiden Diskretisierungsverfahren angewendet und validiert anhandder turbulenten Kanalstr¨omung.Alle numerischen Simulationen wurden mit Hilfe der finiten Elementen Bibliothek Net-gen/NGSolve durchgef¨uhrt. bstract This thesis deals with the investigation of a H (div)-conforming hybrid discontinuous Galerkindiscretization for incompressible turbulent flows.The discretization method provides many physical and solving-oriented properties, whichmay be advantageous for resolving computationally intensive turbulent structures. A stan-dard continuous Galerkin discretization for the Navier-Stokes equations with the well-known Taylor-Hood elements is also introduced in order to provide a comparison. Thefour different main principles of simulating turbulent flows are explained: the Reynolds-averaged Navier-Stokes simulation, large eddy simulation, variational multiscale methodand the direct numerical simulation. The large eddy simulation and variational multiscalehave shown good promise in the computation of traditionally difficult turbulent cases. Thisaccuracy can be only surpassed by directly solving the Navier-Stokes equations, but comeswith excessively high computational costs. The very common strategy is the Reynolds-average approach, since it is the most cost-effective. Those modelling principles have beenapplied to the two discretization techniques and validated through the basic plane channelflow test case.All numerical tests have been conducted with the finite element library Netgen/NGSolve. ontents K − (cid:15) model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.2.3 The K − ω model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294.2.4 The K − ω SST model . . . . . . . . . . . . . . . . . . . . . . . . . . 304.2.5 The scalar convection diffusion problem . . . . . . . . . . . . . . . . 314.2.6 Near-wall treatments . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.3 Large eddy simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.3.1 Sub-grid scale modelling . . . . . . . . . . . . . . . . . . . . . . . . . 374.4 Variational multiscale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 ontents H (div)-conforming HDG method for the Navier-Stokes problem 556.3 Low order H (div)-conforming finite elements . . . . . . . . . . . . . . . . . 586.4 HDG discretization of the turbulence models . . . . . . . . . . . . . . . . . 60 List of Figures 85List of Tables 87References 88 ii Introduction
The phenomena of turbulence in fluid flow is one of the most impressing problems of classicmechanics. Since it has been initially observed and described by Leonardo da Vinci in the16th century, a lot of effort was successively dedicated to understand the emergence ofturbulence and their structures. The very chaotic, irregular and apparently unpredictablebehavior of turbulent flows leads to a challenging subject of study. Nevertheless, this phe-nomena is by far not unknown and occurs very frequently in our daily life, from pouringmilk in a cup of coffee to circulation of air in the atmosphere. Especially for many scien-tific and engineering purposes, turbulence plays a vital role and its prediction is of greatinterest.In the past, due to the complexity of the governing equations of fluid motion, the ana-lytical approach is highly restricted to simple flow cases and can by far not be applied toturbulence. As a result of this problem, the experimental analysis was the only methodto deal with this phenomena. By the rapid growth of computational power and reductionof its operational cost, the numerical approach is becoming more and more advantageouscompared to expensive experiments. Since the past few decades, this benefit drives the sci-entific community to establish methods for studying turbulent flows. The most importantaspects of the development of methods in computational fluid dynamics (CFD) is accuracyof the solution and cost effective algorithms. In order to accurately simulate turbulenceover an appropriate time period, different simulation techniques and discretization methodsneed to be developed.
In the first chapter, the focus lies on the theory behind fluid dynamics and the derivationof the governing equations, the Navier-Stokes equation. This section should reacquaint thebasics, which are indispensable for the upcoming chapters.In chapter two, an overview of turbulence is given. The most relevant fields of interest aredescribed in more detail and should provide sufficient knowledge. The turbulence sectionis divided in stochastic, spectral and theory description.The third chapter is dedicated to the different simulation principles and modelling. Itprovides the three main approaches how to numerically compute turbulent flow. In thelast years, a rather new method has emerged in the field of CFD, which is described in thissection as well. All of the four principles are analyzed and have been conducted with thedifferent types of discretization techniques.The two different spatial discretizations of the famous saddle-point problem are expressed inthe fourth and fifth chapter of the thesis. One is the standard Continuous Galerkin (CG)method with the well-known elements from the Taylor-Hood family, which ensures only1
Introduction discrete divergence-free velocity property. The other technique is composed of a HybridDiscontinuous Galerkin (HDG) method, developed by Christoph Lehrenfeld and JoachimSch¨oberl in [LS16]. This mixed method guarantees an exact divergence-free flow and leadsto an appropriate physical description.The last chapter briefly shows the results of the performed simulations and further ondiscusses the differences.
All numerical examples were implemented and tested in the finite element library Net-gen/NGSolve, see [Sch97] and [Sch14]. 2
Introduction
List of Abbreviations
The following abbreviations are used in this thesis:CFD computational fluid dynamicsCG continuous GalerkinDG discontinuous GalerkinDNS direct numerical simulationDOF degree of freedomGS grid scaleHDG hybrid discontinuous GalerkinIMEX implicit-explicitLBB Ladyshenskaja-Babuˇska-BrezziLES large eddy simulationPDE partial differential equationRANS Reynolds-averaged Navier-StokesRST Reynolds stress tensorRTT Reynolds-transport-thoeremSGS sub-grid scaleVMS variational multiscale 3
Derivation of the Navier-Stokes equations
A fluid element represents an accumulation of fluid molecules within an infinitesimal smallvolume dV , where its averaged motion over an infinitesimal time interval dt . Each elementhas its density ρ and velocity u . This macroscopic view leads to the continuum description,where the characteristics of the fluid element are prescribed via partial differential equations(PDE).In the following, we assume the physical domain Ω ⊂ R and the fixed time interval denotedby [0 , T ]. Further on, we assume the fluid to be Newtonian, pure and viscous. The followingrelevant quantities are shown in Table 2.1.Description Quantity Unitvelocity u ( x, t ) ∈ C (Ω × [0 , T ] , R ) m s − density ρ ( x, t ) ∈ C (Ω × [0 , T ] , R ) kg m − pressure p ( x, t ) ∈ C (Ω × [0 , T ] , R ) kg m − s − force F ( x, t ) ∈ C (Ω × [0 , T ] , R ) kg m s − kinematic viscosity ν ∈ R \ { } m s − Table 2.1: Physical quantitiesWe proceed with the following derivations as in [OBR15] and [Kuh14].
In classical field theories, the kinematic of fluid flow can be basically described by twodifferent point of views.The Lagrangian specification of the fluid field is one way of looking at a fixed fluid elementas it moves through space and time. The observer of the fluid parcel monitors the changeof its properties (e.g. velocity). At a specific time t , the fluid element has the position ξ ( t ) = ξ . Since only one fluid parcel may stay at a certain position, the element is”labeled” ξ . Therefore ξ ( ξ ; t ) is the position vector and describes the trajectory of thefluid parcel u ( ξ ; t ) = dξ ( ξ ; t ) dt . (2.1.1)The Eulerian description is a way of focusing at fluid flow on a specific location x in spaceas the fluid passes through time. As a result, the change of quantities at a fixed positionin space depends on the change of each individual fluid element as well as the change offluid elements at position x .The link between the Lagrangian and Eulerian specifications is given by the material deriva-tive. Suppose the quantity φ ( x, t ) ∈ C (Ω × [0 , T ] , R ) from the Eulerian point of view andsubstitute x = ξ ( ξ ; t ) from the Lagrangian approach, it follows φ ( ξ , ξ , ξ , t ). The total4 Derivation of the Navier-Stokes equations rate of change with respect of time of φ is dφdt = ∂φ∂t + ∂φ∂x ∂ξ ∂t + ∂φ∂x ∂ξ ∂t + ∂φ∂x ∂ξ ∂t = ∂φ∂t + ∇ φ · ∂ξ∂t . Using Equation (2.1.1), it follows:
Definition 1 ( Material derivative ) . Let the quantitiy φ ( x, t ) ∈ C (Ω × [0 , T ] , R ) bean arbitrary function in a fluid transported with velocity u. Then the material derivativeis defined by DφDt = ∂φ∂t + ( u · ∇ ) φ. (2.1.2)So far, a single fluid element with infinitesimal expansion was considered. Now we extentour consideration to a dense pack of many fluid elements, a fluid volume. The relation be-tween Lagrangian and Eulerian view for a fluid volume is given by the Reynolds-transport-theorem (RTT). Definition 2 ( Reynolds-transport-theorem ) . Let V ( t ) ⊂ Ω be an arbitrary vol-ume with fixed number of fluid elements (fixed mass) and the quantitiy φ ( x, t ) ∈ C (Ω × [0 , T ] , R ) be an arbitrary function in a fluid transported with velocity u. Thenthe Reynolds-transport-theorem is defined by DDt (cid:90) V ( t ) φ ( x, t ) dx = (cid:90) V ( t ) ∂φ ( x, t ) ∂t dx + (cid:90) ∂V ( t ) φ ( x, t ) u ( x, t ) · n ds, (2.1.3) while the surface ∂V ( t ) defines the boundary of V ( t ) and n the outward directed normalvector of the surface ∂V ( t ) . Using Gauß’ theorem for Equation (2.1.3) yields
DDt (cid:90) V ( t ) φ ( x, t ) dx = (cid:90) V ( t ) (cid:18) ∂φ∂t + ∇ · ( φu ) (cid:19) dx. (2.1.4) Again let define V ( t ) ⊂ Ω be an arbitrary volume with fixed mass and take density as ourquantity, the total mass at time t is m ( t ) = (cid:90) V ( t ) ρ ( x, t ) dx. As the mass has has to be conserved, the rate of change with respect to time has to be
DmDt = DDt (cid:90) V ( t ) ρ ( x, t ) dx = 0 . With use of Equation (2.1.3), the conservation of mass reads
DmDt = (cid:90) V ( t ) ∂ρ∂t dx + (cid:90) ∂V ( t ) ρu · n ds = 0 . (2.2.1)5 Derivation of the Navier-Stokes equations
The first term of Equation (2.2.1) on the right-hand side represents the change in massdue to rate of change of density with respect to time and the second term gives balanceof in- and outcoming mass flux. Further using Equation (2.1.4) yields to the integral anddifferential form of the continuity equation (cid:90) V ( t ) (cid:18) ∂ρ∂t + ∇ · ( ρu ) (cid:19) dx = 0 , (2.2.2) ∂ρ∂t + ∇ · ( ρu ) = 0 . (2.2.3)In this thesis, we strictly focus on incompressible fluid flow and ρ = const , therefore DρDt = 0 , (2.2.4)substituting Equation (2.2.4) into Equation (2.2.3) leads to the incompressibility constraintfor the velocity ∇ · u = 0 . (2.2.5) By using Newton’s second law, the conservation of momentum equation can be deducted.It follows that the change of momentum of a fluid volume (with fixed mass) equals theresulting forces affecting the fluid volume. The momentum vector is given by P = (cid:90) V ( t ) ρ ( x, t ) u ( x, t ) dx. Therefore
DPDt = DDt (cid:90) V ( t ) ρ ( x, t ) u ( x, t ) dx = (cid:88) m F m . Possible forces are: • Volume forces (e.g. gravity) affect the control volume V ( t ) by F V = (cid:90) V ( t ) ρf dx, (2.3.1) f presents the acceleration in vector form. • Surface forces (e.g. shear stress) acting on the surface of the control volume V ( t ) by F S = (cid:90) S ( t ) σ · n ds, (2.3.2) σ is a symmetric tensor of second order σ = − pI + τ = − p + τ τ τ τ − p + τ τ τ τ − p + τ . Derivation of the Navier-Stokes equations
The diagonal components of τ present the normal stresses, while the other compo-nents describe the shear stresses. By the use of the Gauß’ theorem, the term inEquation (2.3.2) can be reforumlated to (cid:90) S ( t ) σ · n ds = (cid:90) V ( t ) ∇ · σ dx. (2.3.3)The divergence of a matrix as shown in Equation (2.3.3) is taken row-wise.The left-hand side of the momentum equation is given DDt (cid:90) V ( t ) ρu dx = (cid:90) V ( t ) ∂ρu∂t dx + (cid:90) S ( t ) ρu ( u · n ) ds = (cid:90) V ( t ) (cid:18) ∂ρu∂t + ∇ · ( ρu ⊗ u ) (cid:19) dx. Here, u ⊗ u indicates the outer (dyadic) product of two vectors.The right-hand side of the momentum equation consists of (cid:90) V ( t ) (cid:18) ∇ · σ + ρf (cid:19) dx, while σ may be decomposed to − pI + τ , resulting in (cid:90) V ( t ) (cid:18) − ∇ p + ∇ · τ + ρf (cid:19) dx. The final momentum equation in integral and differential form is then (cid:90) V ( t ) (cid:18) ∂ρu∂t + ∇ · ( ρu ⊗ u ) (cid:19) dx = (cid:90) V ( t ) (cid:18) − ∇ p + ∇ · τ + ρf (cid:19) dx, (2.3.4) ∂ρu∂t + ∇ · ( ρu ⊗ u ) = −∇ p + ∇ · τ + ρf . (2.3.5)Since we are assuming incompressible flows, the constraint from Equation (2.2.4) and (2.2.5)further simplifies the conservation of momentum equation. Exploiting the identity ∇ · ( u ⊗ u ) = ( u · ∇ ) u + u ( ∇ · u ) . (2.3.6)As a result, the incompressible momentum equation in differential form reads as follows ∂u∂t + ( u · ∇ ) u = − ρ ∇ p + 1 ρ ∇ · τ + f . (2.3.7)7 Derivation of the Navier-Stokes equations
In the following, further assumptions have to be made to generalize the stress tensor σ . Inthis thesis, we assume that the fluid is a Stokes fluid. For this sake, the properties of suchfluid is given: • Conservation of angular momentum, means that σ is symmetric σ = σ T . • σ is isotropic. • If the fluid rests or moves due to rigid-body motion, it holds σ = − pI .The general ansatz for the stress tensor, assuming viscous Newtonian fluid and incompress-ible flow, is τ = µ ( ∇ u + ∇ u T ) , while ∇ u denotes the vector gradient ∇ u = ( ∇ ⊗ u ) T .The divergence of the stress tensor can be further simplified if µ = const to ∇ · τ = µ ∆ u. Additionally, for the Navier-Stokes model one needs suitable initial and boundary condi-tions for u . The initial condition is u ( x,
0) = u ( x ), which has to be physically correct aswell as divergence free. The boundary conditions are distinguished between Dirichlet, Neu-mann or Robin boundary. Let ∂ Ω be subdivided into three parts ∂ Ω = ∂ Ω D ∪ ∂ Ω N ∪ ∂ Ω R with ∂ Ω D ∩ ∂ Ω N ∩ ∂ Ω R = ∅ . Dirichlet conditions are of the form u ( x, t ) = u D ( x ) for x ∈ ∂ Ω D and in applications often used for inflow or no-slip boundary conditions. In orderto describe outflow conditions, the Neumann boundary conditions of the form ∂u∂n = g ( x )for x ∈ ∂ Ω N are appropriate. The Robin boundary condition is a linear combination ofthe Dirichlet and Neumann condition. Since it is not a major issue in this thesis, this thirdtype boundary condition will be neglected.Finally, the Navier-Stokes problem is complete and can be read as: Problem 1 ( Incompressible Navier-Stokes equations ) . Let ν, ρ ∈ R \ { } , f ∈ C (Ω , R ) , u D ∈ C ( ∂ Ω D , R ) , g ∈ C ( ∂ Ω N , R ) and u ∈ C (Ω , R ) find u ∈ C (Ω × [0 , T ] , R ) and p ∈ C (Ω × [0 , T ] , R ) such that ∇ · u = 0 in Ω × [0 , T ] , (2.4.1) ∂u∂t + ( u · ∇ ) u = − ρ ∇ p + ν ∆ u + f in Ω × [0 , T ] , (2.4.2) u = u in Ω , t = 0 ,u = u D in ∂ Ω D × [0 , T ] ,∂u∂n = g in ∂ Ω N × [0 , T ] . Derivation of the Navier-Stokes equations
As it is of indispensible importance, a specfic dimensionaless parameter can be deductedvia the dimensionaless Navier-Stokes equations. The following variables are: u ∗ = uU , p ∗ = pρU ,x ∗ = xL , t ∗ = t UL , while U and L are the corresponding characteristic length and velocity scales. If we insertthe dimensionless variables into the continuity and momentum equation from Problem 1(neglecting the volume force term), we get the Navier-Stokes equations in dimensionlessform ∇ · u ∗ = 0 , (2.4.3) ∂u ∗ ∂t ∗ + ( u ∗ · ∇ ) u ∗ = −∇ p ∗ + 1 Re ∆ u ∗ . (2.4.4)The new parameter Re = ULν is called Reynolds number and gives the ratio between inertialand viscous forces and is a measure for turbulence.9
Turbulence phenomena
Da Vinci’s observations and drawings primarily describe that turbulent flow regime con-tains of eddying motion and structures of whirls in the 16th century.The experiment of Osborne Reynolds was one of the first attempts to experimentally quan-tify turbulence and his work laid the foundation of turbulence theory. Reynolds showed thattwo flow regimes exist, laminar and turbulent, and introduced a parameter (the Reynoldsnumber) to distinguish between those states of fluid flow. The transition of laminar toturbulent flow occurs only if a certain critical Reynolds number has been exceeded.Lewis F. Richardson firstly expressed the concept of energy cascade in well-developed tur-bulence and captured his observations in his famous poem in the 1920s. In the regime ofturbulent flow, a wide range of different length scales exists, from large whirls to smallerones. The biggest eddies of size of almost the characteristic length scale contain the most ofkinetic energy, while the smaller whirls get supplied by the larger ones. Evoking a cascadingwaterfall until a certain length scale is reached, where the remaining energy dissipates dueto viscosity.Around 20 years later, Andrei N. Kolmogorov extended the work of Richardson and evolvedthe theory of turbulence and the concept of energy spectrum. This spectrum gives the dis-tribution of energy among turbulence vortices as function of vortex size. In his analysis, hefound that these scales are well separated, where the intermediate subrange of scales arestatistically isotropic and Kolmogorov hypothesized a universal form of the energy spec-trum in this region.Over the last sixty years since the publishing of Kolmogorov’s theory, much progress hasbeen made in the field of turbulence flow. Particularly, the defining and identifying ofcoherent turbulent structures, vortex structures which persist in the flow for a relativelylong time, through experiments have given great insight. A significant contribution to thisprogress is due to the development of advanced numerical simulation methods. These newtechniques have made available data, which could not been measured in any experimentalinvestigation.
However, defining turbulence is by no means easy. It is common to describe it by listingits characteristics. Further detailed description may be found in [Hin75] and [Pop00]. • Turbulent flows are chaotic . In the sense of small disturbances in the initial field,which will be amplified and leading to an uncorrelated flow field. This makes a de-terministic approach impractical and intractable to describe its motion in full detailsas function of space and time. The random fluctuations may have amplitudes of tento thirty percent of its mean value. 10
Turbulence phenomena • Turbulent flows are unsteady . In the regime of turbulence, a true stationary solutiondoes not exist due to its high irregularity. Only stochastic steady form can be reached. • Turbulent flows are rotational . It has been shown many times that turbulence onlyarises and persists in rotational flows, in the presence of shear. An initially irrotationalflow may become rotational by laminar-turbulent transition. Such process may onlyhappen when inertial forces dominate viscous effects (at high Reynolds numbers) andsmall perturbations are no longer be damped by molecular viscosity. • Turbulent flows are diffusive . Such flows cause very rapid mixing of any transportedquantity like momentum or heat. The turbulent diffusion allows much faster mixingof quantities than if only molecular diffusion processes were involved. Diffusivity maybe important from a practical point of view, for instance, the turbulent drag on anairfoil. • Turbulence is still a phenomenon of continuum mechanics . Even the smallest eddiesoccurring in turbulent flow are far greater than any single fluid element (molecularscale).
In general, a detailed description of flow quantities in time and space is neither advisable nordesirable. In order to be able to compare two different turbulent flows, it only makes senseif initial and boundary conditions match, therefore an appropriate statistical descriptionis desired. We usual refer to a statistical representation of the fluctuations. A turbulentfluid flow is called statistically steady, if the averaged quantities of two far separated timeinstants are equal. Osborne Reynolds firstly occupies this issue and introduced the Reynoldsdecomposition and the Reynolds operator. He decomposed the flow quantities in its meanvalue and fluctuating part u ( x, t ) = (cid:104) u (cid:105) ( x ) + u (cid:48) ( x, t ) , p ( x, t ) = (cid:104) p (cid:105) ( x ) + p (cid:48) ( x, t ) . (3.3.1)The mean value is defined via the time-average operator (cid:104) φ (cid:105) ( x ) = lim T →∞ T T (cid:90) φ ( x, t ) dt, (3.3.2)which is a Reynolds operator.An ensemble average (stochastic mean) of a random variable u ( x, t ) calculated from N independent realizations of the same phenomenon is defined as (cid:104) u (cid:105) ( x, t ) = lim N →∞ N N (cid:88) n =1 u ( n ) ( x, t ) . (3.3.3)An ensemble { u ( n ) ( x, t ) } Nn =1 is a collection of notionally identical experiments. Since theflow is turbulent, the fluid motion differs from each instance of the ensemble, because mi-croscopically differences in the (experimental) setup become significant as time progresses.11 Turbulence phenomena
Nevertheless, it can be shown that for a very long time period T → ∞ each realizationof the ensemble can be considered as a representative of all possible realizations of anensemble. In equilibrium systems, time and ensemble averages of physical quantities areequivalent due to ergodic principles.In the following, a Reynolds operator is defined by satisfiying certain properties. Definition 3 ( Reynolds operator ) . Let φ, ψ ∈ C (Ω × [0 , T ] , R ) and c ∈ R , then thefollowing properties are satisfied by the Reynolds operator:(i) (cid:104) φ + ψ (cid:105) = (cid:104) φ (cid:105) + (cid:104) ψ (cid:105) (ii) (cid:104) cφ (cid:105) = c (cid:104) φ (cid:105) (iii) (cid:104) c (cid:105) = c (iv) (cid:104)(cid:104) φ (cid:105) ψ (cid:105) = (cid:104) φ (cid:105)(cid:104) ψ (cid:105) (v) (cid:104) ∂φ∂x (cid:105) = ∂ (cid:104) φ (cid:105) ∂x (vi) (cid:104) ∂φ∂t (cid:105) = ∂ (cid:104) φ (cid:105) ∂t Hence, these conditions ensue:(vii) (cid:104)(cid:104) φ (cid:105)(cid:105) = (cid:104) φ (cid:105) (viii) (cid:104) φ (cid:48) (cid:105) = 0 , φ (cid:48) = φ − (cid:104) φ (cid:105) (ix) (cid:104)(cid:104) φ (cid:105)(cid:104) ψ (cid:105)(cid:105) = (cid:104) φ (cid:105)(cid:104) ψ (cid:105) However, fluctuation moments of second or higher order are not necessarily zero. The stan-dard deviation is often called root-mean-square velocity in turbulence flow and is definedas u i,rms = (cid:113) (cid:104) u (cid:48) i u (cid:48) i (cid:105) = (cid:112) (cid:104) ( u i − (cid:104) u i (cid:105) ) (cid:105) , u rms = (cid:115) (cid:88) i (cid:104) u (cid:48) i u (cid:48) i (cid:105) , (3.3.4)while u i for i = 1 , , (cid:28) ∇ · ( (cid:104) u (cid:105) + u (cid:48) ) (cid:29) = 0 , (cid:28) ∂ ( (cid:104) u (cid:105) + u (cid:48) ) ∂t + (cid:0) ( (cid:104) u (cid:105) + u (cid:48) ) · ∇ (cid:1) ( (cid:104) u (cid:105) + u (cid:48) ) (cid:29) = (cid:28) − ρ ∇ ( (cid:104) p (cid:105) + p (cid:48) ) + ν ∆( (cid:104) u (cid:105) + u (cid:48) ) + f (cid:29) . Using the properties from Definition 3, it can be further simplified to ∇ · (cid:104) u (cid:105) = 0 ,∂ (cid:104) u (cid:105) ∂t + (cid:28)(cid:0) ( (cid:104) u (cid:105) + u (cid:48) ) · ∇ (cid:1) ( (cid:104) u (cid:105) + u (cid:48) ) (cid:29) = − ρ ∇(cid:104) p (cid:105) + ν ∆ (cid:104) u (cid:105) + (cid:104) f (cid:105) . The second term of the left-hand side of the momentum equation is the nonlinear convectiveterm, which needs particular considerations. It follows (cid:28)(cid:0) ( (cid:104) u (cid:105) + u (cid:48) ) · ∇ (cid:1) ( (cid:104) u (cid:105) + u (cid:48) ) (cid:29) = (cid:28) ∇ · (cid:0) ( (cid:104) u (cid:105) + u (cid:48) ) ⊗ ( (cid:104) u (cid:105) + u (cid:48) ) (cid:1)(cid:29) = ∇ · (cid:28) (cid:104) u (cid:105) ⊗ (cid:104) u (cid:105) + (cid:104) u (cid:105) ⊗ u (cid:48) + u (cid:48) ⊗ (cid:104) u (cid:105) + u (cid:48) ⊗ u (cid:48) (cid:29) = ∇ · ( (cid:104) u (cid:105) ⊗ (cid:104) u (cid:105) + (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) )= ( (cid:104) u (cid:105) · ∇ ) (cid:104) u (cid:105) + ∇ · (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) . Turbulence phenomena
Thus we obtain the incompressible RANS equations ∇ · (cid:104) u (cid:105) = 0 , (3.3.5) ∂ (cid:104) u (cid:105) ∂t + ( (cid:104) u (cid:105) · ∇ ) (cid:104) u (cid:105) = − ρ ∇(cid:104) p (cid:105) + ν ∆ (cid:104) u (cid:105) + (cid:104) f (cid:105) − ∇ · (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) . (3.3.6)Note that the Navier-Stokes equations and the RANS equations do not formally differ alot, except the additional term ∇ · (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) . This term is the divergence of the so calledReynolds stress tensor (RST) (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) = (cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105)(cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105)(cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105) (cid:104) u (cid:48) u (cid:48) (cid:105) . The application of the RANS equations only makes sense for turbulent flows, where itis assumed that highly fluctuating quantities occur. In the laminar case, u (cid:48) = 0 andtherefore the Navier-Stokes equations are reobtained. From Equation (3.3.6), we see thatthe RST affects the flow in form of additional stresses. It can be shown that this system ofequations is not closed. No matter how many manipulations we perform, there are alwaysmore unknowns than equations relating them. This is known as the closure problem ofturbulence and it arises because of the nonlinearity of the Navier-Stokes equations.The RST term needs to be modeled to circumvent this problem. In the field of RANSsimulation, the focus is primarily on modelling the RST and solving the RANS equations.A detailed overview is given in the upcoming chapter. The statistical moments defined in the previous subchapter are single point moments. Thatis, they contain only information about a variable at a point.In homogeneous turbulent flow, it only makes sense to have some statistical measure ofspatial information about the flow. For example, to draw conclusions about length scaleinformation, two point statistics are needed. The autocorrelation function is the correlationbetween velocity components at two different times. The normalized correlation function(autocorrelation function) is defined as D i ( τ ) = (cid:104) u i ( t ) u i ( t + τ ) (cid:105)(cid:104) u i ( t ) (cid:105) . (3.4.1)Note that D i (0) = 1. Also from Schwartz’s inequality, D i ( τ ) < τ (cid:54) = 0. Theautocorrelation coefficient is often used to define an integral scale of turbulence L i = (cid:90) ∞ D i ( τ ) dτ. (3.4.2)The integral length scale gives an estimate of the time interval, over which the velocitycomponent is correlated.The spatial correlation tensor gives the correlation between velocity components at two13 Turbulence phenomena different spatial locations and has an important interpretation in turbulent flows. It isdefined by R ( r ) = (cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105)(cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105)(cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105) (cid:104) u ( x ) u ( x + r ) (cid:105) . (3.4.3)To describe the various scales of spatial motion in a turbulent flow, it is more instructiveto work with the Fourier transform of the correlation tensor rather than the correlationtensor itself.We assume that the velocity field may be Fourier transformed under the certain require-ments. The Fourier transform of u i ( x, t ) isˆ u i ( k ) = (cid:90) ∞−∞ u i ( x ) e − jk · x dx, (3.4.4)herein k is the wave number and j the imaginary unit. The inverse Fourier transform is u i ( x ) = 1(2 π ) (cid:90) ∞−∞ ˆ u i ( k ) e jk · x dk. (3.4.5)The Fourier transformed spatial correlation tensor ˆ R is appropriately called the spectrumtensor or spectral density as it represents the contribution of a wave number k ˆ R ( k ) = (cid:90) ∞−∞ R ( r ) e − jk · r dr (3.4.6)= ˆ u ⊗ ˆ u ∗ , (3.4.7)where ∗ indicates the conjugate complex value.In other words, ˆ R ( k ) gives the wavenumber distribution of the correlation tensor. Eachwave number k corresponds to a physical space structure with a wavelength of πk .Of particular significance is the sum of the diagonal components of R ( r ) for r = 0. For thiscase we have tr ( R (0)) = u · u = 2 E, which is twice the kinetic energy E . The operator tr () indicates the trace of a squarematrix. In the spectral space, we have tr ( ˆ R ( k )) = ˆ u · ˆ u ∗ = 2 ˆ E ( k ) . And so E = 1(2 π ) (cid:90) ∞−∞ tr ( ˆ R ( k )) dk. (3.4.8)We are interested in the spectral kinetic energy in dependence of the magnitude of k , named k . Therefore, we integrate the spectrum tensor over a spherical shell with radius k ˆ R ( k ) = (cid:90) π (cid:90) π ˆ R ( k cos θ, k sin θ cos φ, k sin θ cos φ ) sin θ dθ dφ. (3.4.9)14 Turbulence phenomena
And its trace is tr ( ˆ R ( k )) = ˆ E ( k ) . We can rewrite Equation (3.4.8) to E = 1(2 π ) (cid:90) ∞ ˆ E ( k ) dk, (3.4.10)ˆ E ( k ) represents the contribution of the kinetic energy at a wavenumber of k , which is calledthe three dimensional energy spectrum.Integration of ˆ E ( k ) over all k gives the total kinetic energy. With an eddy of a particularsize, associated with a wavenumber of certain magnitude, the energy spectrum can beinterpreted to give the distribution of energy among the different eddy sizes. As discussedabove, the contribution of a wavenumber k corresponds to a structure with a wavelengthof πk . A large portion of the theoretical work on turbulent flows (including modeling) isconcerned with the description of energy in the wavenumber spectrum and the transfer ofenergy among the different wavenumbers and frequencies. As shortly described in the beginning of this chapter, another process in turbulent flow isthe breakdown of the large coherent structures into small eddies, while energy is transferredvia the cascadic breakdown of such vortices. In comparison to the more structured parentalvortices that sometimes seem to be ”predictable” or of periodic behavior, the little noisyeddies are completely three-dimensional, very random and tend to be completely indepen-dent of the big coherent structures. The observed scale distribution of these whirls matchquite well with the theory of Kolmogorov. He assumes that the energy spectrum can besplit into three sections: • The energy-containing, or also called integral scales. These are motions of permanentcharacter, mainly dominate the flow regime over many periods. More important,those scales contain by far the most of kinetic energy and are responsible of intro-ducing turbulent kinetic energy to the system. • In the second region, known as the inertial sub-range, the transitive scales are takingplace. These scales obey Kolmogorov’s famous law. They are completely independentof the forcing scale, not dominated by inertial forces rather than viscosity. Their mainaction is to transfer energy from the large scales to the very small ones. In this sectionthe energy spectrum only depends on k and (cid:15) , the dissipation rate. It can be derivedusing dimensional analysis ˆ E ( k ) = C K (cid:15) / k − / , (3.4.11)where C K is the Kolmogorov constant. • The dissipation region, which comprises the smallest scales, is the place where thekinetic energy is dissipated by the viscous effects. These are scales of motion which aresmaller than the Kolmogorov scale η , the length at which viscosity starts to strongly15 Turbulence phenomena damp the turbulent motion. The end of the curve is characterized by a rapid dropoffin energy content. The scale η is defined as η = (cid:18) ν (cid:15) (cid:19) . (3.4.12)Figure 3.1: Energy spectrum over k . The case of the fully developed turbulent plane channel flow is widely used in turbulencesimulation for validation. Due to its geometric simplicity, many advantages are includedas well as a lot of experimental and numerical data is available for comparison.We assume a fully developed turbulent flow in a channel with a large extension in x - and x -direction. And so, the flow is assumed to be statistically stationary (homogeneous) instreamwise and spanwise direction. Under those assumptions, the following regularitiesmay be deducted. In Figure 3.2 the schematic representation of the turbulent channel flowmay be seen. The bulk velocity and the corresponding bulk Reynolds number are definedas U b = 12 δ (cid:90) δ (cid:104) u (cid:105) dx , Re b = U b δν . (3.5.1)Due to symmetry in x , it follows that (cid:104) u (cid:105) = 0. By using the continuity equation andhomogeneity in x and x , it follows that ∂ (cid:104) u (cid:105) ∂x = 0. The x -momentum equation of theRANS equations simplifies to0 = − ρ d (cid:104) p (cid:105) dx + ν ∂ (cid:104) u (cid:105) ∂x − ∂ (cid:104) u (cid:48) u (cid:48) (cid:105) ∂x . With the relation ∂τ ges ∂x = d (cid:104) p (cid:105) dx , (3.5.2)16 Turbulence phenomena
Figure 3.2: Schematic representation of the turbulent channel flow.the total shear stress τ ges is obtained τ ges = µ ∂ (cid:104) u (cid:105) ∂x − ρ (cid:104) u (cid:48) u (cid:48) (cid:105) . (3.5.3)Due to the symmetry of the problem, the wall shear stress τ w is defined as τ ges ( x = 0) = µ ∂ (cid:104) u (cid:105) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 = − µ ∂ (cid:104) u (cid:105) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =2 δ = τ w . (3.5.4)Hence, the no-slip condition at the wall results in − d (cid:104) p (cid:105) dx = τ w δ , (3.5.5)giving τ ges as a linear function of x τ ges = τ w (1 − x δ ) . (3.5.6)The profiles of the molecular and turbulent shear stress are shown in Figure 3.3. In vicinityof the wall, the viscous diffusion is clearly dominating. The turbulent diffusion is zero atthe wall, but increases rapidly and dominates over almost the whole channel height.The turbulent velocity profile is determined by ν, τ w , ρ and δ . Due to the self similar be-havior of the mean velocity profile, the following Reynolds number independent parametersare introduced. The near-wall region is scaled by use of the friction velocity u τ and thewall unit l + u τ = (cid:114) τ w ρ , (3.5.7) l + = νu τ . (3.5.8)The corresponding friction Reynolds number Re τ is defined as Re τ = u τ δν = δl + . (3.5.9)17 Turbulence phenomena
Figure 3.3: Normalized viscous shear stress profile and turbulent shear stress profile.The non-dimensional distance to the wall y + is scaled with the wall unit l + y + = x l + = u τ ν x . (3.5.10)The non-dimensional mean velocity is scaled with the friction velocity u τ and reads u + = (cid:104) u (cid:105) u τ . (3.5.11)The law of the wall for turbulent flows is derived by assuming that turbulence near theboundary is a function only of the flow conditions pertaining at that wall and is independentof the flow conditions further away. The ansatz of so called universal velocity profile forturbulent flow near a wall is du + dy + = 1 y + Φ( y + ) , while Φ( y + ) is non-dimensional function of y + . In order to approximate this ordinarydifferential equation, it may be divided into several sections. • Viscous sublayer y + <
5: Since the turbulent fluctuations must go to zero at the wall,it follows that there is always a very small layer next to the wall, in which the flow isessentially laminar. Due to the no-slip condition at the wall, we get τ w = µ ∂ (cid:104) u (cid:105) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 ⇒ ρu τ = ρν u τ l + du + dy + (cid:12)(cid:12)(cid:12)(cid:12) y + =0 . For small y + it follows du + dy + ≈ ⇒ u + = y + . • Logarithmic layer y + > y/δ < .
3: Further from the wall, where the turbulentfluctuations dominate, the molecular diffusion may be neglected. In this region,18
Turbulence phenomena Φ( y + ) = const since it does not depend on ν and therefore y + anymore. There itholds du + dy + = 1 κy + .κ indicates the Von K´arm´an constant. Integration over y + leads to the logarithmiclaw of the wall u + = 1 κ ln ( y + ) + C log . For the channel flow, κ = 0 .
41 and C log = 5 was determined via experimental andnumerical tests. • Buffer layer 5 < y + <
30: In the buffer layer neither of the two laws hold.This universal velocity profile is sketched in Figure 3.4 and compared with numerical data.Figure 3.4: Comparison of DNS data ( Re τ = 395) with the universal velocity profile. The process of a laminar flow becoming turbulent is called laminar-turbulent transition. Itis mainly used in the context of boundary layers, but applies to any fluid flow. Generally,a laminar flow will transit to a turbulent flow if a certain limit of instabilities is exceededand amplifying magnitudes cannot be damped anymore. The instability characteristics ofa laminar flow mainly depend on its Reynolds number and on the intensity of mechanicalexcitation (disturbances in the flow, wall roughness, vibration,..). The interaction of allthe influences determine the degree of instability.The initial stage of the natural transition process is known as the receptivity phase. Thereceptivity gives a measure how the environmental disturbances are transformed to pertur-bations in the flow, deciding which form of transitional phase is taken.If the initially generated disturbances are small enough, the next stage of the laminar-turbulent transition process is that of primary mode growth. The growth rate can be19
Turbulence phenomena described by the linear stability theory assuming small amplitudes. In the case of bound-ary layers, the initial dominate instability mode will occur as a two-dimensional wave,called Tollmien-Schlichting wave, travelling in streamwise direction and its rotational axiscrosswise to the flow.The primary modes are for its part sensitive to disturbances, which leads to the secondphase of transition so called secondary instabilities. As the linear modes grow and begin toslightly distort the mean flow, they start to exhibit nonlinearities and the linear stabilitytheory no longer holds. The secondary instabilities often lead to coherent structures suchas harpin vortices. This stage is rapidly followed by tertiary instabilities and the finalbreakdown into fully developed turbulent regime.Figure 3.5: Illustration of laminar-turbulent transition process.
Solid boundaries interact with fluid flows by retarding motion tangential to the surfacevia viscous shear and by blocking the motion of fluid normal to the interface. Typically,near-wall structures are streaks of relatively small velocity, which are diving in regions ofhigher velocity. These coherent structures commonly appear in turbulent flows and arenamed low-speed streaks. Such vortex structure can be seen in Figure 3.6. The streaksare generated by the lifting of low-speed fluid near the wall induced from the streamwise20
Turbulence phenomena vortices. Until they are sufficiently lifted ( y + > (a) View in x , x -plane. Streamwise direction from left to right.(b) View in x , x -plane. Figure 3.6: Low-speed streaks.Despite of its prevalence, the shape of the turbulent energy spectrum shown previouslyis by far not universal. The transformation of bulk kinetic energy into turbulent energycan be mainly divided into two categories, both involving the existence of shear. The firstare free-shear flows such as mixing layers or jets and the second are wall-bounded flowssuch as boundary layers or channel flow. In free-shear flows, the growth of instabilities ismainly an inviscid process. These instabilities generally start to develop into big unsteadyquasi two-dimensional vortices and then degenerate into smaller eddies through the cas-cadic process described earlier. This implies that for equilibrium state, the small scales willgenerally tend to obey the law of Kolmogorov.Wall-bounded flows do not allow the inviscid growth of instabilities. In the existence ofsolid boundaries, the primary instabilities develop through viscous processes instead. Thismeans that turbulent energy is introduced into the system at small scales near the wall. Ex-periments have shown that the near-wall eddies are just as anisotropic and inhomogeneousas the large whirls in free-shear flows. This means that there must be some mechanism totransfer energy from small to large scales away from the wall, which is called the inverseenergy cascade.In the near-wall boundary layer, viscosity acts as a momentum sink to the core flow, in21
Turbulence phenomena a similar way as the dissipative effect of the small scale end of the kinetic energy spec-trum. The mean flow momentum and therefore mean kinetic energy is transfered to thesurface layer by Reynolds stresses and there converted into turbulent kinetic energy andheat (through viscous dissipation). Most of the turbulent energy generated near the wall islost to dissipation because of large velocity gradients in this region. A significant portionis however transported to the outer flow through turbulent diffusion before it dissipates.Since turbulent production in the outer flow is generally small, this makes the surface layerthe main source of turbulent energy of the entire flow.In order to obtain a clearer picture of the balance of the turbulent kinetic energy, considerits transport equation. The turbulent kinetic energy K is defined as K = 12 (cid:104) u (cid:48) · u (cid:48) (cid:105) . (3.6.1)Hence, the transport equation of K may derived in the way that firstly substituting theReynolds decomposition into the Navier-Stokes equations (Equation (2.4.2)), then sub-tracting it by the RANS equations (Equation (3.3.6)) and multiply with u (cid:48) .The following equation is obtained u (cid:48) · (cid:18) ∂u (cid:48) ∂t + ∇ · (cid:0) (cid:104) u (cid:105) ⊗ u (cid:48) + u (cid:48) ⊗ (cid:104) u (cid:105) + u (cid:48) ⊗ u (cid:48) − (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) (cid:1)(cid:19) = u (cid:48) · (cid:18) − ρ ∇ p (cid:48) + ν ∆ u (cid:48) (cid:19) . This equation may be further simplified with the relation u (cid:48) · ∂u (cid:48) ∂t = ∂ ( u (cid:48) · u (cid:48) ) ∂t and theincompressibility constraint to u (cid:48) · (cid:18) ∇ · (cid:0) (cid:104) u (cid:105) ⊗ u (cid:48) + u (cid:48) ⊗ (cid:104) u (cid:105) (cid:1)(cid:19) = u (cid:48) · (cid:0) ∇ · ( (cid:104) u (cid:105) ⊗ u (cid:48) ) (cid:1) + u (cid:48) · (cid:0) ∇ · ( u (cid:48) ⊗ (cid:104) u (cid:105) ) (cid:1) = u (cid:48) · ( (cid:104) u (cid:105) · ∇ ) u (cid:48) + u (cid:48) · ( u (cid:48) · ∇ ) (cid:104) u (cid:105) = 12 (cid:104) u (cid:105) · ∇ ( u (cid:48) · u (cid:48) ) + ( u (cid:48) ⊗ u (cid:48) ) : ∇(cid:104) u (cid:105) , while : denotes the Frobenius inner product of two matrices.This leads to the equation12 ∂ ( u (cid:48) · u (cid:48) ) ∂t + 12 (cid:104) u (cid:105) · ∇ ( u (cid:48) · u (cid:48) ) + ( u (cid:48) ⊗ u (cid:48) ) : ∇(cid:104) u (cid:105) + 12 ∇ · (cid:0) ( u (cid:48) · u (cid:48) ) u (cid:48) (cid:1) − u (cid:48) · (cid:0) ∇ · (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) (cid:1) = − ρ ∇ · ( u (cid:48) p (cid:48) ) + νu (cid:48) · ∆ u (cid:48) . Reynolds averaging over the whole equation leads to ∂K∂t + (cid:104) u (cid:105) · ∇ K = −(cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) : ∇(cid:104) u (cid:105) − ∇ · (cid:104) ( u (cid:48) · u (cid:48) ) u (cid:48) (cid:105) − ρ ∇ · (cid:104) u (cid:48) p (cid:48) (cid:105) + ν (cid:104) u (cid:48) · ∆ u (cid:48) (cid:105) . The term (cid:104) u (cid:48) · ∆ u (cid:48) (cid:105) may be further reduced to12 (cid:104) ∆( u (cid:48) · u (cid:48) ) (cid:105) = 12 (cid:10) ∇ · (cid:0) ∇ ( u (cid:48) · u (cid:48) ) (cid:1)(cid:11) = (cid:10) ∇ · (cid:0) ( ∇ u (cid:48) ) u (cid:48) (cid:1) (cid:105) = (cid:104) u (cid:48) · ∆ u (cid:48) (cid:105) + (cid:104)∇ u (cid:48) : ∇ u (cid:48) (cid:105) . Turbulence phenomena
Finally, the turbulent kinetic energy equation is ∂K∂t + (cid:104) u (cid:105) · ∇ K = Π − (cid:15) + ν ∆ K + ∇ · Γ . (3.6.2)The physically interpretation of each single term of Equation (3.6.2) is: • Material derivative of K : DKDt = ∂K∂t + (cid:104) u (cid:105) · ∇ K. • Production of K (generally Π > −(cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) : ∇(cid:104) u (cid:105) . • Molecular dissipation ( (cid:15) ≥ (cid:15) = ν (cid:104)∇ u (cid:48) : ∇ u (cid:48) (cid:105) . • Viscous diffusion: ν ∆ K. • Redistribution terms: Sum of pressure-diffusion and turbulent transport term ∇ ·
Γ = ∇ · (cid:0) − ρ (cid:104) u (cid:48) p (cid:48) (cid:105) − (cid:104) ( u (cid:48) · u (cid:48) ) u (cid:48) (cid:105) (cid:1) . If the flow is parallel and in equilibrium (e.g. channel flow), the time-averaged streamwisederivatives and cross-stream velocities tend to zero. Assuming the contribution of viscousdiffusion is small compared to the turbulent contribution, Equation (3.6.2) reduces to ∂ Γ ∂x = Π − (cid:15), so that the only significant mean flux is the wall-normal component Γ . Figure 3.7 showsthe balance between Π and (cid:15) produced by numerical data from [MKM99] and Re τ = 395.It clearly shows that net turbulent production occurs only in the outer viscous and bufferregion, while dissipation exceeds production in the viscous sublayer. This interpretationis reinforced by the energy flux, which is positive through most of the domain (exceptvery near the wall), signifying a movement of energy away from the surface. Hardly anyturbulence is produced in the central part of the channel. Hence, most of the turbulentenergy is provided by the flux away from the walls.Since most of the turbulent energy is contained in the largest scale eddies and the eddysize is limited by the distance from the wall, the near-wall energy containing eddies willbe small while those centered in the core will be large. Therefore, the transfer of energyfrom the small to large scales is an example of an inverse energy cascade as opposed to theclassical Kolmogorov cascade. Since the nature of the fluxes that are being transferred aredifferent from that in the Kolmogorov theory, the spectral slope should be different too.In the work of [DV06], dimensional considerations have brought that the inverse energyspectrum has the following shape ˆ E ( k ) ∼ k − . (3.6.3)23 Turbulence phenomena (a) (b)
Figure 3.7: (a): Normalized energy balance Π − (cid:15) . (b): Normalized energy flux Γ .24 Simulation principles and modelling
This chapter gives a brief overview on the several principles to numerically simulate tur-bulence. A CFD simulation generally consists of many steps, from modelling to post-processing. One decisive step is the translation from a mathematical model to an algebraicsystem of equations, so called discretization. The discretization methods used in this thesisare given in more detail in Chapter 5 and 6. The most well-known principles of simulatingincompressible turbulent flows will be discussed in the upcoming subsections.Basically, three different approaches are considered: • Direct numerical simulation (DNS): Resolving all relevant turbulent scales. A hugecomputational effort is necessary therefore. • RANS simulation: Modelling all scales. Hence, the computational cost is significantlylow. • Large eddy simulation (LES): Resolve the large turbulent scales and modelling thesmall scales.A comparison of the velocity signal of the different approaches is shown in Figure 4.1.Figure 4.1: Comparison of velocity signal u of different turbulence simulation principlesover time t . The DNS of a turbulent fluid flow directly solves the Navier-Stokes equations of the formof Equation (2.4.1) and (2.4.2) without any modelling concept or assumptions. To obtaina physically correct flow, all relevant scales from the Kolmogorov scale up to the integralscale have to be fully numerically resolved. For example, in atmospheric flows the scale25
Simulation principles and modelling range includes more than nine orders of magnitude. Hence, particular requirements of theefficiency of solving algorithms have to be made. In most of the engineering applications,no such temporal and spatial richness of detail is desired. Therefore, DNS is more or lessused for validation and foundational research.An estimation of the required resolution of computational domain may be given, assumingthe mesh size of h ≈ η . Then, the number of grid points in one dimension is given by theratio of biggest to smallest length scale N D ≈ Lη . Through dimensional analysis it follows N D ≈ Lη = L ( ν (cid:15) ) / = L (cid:0) U b L (cid:1) / ν / = (cid:18) U b Lν (cid:19) / = Re / b , while the dissipation rate is (cid:15) = U b L .One similar ansatz is made for the time scale N T ≈ TT η = LU b (cid:18) U b Lν (cid:19) / = (cid:114) U b Lν = (cid:112) Re b . Furthermore, the total computational cost may scale with N T N D ≈ Re / b . Since it grows almost cubic with respect to the Reynolds number, high turbulent flow incomplex geometries are computationally unfeasible in the near future. Herein, it has to bementioned that it only gives a rough estimation, since the total complexity depends on avariety of things (solving algorithm, discretization,..). Additionally, in order to reduce thecomputational effort, for several cases a full resolution of the flow domain is often not ofmajor importance even for a DNS.
As the DNS of turbulent flows is still prohibitively expensive, models based on the RANSequations are very commonly employed in CFD nowadays. For engineering purposes, thedetailed resolved scales are not desired and often only statistically values are of main inter-est. In RANS simulation, the whole range of turbulent scales is modeled by decomposingand averaging the governing equations of fluid motion, leading to the previous derivedRANS equations of the from of Equation (3.3.5) and (3.3.6).Hence, its main goal is the numerical solution of its equations. The aim of the statisticallyturbulence modelling is to solve the closure problem of Equation (3.3.6). In order to obtaina physically rightful model, it has to fulfill several physical and mathematical constraints(e.g. tensorial consistency) and development of such models is still a field of research.It may be divided into two classes, one are the eddy-viscosity models and the other areReynolds-stress models. In this thesis, only models of the first mentioned class will bediscussed.
Boussinesq [Sch07] firstly introduced the concept of eddy viscosity. It is based on theanalogy of turbulent to gas kinetic processes and he proposed a relation of the turbulence26
Simulation principles and modelling stresses to the mean flow. The transfer of momentum in between molecules results in fric-tion. In macroscopic scales, the colliding of ”turbulence clusters” somehow emerges tensionsin rather similar way as in the microscopic scale. The principle component of modelling isto determine a velocity scale, which is characteristic for the intensity of turbulent mixtureand an appropriate typically length scale, where turbulence takes place.Equation (3.3.6) may be rewritten in the form ∂ (cid:104) u (cid:105) ∂t + ( (cid:104) u (cid:105) · ∇ ) (cid:104) u (cid:105) = − ρ ∇(cid:104) p (cid:105) + ∇ · (cid:0) ν (cid:104) S (cid:105) − (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) (cid:1) , neglecting the volume force term and introducing the mean strain rate tensor (cid:104) S (cid:105) = 12 ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) . According to the linear proportionality of the molecular diffusion term to the viscosity andmean strain rate tensor, the same proportionality is assumed for the Reynolds stresses.The ansatz is (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) − KI = − ν T (cid:104) S (cid:105) , (4.2.1)where ν T ∈ C (Ω × [0 , T ] , R ), ν T > ν T ∼ l T U T . In incompressible flows, the trace of the strain rate tensor is zero (divergence-free velocity) and therefore the viscous stresses are only described via the deviatoric partof the tensor. The left-hand side of Equation (4.2.1) represents the deviatoric part of theRST (its trace is generally not zero), in order to obtain tensorial consistency.Find ν T = ν T ( x, t ) and Equation (4.2.1) closes the RANS equations of form of Equa-tion (3.3.5) and (3.3.6). Problem 2 ( Reynolds-averaged Navier-Stokes equations ) . Let ν, ρ ∈ R \ { } , ν T ∈ C (Ω × [0 , T ] , R ) , (cid:104) u (cid:105) D ∈ C ( ∂ Ω D , R ) , g ∈ C ( ∂ Ω N , R ) and (cid:104) u (cid:105) ∈ C (Ω , R ) find (cid:104) u (cid:105) ∈ C (Ω × [0 , T ] , R ) and ˆ p ∈ C (Ω × [0 , T ] , R ) such that ∇ · (cid:104) u (cid:105) = 0 in Ω × [0 , T ] , (4.2.2) ∂ (cid:104) u (cid:105) ∂t + ( (cid:104) u (cid:105) · ∇ ) (cid:104) u (cid:105) = − ρ ∇ ˆ p + ∇ · (cid:0) ν + ν T ) (cid:104) S (cid:105) (cid:1) in Ω × [0 , T ] , (4.2.3) (cid:104) u (cid:105) = (cid:104) u (cid:105) in Ω , t = 0 , (cid:104) u (cid:105) = (cid:104) u (cid:105) D in ∂ Ω D × [0 , T ] ,∂ (cid:104) u (cid:105) ∂n = g in ∂ Ω N × [0 , T ] . Since the pressure has no thermodynamic significance in incompressible flows, the isotropicpart of the RST is incoporated by the modified pressure ˆ p = (cid:104) p (cid:105) + ρK .The weakness of the Boussinesq assumption is that it is not valid in general. There is noevidence that the Reynolds stress tensor must be proportional to the strain rate tensor.In most of the simple flow cases, the linear proportionality gives a good approximation.For flow scenarios with strong curvature or strongly accelerated or decelerated flows, theBoussinesq assumption is simply not valid. 27 Simulation principles and modelling
Contrary to the molecular viscosity, the turbulent viscosity is no material property anddepends on the flow itself. In order to be able to determine ν T , further equations haveto be provided. Nowadays, a large amount of different turbulence models exist, eachwith its advantages and drawbacks. In this thesis, the three most important two-equationturbulence models for wall-bounded flows are explained. K − (cid:15) model Probably the most well-known turbulence model, the K − (cid:15) model and its variations werefirstly introduced by Harlow and Nakayama [HN68]. Through dimensional analysis, thefollowing proportionalities for the characteristic velocity and length scale are given l T ∼ K / (cid:15) ,U T ∼ √ K. Core assumption of this model is the following ansatz for the eddy viscosity ν T = C µ K (cid:15) , (4.2.4)where C µ is a constant.To be able to compute the turbulent viscosity as in Equation (4.2.4), two additional trans-port equations for K and (cid:15) need to be deducted. For the kinetic turbulent energy, Equa-tion (3.6.2) is used and its redistribution term is modelled asΓ ≈ ν T σ K ∇ K, with σ K to be a constant. This new modeled term should act as a turbulent diffusion. Theproduction term Π may be rewritten and closed with Equation (4.2.1)Π = −(cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) : ∇(cid:104) u (cid:105) = − ( (cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) − KI ) : (cid:104) S (cid:105) = 2 ν T (cid:104) S (cid:105) : (cid:104) S (cid:105) . Secondly, an equation for (cid:15) has to be determined. Although a transport equation maybe derived for (cid:15) as well, it contains even more unclosed terms and is subtantially morecomplicated than the transport equation for K . Therefore, analogous to the K transportequation, a formula for (cid:15) is specified in the following ∂(cid:15)∂t + (cid:104) u (cid:105) · ∇ (cid:15) = C (cid:15) Π (cid:15)K − C (cid:15) (cid:15) K + ∇ · (cid:0) ( ν + ν T σ (cid:15) ) ∇ (cid:15) (cid:1) , where C (cid:15) , C (cid:15) and σ (cid:15) are model constants.The coupled system of equations is summarized in the following problem.28 Simulation principles and modelling
Problem 3 ( K − (cid:15) transport equations ) . Let σ K , σ (cid:15) , C (cid:15) , C (cid:15) , ν ∈ R \ { } , ν T ∈ C (Ω × [0 , T ] , R ) , (cid:104) u (cid:105) ∈ C (Ω × [0 , T ] , R ) , K D ∈ C ( ∂ Ω D , R ) , g k ∈ C ( ∂ Ω N , R ) , (cid:15) D ∈ C ( ∂ Ω D , R ) , g (cid:15) ∈ C ( ∂ Ω N , R ) , K ∈ C (Ω , R ) and (cid:15) ∈ C (Ω , R ) find K ∈ C (Ω × [0 , T ] , R ) and (cid:15) ∈ C (Ω × [0 , T ] , R ) such that ∂K∂t + (cid:104) u (cid:105) · ∇ K = Π − (cid:15) + ∇ · (cid:0) ( ν + ν T σ K ) ∇ K (cid:1) in Ω × [0 , T ] , (4.2.5) ∂(cid:15)∂t + (cid:104) u (cid:105) · ∇ (cid:15) = C (cid:15) Π (cid:15)K − C (cid:15) (cid:15) K + ∇ · (cid:0) ( ν + ν T σ (cid:15) ) ∇ (cid:15) (cid:1) in Ω × [0 , T ] , (4.2.6)Π = 12 ν T ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) : ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) in Ω × [0 , T ] ,K = K in Ω , t = 0 ,(cid:15) = (cid:15) in Ω , t = 0 ,K = K D in ∂ Ω D × [0 , T ] ,(cid:15) = (cid:15) D in ∂ Ω D × [0 , T ] ,n · ∇ K = g k in ∂ Ω N × [0 , T ] ,n · ∇ (cid:15) = g (cid:15) in ∂ Ω N × [0 , T ] . K − ω model The K − (cid:15) model suffers from some major drawbacks, one of them is the lack of sensitivityto adverse pressure gradients as it is observed that under such conditions the model tendsto overestimate the shear stress and by that delay flow separation.The K − ω model as devised by Wilcox [Wil94] overcomes this weakness of the K − (cid:15) modelby the definition of the specific dissipation rate ωω = (cid:15)C µ K . (4.2.7)The eddy-viscosity reads ν T = Kω . (4.2.8)Again, the transport equation for ω is derived in a very similar way as for Equation (4.2.6).29 Simulation principles and modelling
Problem 4 ( K − ω transport equations ) . Let β ∗ , σ K , σ ω , C ω , C ω , ν ∈ R \ { } , ν T ∈ C (Ω × [0 , T ] , R ) , (cid:104) u (cid:105) ∈ C (Ω × [0 , T ] , R ) , K D ∈ C ( ∂ Ω D , R ) , g k ∈ C ( ∂ Ω N , R ) , ω D ∈ C ( ∂ Ω D , R ) , g ω ∈ C ( ∂ Ω N , R ) , K ∈ C (Ω , R ) and ω ∈ C (Ω , R ) find K ∈ C (Ω × [0 , T ] , R ) and ω ∈ C (Ω × [0 , T ] , R ) such that ∂K∂t + (cid:104) u (cid:105) · ∇ K = Π − β ∗ Kω + ∇ · (cid:0) ( ν + ν T σ K ) ∇ K (cid:1) in Ω × [0 , T ] , (4.2.9) ∂ω∂t + (cid:104) u (cid:105) · ∇ ω = C ω ωK Π − C ω ω + ∇ · (cid:0) ( ν + ν T σ ω ) ∇ ω (cid:1) in Ω × [0 , T ] , (4.2.10)Π = 12 ν T ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) : ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) in Ω × [0 , T ] ,K = K in Ω , t = 0 ,ω = ω in Ω , t = 0 ,K = K D in ∂ Ω D × [0 , T ] ,ω = ω D in ∂ Ω D × [0 , T ] ,n · ∇ K = g k in ∂ Ω N × [0 , T ] ,n · ∇ ω = g ω in ∂ Ω N × [0 , T ] . This model has superior performance for wall-bounded flows with relatively low Reynoldsnumbers. K − ω SST model
The shear stress transport (SST) K − ω model, firstly revised by Menter [Men94], combinesthe two previous models such that the K − ω model is used in the proximity of walls andswitches to the K − (cid:15) in the free shear region. Menter recognized that the (cid:15) transportequation from Equation (4.2.6) may be transformed into a new ω transport equation viasubstituting ω from Equation (4.2.7). This new transformed equation looks very similarto the one from Equation (4.2.10), expect it includes an additional non-conservative cross-diffusion term 2 C ω ω ∇ K : ∇ ω, introducing a new constant C ω .The inclusion of this term will potentially make it able to smoothly switch between themodels via a blending functions (highly nonlinear functions of the several parameter, in-cluding the wall distance). The additional functions are given by F = tanh( ζ ) , (4.2.11) ζ = (cid:18) min (cid:20) max (cid:0) √ Kβ ∗ ωd , νωd (cid:1) , C ω KCDd (cid:21)(cid:19) , (4.2.12) CD = max (cid:0) C ω ω ∇ K : ∇ ω, − (cid:1) , (4.2.13) F = tanh( ζ ) , (4.2.14) ζ = (cid:18) max (cid:0) √ Kβ ∗ ωd , νωd (cid:1)(cid:19) . (4.2.15)30 Simulation principles and modelling
The variable d describes the distance to the nearest wall. If F = 1 the cross-diffusionterm disappears and the K − ω model is recovered, otherwise if F = 0 the K − (cid:15) model isreobtained.The model parameters (e.g. σ K , σ ω ) are blended via χ = F χ + (1 − F χ , (4.2.16)where χ represents the K − ω model constants and χ the ones from the K − (cid:15) model.The eddy viscosity is computed from ν T = a Kmax ( a ω, SF ) , (4.2.17)while a is a model constant and S is the magnitude of the strain rate tensor S = (cid:113) (cid:104) S (cid:105) : (cid:104) S (cid:105) . (4.2.18) Problem 5 ( K − ω SST transport equations ) . Let β ∗ , σ K , σ ω , C ω , C ω , C ω , ν ∈ R \ { } , ν T ∈ C (Ω × [0 , T ] , R ) , (cid:104) u (cid:105) ∈ C (Ω × [0 , T ] , R ) , K D ∈ C ( ∂ Ω D , R ) , g k ∈ C ( ∂ Ω N , R ) , ω D ∈ C ( ∂ Ω D , R ) , g ω ∈ C ( ∂ Ω N , R ) , K ∈ C (Ω , R ) and ω ∈ C (Ω , R ) find K ∈ C (Ω × [0 , T ] , R ) and ω ∈ C (Ω × [0 , T ] , R ) such that ∂K∂t + (cid:104) u (cid:105) · ∇ K = Π − β ∗ Kω + ∇ · (cid:0) ( ν + ν T σ K ) ∇ K (cid:1) in Ω × [0 , T ] , (4.2.19) ∂ω∂t + (cid:104) u (cid:105) · ∇ ω = C ω ν T Π − C ω ω + ∇ · (cid:0) ( ν + ν T σ ω ) ∇ ω (cid:1) + 2(1 − F ) C ω ω ∇ K : ∇ ω in Ω × [0 , T ] , (4.2.20)Π = 12 ν T ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) : ( ∇(cid:104) u (cid:105) + ∇(cid:104) u (cid:105) T ) in Ω × [0 , T ] ,K = K in Ω , t = 0 ,ω = ω in Ω , t = 0 ,K = K D in ∂ Ω D × [0 , T ] ,ω = ω D in ∂ Ω D × [0 , T ] ,n · ∇ K = g k in ∂ Ω N × [0 , T ] ,n · ∇ ω = g ω in ∂ Ω N × [0 , T ] . The SST model exhibits the closest agreement with the test cases and works very wellin wall-bounded shear flows with adverse pressure gradients. Although, its computationalcost is slightly higher than the most linear eddy-viscosity models (due to the nonlinearitiesin its evaluation), the effort is reasonable in comparison with a LES.
These two-equation eddy-viscosity models mainly describe two coupled scalar convectiondiffusion problems with (nonlinear) source and sink terms. The scalar quantity φ ∈ C (Ω × Simulation principles and modelling [0 , T ] , R ) gets transported along a given stream (cid:104) u (cid:105) with the property ∇ · (cid:104) u (cid:105) = 0. Thecorresponding convective term is then simply ( (cid:104) u (cid:105) · ∇ φ ). Convection is superimposed bydiffusion which describes the exchange of the transported quantity on a molecular andturbulent level. The corresponding diffusive flux term is ∇ · (cid:0) ( ν + ν T ) ∇ φ (cid:1) . The term f describes the sources and sinks and may consist of nonlinear terms.Generally, the scalar convection diffusion problem reads: Problem 6 ( The scalar convection diffusion equation ) . Let ν ∈ R \ { } , ν T ∈ C (Ω × [0 , T ] , R ) , (cid:104) u (cid:105) ∈ C (Ω × [0 , T ] , R ) , φ D ∈ C ( ∂ Ω D , R ) , g φ ∈ C ( ∂ Ω N , R ) and φ ∈ C (Ω , R ) find φ ∈ C (Ω × [0 , T ] , R ) such that ∂φ∂t + (cid:104) u (cid:105) · ∇ φ = ∇ · (cid:0) ( ν + ν T ) ∇ φ (cid:1) + f in Ω × [0 , T ] , (4.2.21) ∇ · (cid:104) u (cid:105) = 0 in Ω × [0 , T ] ,φ = φ in Ω , t = 0 ,φ = φ D in ∂ Ω D × [0 , T ] ,n · ∇ φ = g φ in ∂ Ω N × [0 , T ] . In the proximity of walls, most of the eddy-viscosity models suffer in performance, sincethe core assumptions of isotropic turbulence and high Reynolds number flow do not hold inthis area. Especially the K − (cid:15) model is not valid close to solid walls. This situation givesrise to a plethora of near-wall conditions. Generally, two approaches are distinguished, thelow Reynolds number (LRN) treatment and high Reynolds number (HRN) treatment.Models using the LRN version integrate every equation up to the wall using the appropriateDirichlet or Neumann boundary conditions for the physical quantities at the solid wall.Therefore, the first computational cell has to be in y + ∼
1, resulting in fine resolved meshesclose to walls. Additionally, some models use damping functions of the model parametersto guarantee asymptotic consistency with the turbulent boundary layer behavior.The HRN approach uses wall functions, which are applied in the nearest cell at the wallinstead of integration. These functions often rely on approximated log-law velocity profiles.No direct boundary conditions are therefore necessary, since the closest cell to the wall iscomputed according to the wall functions. This method enhances the computational effort,however it is not suitable for more complex scenarios.In this thesis, we stick to the LRN approach for all RANS simulations.
The LES is a technique intermediate between the DNS and the solution of RANS simulation.In LES, the contribution of the largest, kinetic energy-carrying eddies is computed exactly,while only the effect of the smaller scale whirls is modeled. Since the smaller structures tendto be more homogeneous and universal and less affected by the global conditions than thelarger scaled eddies, it seems to be more attractive to model the small scale part. Similarto DNS, LES provides a three-dimensional, time-dependent solution of the Navier-Stokes32
Simulation principles and modelling equations.To seperate the large scales of motion from the smaller ones, some kind of averaging has tobe done. In contrast to the ensemble average in RANS, in LES the averaging operator isgenerally a spatial and temporal low-pass filter. Formally, each flow variable is decomposedin the large and small scale part u ( x, t ) = u ( x, t ) + u (cid:48)(cid:48) ( x, t ) , p ( x, t ) = p ( x, t ) + p (cid:48)(cid:48) ( x, t ) . (4.3.1)The overbar donates the resolved, large components (grid scales (GS)) and the doubleprime the unresolved part (sub-grid scales (SGS)).To be able to extract the low frequency components of the quantity, the filtering operationis defined as φ ( x, t ) = (cid:90) Ω G ( x − x (cid:48) , t − t (cid:48) ; ∆ δ ) φ ( x (cid:48) , t (cid:48) ) dx (cid:48) dt (cid:48) , (4.3.2)where G is the filter convolution kernel, associated with the cutoff length scale ∆ δ alsocalled filter width. An important note is that the LES filtering operation does not satisfythe properties of a Reynolds operator as described in Definition 3. A LES filter must satisfythe following set of properties. Definition 4 ( LES filter operator ) . Let φ, ψ ∈ C (Ω × [0 , T ] , R ) , c ∈ R and G theLES filter operator as defined in Definition 4.3.2, then the following properties have tobe satisfied: • Linearity: φ + ψ = φ + ψ • Commutation with derivatives: ∂φ∂x = ∂φ∂x ∂φ∂t = ∂φ∂t • Constants: c = c which implies that, (cid:90) Ω G ( x (cid:48) ; ∆ δ ) dx (cid:48) = 1 Generally, this filter operator does not satisfy the following properties: • φ (cid:54) = φ • φ (cid:48)(cid:48) (cid:54) = 0 Meaning that it is not a Reynolds operator.
The Fourier transform of the filter operator is defined asˆ G ( k ; ∆ δ ) = (cid:90) ∞−∞ G ( x ; ∆ δ ) e − jk · x dx. (4.3.3)In frequency domain, the filtering operation is simply obtained byˆ φ ( k ) = ˆ G ( k ; ∆ δ ) ˆ φ ( k ) . (4.3.4)33 Simulation principles and modelling (a)(b)
Figure 4.2: Comparison unfiltered and filtered (top-hat, sharp spectral and Gaussian filterkernel) velocity signal with filter width ∆ δ = . (a): Time domain. (b):Frequency domain. 34 Simulation principles and modelling
The most common filter kernels that have been applied to LES are the Gaussian filter,sharp spectral filter or top-hat filter. The one-dimensional filter kernels are defined as:Filter Function Fourier transformGaussian G G ( x ; ∆ δ ) = (cid:113) π ∆ δ e − x / ∆ δ ˆ G G ( k ; ∆ δ ) = e − k ∆ δ / Sharp spectral G S ( x ; ∆ δ ) = sin xπ ∆ δ xπ ˆ G S ( k ; ∆ δ ) = (cid:40) , k ≤ π ∆ δ , elseTop-hat G T ( x ; ∆ δ ) = (cid:40) δ , − ∆ δ ≤ x ≤ ∆ δ , else ˆ G T ( k ; ∆ δ ) = sin k ∆ δ k ∆ δ If this filtering operator is applied to the Navier-Stokes equations of the form of Equa-tion (2.4.1) and (2.4.2) (again neglecting the volume force term), we obtain the spatial-filtered Navier-Stokes equations (SFNS) ∇ · u = 0 , (4.3.5) ∂u∂t + ∇ · ( u ⊗ u ) = − ρ ∇ p + ∇ · ν ( ∇ u + ∇ u T ) . (4.3.6)The effect of the spatial filtering can be also seen in a damped energy spectrum in Figure 4.3.Figure 4.3: Comparison unfiltered and filtered energy spectrum.Although the definition of the quantities differs from that in the RANS equation, the closureproblem is conceptually very similar. Since u ⊗ u (cid:54) = u ⊗ u , a model approximation has tobe taken into account. The common way to address this issue of closure is to introduce theso called SGS stress tensor, T = u ⊗ u − u ⊗ u. (4.3.7)The symmetric tensor T has to have the property that | T | → δ →
0, so that in thelimit of mesh spacing the DNS solution is recovered. The SGS stress tensor is functionallyvery similar to the RST, but the physics of the problem is somewhat different.35
Simulation principles and modelling
Inserting T in Equation (4.3.6), it follows ∂u∂t + ( u · ∇ ) u = − ρ ∇ p + ∇ · ν ( ∇ u + ∇ u T ) − ∇ · T . (4.3.8)The filtered kinetic energy E can be decomposed into E = 12 u · u = 12 ( u · u + u · u − u · u )= E GS + 12 tr ( T ) , while E GS is the kinetic energy of the resolved filtered scales and tr ( T ) the SGS energy.The conservation equation for E GS can be obtained by multiplying the filtered momentumtransport equation (Equation (4.3.8)) by the filtered velocity u to yield u · (cid:18) ∂u∂t + ( u · ∇ ) u (cid:19) = u · (cid:18) − ρ ∇ p + ν ∆ u − ∇ · T (cid:19) . The derivation is quite similar to the one in Chapter (3.6.2). Further rearrangements leadto u · ( ν ∆ u ) = 12 ν ∆( u · u ) − ν ∇ u : ∇ u = ν ∆ E GS − (cid:15) GS , and u · ( ∇ · T ) = ∇ · ( T u ) − T : ∇ u. The final transport equation is then ∂E GS ∂t + ( u · ∇ ) E GS = ∇ · (cid:0) − ρ p u − T u (cid:1) + ν ∆ E GS − (cid:15) GS + T : ∇ u. (4.3.9)The pyhisical interpretation of the first term on the right-handside of Equation (4.3.9) isthe redistribution term and the second one the viscous diffusion. The term (cid:15) GS > (cid:15) SGS = T : ∇ u. The SGS dissipation may be positive or negative, meaning if it is negative, energy dissipatesfrom the resolved scales to the sub-grid scales, which is called forwardscatter. However,if it is positive, energy transfers from the sub-grid scales to the resolved ones, so calledbackscatter.A much smaller part of the turbulent energy spectrum is covered by the SGS energy thanthe filtered kinetic energy, meaning that accuracy of the SGS model may be less crucialthan in RANS.In the following, a method to successfully model the SGS stress tensor is shown.36
Simulation principles and modelling
The introduced model in this subsection have many parallels with the RANS counterparts.Nevertheless, the fact that a much smaller part of the turbulent energy spectrum has tobe modeled, contributes to a smaller error potential and simple models may produce goodresults.In LES, the dissipative scales are generally not resolved, therefore the main role of the SGSmodel is to extract energy from the resolved scales. This can be accomplished with aneddy-viscosity model similar to the RANS model. To this end we assume T − tr ( T ) I = − ν SGS ( ∇ u + ∇ u T ) = − ν SGS S, (4.3.10)where S is the filtered strain rate tensor S = 12 ( ∇ u + ∇ u T ) . (4.3.11)The isotropic part of the SGS stress tensor is incoporated by the filtered pressure. TheSFNS problem then reads in the following: Problem 7 ( Spatial-filtered Navier-Stokes equations ) . Let ν , ρ ∈ R \ { } , ν SGS ∈ C (Ω × [0 , T ] , R ) , u D ∈ C ( ∂ Ω D , R ) , g ∈ C ( ∂ Ω N , R ) and u ∈ C (Ω , R ) find u ∈ C (Ω × [0 , T ] , R ) and ˆ p ∈ C (Ω × [0 , T ] , R ) such that ∇ · u = 0 in Ω × [0 , T ] , (4.3.12) ∂u∂t + ( u · ∇ ) u = − ρ ∇ ˆ p + ∇ · (cid:0) ( ν + ν SGS )( ∇ u + ∇ u T ) (cid:1) in Ω × [0 , T ] , (4.3.13) u = u in Ω , t = 0 ,u = u D in ∂ Ω D × [0 , T ] ,∂u∂n = g in ∂ Ω N × [0 , T ] . Two fundamental topics in LES have to be mentioned.Firstly, there are mainly three different approaches, implicit, implicitly filtered and ex-plicitly filtered LES. The common method is implicit LES, where the system of equationsare never acted upon a filtering. The filtering is provided implicitly by two causes, thecomputational grid and the discretization used. Traditionally, the filtering is done by thegrid itself and the inherent numerical diffusion acts implicitly as sub-grid scale model. Itis also known as quasi or coarse DNS. In implicitly filtered LES, the filtered Navier-Stokesequations given in Problem 7 with an appropriate SGS model are solved numerically. For-mally, the governing system of equations do not differ from the unsteady RANS equations.Nonetheless, an explicit filter can still be applied in order to derive variables used in theSGS model. In explicitly filtered LES, the application of an numerical filter on the equa-tions is performed at each time step. Therefore, it is possible to control the shape and typeof the filter. In this thesis, we investigate the first two procedures.Secondly, the derivation of the SFNS equations takes the advantage that the filter operatorcommutes with differentiation, which holds in absence of boundaries. Briefly, in presence37
Simulation principles and modelling of boundaries, it does not commute and leads to the so called commutation error. In sce-narios with periodic boundary conditions and nearly homogeneous turbulence, this error isnegligible. The numerical analysis of the commutation error in LES is given in [BIL06].A rich variety of SGS models has been developed. The first proposed eddy-viscosity modelis the Smagorinsky model [Joh04]. Via dimensional analysis, it follows for the dissipationrate that (cid:15) ∼ U I L I . The same relation also holds for the filter width ∆ δ and the corresponding characteristicvelocity of the unresolved scales U δ (cid:15) ∼ U δ ∆ δ . It follows that ν SGS ∼ U δ ∆ δ ∼ U I L − / I ∆ / δ . The mixing length assumption is then U I ∼ L I S, where S is the magnitude of the filtered strain rate tensor.One gets by replacing the mixing length assumption into the viscosity relation ν SGS = CL / I ∆ / δ S. The integral length scale is approximated by L I ∼ ∆ δ .Then, the artificial viscosity is defined as ν SGS = ( C S ∆ δ ) S, (4.3.14)where C S is the Smagorinsky coefficient. The filter width is thereby computed via thegeometrical mean ∆ δ = (cid:112) ∆ x ∆ x ∆ x , (4.3.15)where ∆ xi , i = 1 , , C S → C S = C S ( y + ) = C S (1 − e − y + / ) . (4.3.16)It improves the performance of the model in situations with simple geometries, where moreor less the boundary layer theory holds.Throughout all LES computations performed in this thesis, the Smagorinsky model withVan Driest damping was used. 38 Simulation principles and modelling
Variational Multiscale (VMS) approach is a comparatively new method to simulate incom-pressible turbulent flow. The fundamental idea is based on scale separation similar to LES,but referring to the variational framework of the underlying equations. VMS concepts forLES were primarily introduced by Hughes [HMJ00]. Instead of using the filtered governingequations, the weak form of the equations and variational projection into subspaces is thebasic concept behind VMS. Due to the variational formulation, the finite element methodframework may be preferable, although it is also suitable for other discretization techniques.Nowadays, many different classes of VMS methods exist, in this thesis the focus will be onthe projection-based VMS method by [JKM05].The first step is the variational formulation of the Navier-Stokes equations of form ofEquation (2.4.1) and (2.4.2). For the case, we only consider homogeneous Dirichlet andhomogeneous Neumann boundary conditions in the following. For this let ∂ Ω = ∂ Ω D ∪ ∂ Ω N with | ∂ Ω N | >
0. As we have a system of partial differential equations for the velocity andpressure, we define two spaces for our solutions and test functions V = H ,∂ Ω D (Ω , R ) = { v ∈ H (Ω , R ) : v = 0 on ∂ Ω D } . (4.4.1) Q = L (Ω , R ) . (4.4.2)We multiply Equation (2.4.2) with test function v ∈ V and Equation (2.4.1) with q ∈ Q respectively, integrate over the whole domain Ω and integrate by parts, to get (cid:90) Ω ∂u∂t · v dx + (cid:90) Ω νS ( u ) : S ( v ) dx − (cid:90) Ω ( ∇ · v ) p dx + (cid:90) Ω ( u · ∇ ) u · v dx = (cid:90) Ω f · v dx ∀ v ∈ V, − (cid:90) Ω ( ∇ · u ) q dx ∀ q ∈ Q,S ( u ) is the strain rate tensor S ( u ) = 12 ( ∇ u + ∇ u T ) . We consequently define the bilinear forms a ( · , · ) : V × V → R , b ( · , · ) : V × Q → R , thetrilinear form c ( · , · , · ) : V × V × V → R and linear form f ( · ) : V → R : a ( u, v ) = (cid:90) Ω νS ( u ) : S ( v ) dx, (4.4.3) b ( u, q ) = − (cid:90) Ω ( ∇ · u ) q dx, (4.4.4) c ( u, u, v ) = (cid:90) Ω ( u · ∇ ) u · v dx, (4.4.5) f ( v ) = (cid:90) Ω f · v dx. (4.4.6)Generally, the L inner product is defined as( φ, ψ ) = (cid:90) Ω φ · ψ dx ∀ φ, ψ ∈ L . (4.4.7)39 Simulation principles and modelling
The variational problem of the Navier-Stokes equations then reads:
Problem 8 ( Weak formulation of the Navier-Stokes equations ) . Find u ∈ V and p ∈ Q satisfying ( ∂u∂t , v ) + a ( u, v ) + b ( v, p ) + c ( u, u, v ) = f ( v ) ∀ v ∈ V, (4.4.8) b ( u, q ) = 0 ∀ q ∈ Q. (4.4.9)The equations from Problem 8 may also be written in short form as, A (( u, p ) , ( v, q ))= ( ∂u∂t , v ) + a ( u, v ) + b ( v, p ) + c ( u, u, v ) + b ( u, q ) ∀ ( v, q ) ∈ V × Q. (4.4.10)Systems of this form are called saddle-point problem and its analysis is seen in the nextchapter.For VMS methods, the corresponding trial and test spaces are decomposed into three parts,large scales, small scales and unresolved scales V = V L ⊕ V S ⊕ V U , Q = Q L ⊕ Q S ⊕ Q U . For the solution and test functions, we obtain u = u L + u S + u U , p = p L + p S + p U ,v = v L + v S + v U , q = q L + q S + q U . Inserting the scale seperation into Equation (4.4.10), it may be written as a system of threevariational equations A (( u L , p L ) , ( v L , q L )) + A (( u S , p S ) , ( v L , q L )) + A (( u U , p U ) , ( v L , q L )) = f ( v L ) , (4.4.11) A (( u L , p L ) , ( v S , q S )) + A (( u S , p S ) , ( v S , q S )) + A (( u U , p U ) , ( v S , q S )) = f ( v S ) , (4.4.12) A (( u L , p L ) , ( v U , q U )) + A (( u S , p S ) , ( v U , q U )) + A (( u U , p U ) , ( v U , q U )) = f ( v U ) . (4.4.13)As it is not intended to explicitly solve the unresolved scales, Equation (4.4.13) is neglected.Another assumption is that the unresolved scales do not influence the large scales directly,therefore A (( u U , p U ) , ( v L , q L )) = 0 . The influence of the unresolved scales onto the resolved has to be modeled. A widely usedway is to use the eddy-viscosity model as previously mentioned in Equation (4.3.14) A (( u U , p U ) , ( v S , q S )) ≈ (2 ν U S ( u S ) , S ( v S )) . As a result of neglecting the unresolved scales, we get the new system of equations.40
Simulation principles and modelling
Problem 9 ( Three-scale VMS formulation of the Navier-Stokes equations ) . Find u L + u S ∈ V L ⊕ V S and p L + p S ∈ Q L ⊕ Q S such that A (( u L , p L ) , ( v L , q L )) + A (( u S , p S ) , ( v L , q L )) = f ( v L ) , (4.4.14) A (( u L , p L ) , ( v S , q S )) + A (( u S , p S ) , ( v S , q S )) + (2 ν U S ( u S ) , S ( v S )) = f ( v S ) , (4.4.15) for all v L + v S ∈ V L ⊕ V S and q L + q S ∈ Q L ⊕ Q S . A crucial point of VMS methods is the definition of the appropriate spaces for the largescales and small scales . The strategy used in this thesis is a coarse space projection-basedmethod.Consider a standard pair of conforming finite element spaces V h × Q h ⊂ V × Q for allscales of velocity and pressure which fulfills several conditions for saddle-point problems(discussed in detail in the next chapter). In addition, let L be a finite element space ofsymmetric d × d tensor-valued functions L h ⊂ L = { l ∈ L (Ω , R d × d ) : l = l T } . (4.4.16)Let V hL ∈ H (Ω , R ) be the discrete space for the large scales such that the condition L h = { S ( v hL ) : v hL ∈ V hL } ⊆ { S ( v h ) : v h ∈ V h } holds. There are two possibilities of choosingthe coarse finite element space V hL . On the one hand, if for V h a higher order finite elementspace is chosen, one may take a lower order finite element space for V hL on the same grid.This approach is the one-level projection-based VMS method. On the other hand, in thecase of the same order for the resolved and large scales spaces, V hL may be defined ona coarser grid, which is called the two-level projection-based method. Due to reasons ofsimplicity, we stick to the first mentioned approach in this thesis.However, the discrete space V hL does not incoporate boundary conditions, thus generallyspeaking V hL is no subset of V h .Define the projection operator Π V : V h → V hL such that, (cid:0) S ( u h − Π V u h ) , S ( v hL ) (cid:1) = 0 ∀ v hL ∈ V hL . (4.4.17)In addition, let Π L : L → L h be the L -projection from L to L h respectively, (cid:0) S ( u h ) − Π L S ( u h ) , l h (cid:1) = 0 ∀ l h ∈ L h . (4.4.18)One important aspect is that the strain rate tensor of the large scales defined in Equa-tion (4.4.17) equals the large scales of the strain rate tensor defined in Equation (4.4.18).It follows that the definition by projection of the large scales and differentiation commutes.As mentioned previously, this does not hold in classic LES generally. Lemma 1.
Let v h ∈ V h and L h = { S ( v hL ) : v hL ∈ V hL } then it holds Π L S ( v h ) = S (Π V v h ) . (4.4.19)A simple proof is given in [JKM05].Taking the system of equations of Problem 9, reunite the decomposition V h = V hL ⊕ V hS Simulation principles and modelling and Q h = Q hL ⊕ Q hS with the small scale part defined by the projection V hS = ( I − Π V ) V h ,we obtain A (( u h , p h ) , ( v h , q h )) + (2 ν U S ( u hS ) , S ( v hS )) = f ( v h ) ∀ ( v, q ) ∈ V h × Q h . The modeled term may be rewritten as (cid:0) ν U S ( u hS ) , S ( v hS ) (cid:1) = (cid:0) ν U S (cid:0) ( I − Π V ) u h (cid:1) , S (cid:0) ( I − Π V ) v h (cid:1)(cid:1) = (cid:0) ν U ( I − Π L ) S ( u h ) , ( I − Π L ) S ( v h ) (cid:1) = (cid:0) ν U S ( u h ) , S ( v h ) (cid:1) − (cid:0) ν U Π L S ( u h ) , S ( v h ) (cid:1) . Let g h ∈ L h such that g h = Π L S ( u h ), we finally obtain the projection-based VMS method. Problem 10 ( Projection-based VMS formulation of the Navier-Stokes equa-tions ) . Let ( u h , p h , g h ) ∈ V h × Q h × L h , such that ( ∂u h ∂t , v h ) + (cid:0) ν + ν U ) S ( u h ) , S ( v h ) (cid:1) (4.4.20)+ c ( u h , u h , v h ) + b ( v h , p h ) − (cid:0) ν U g h , S ( v h ) (cid:1) = f ( v h ) ∀ v h ∈ V h , (4.4.21) b ( u h , q h ) = 0 ∀ q ∈ Q h , (4.4.22) (cid:0) g h − S ( u h ) , l h (cid:1) ∀ l h ∈ L h . (4.4.23)Refering to the numerical analysis of Problem 10 in [ARJR15]. The principal way ofperforming the analysis is the same as for the Galerkin discretization of the Navier-Stokesequations.In order to obtain an efficient implementation for solving Problem 10, the space L h hasto be a discontinuous finite element space with a L -orthogonal basis. This ensures thatthe mass matrix is a diagonal matrix and its inverse may be computed easily. Using adiscontinuous space for L h makes also sense from the point of view that the functions of L h are L -projections of strain rate tensor of finite element functions, which are usuallydiscontinuous functions too.All in all, the combination of choosing the space L h and the eddy-viscosity model resultsin the projection-based VMS method. Nevertheless, this VMS method is less sensitive tochoice of the eddy-viscosity model than in traditional LES. This is due to the fact that themodeled scales influence much less scales directly in VMS than in LES. On the one hand, if L h = { } the eddy-viscosity model influences all scales, recovering the classic LES model.On the other hand, if L h = { S ( v h ) : v h ∈ V h } the modeled viscosity is switched off andthe Navier-Stokes equations are reobtained. Therefore, a low-order space of L h means thatthe turbulence model has a larger influence and for a higher order space of L h the modelhas less influence.The local turbulence intensity may be estimated with the size of the local small resolvedscales, η T = (cid:107) g h − S ( u h ) (cid:107) L ( T ) , (4.4.24)where T is a triangulation of a domain Ω and T ∈ T . If the size of the small scales is large,many unresolved scales can be expected and vice versa. A method choosing the adaptive42 Simulation principles and modelling projection space may be found in [JK10]. In this thesis, we stick to the conventionalnon-adaptive version of the three-scale projection-based VMS method.43
Continuous Galerkin method for the
Navier-Stokes equations
In this chapter, we present a continuous Galerkin discretization technique for the unsteadyNavier-Stokes equations. Firstly, we consider the discretization of the steady Stokes prob-lem and discuss the characteristic conditions to obtain a stable method. Further on, a timediscretization of the governing equations of motion is introduced.For ease of simplicity and practical reasons, inhomogeneous Dirichlet and homogeneousNeumann boundary conditions are assumed in all derivations.
At first, we consider the weak formulation of Stokes problem and the unsteady Navier-Stokes problem (as already derived in the previous chapter in Problem 8). The appropriatespaces for the velocity and pressure are of the form of Equation (4.4.1) and (4.4.2). Ad-ditionally, the Dirichlet boundary condition is incoporated in the solution space V D asdefined V D = { v ∈ H (Ω , R ) : v = v D on ∂ Ω D } . We obtain the weak formulation of the steady Stokes equations by neglecting the timederivative and convective term of the Navier-Stokes equations.
Problem 11 ( Weak formulation of the steady Stokes equations ) . Find u ∈ V D and p ∈ Q satisfying a ( u, v ) + b ( v, p ) = f ( v ) ∀ v ∈ V, (5.1.1) b ( u, q ) = 0 ∀ q ∈ Q. (5.1.2)Systems of the form seen in Problem 11 are called saddle-point problems, since the solution u, p ∈ V × Q is also a minimizer of J ( u ) = a ( u, u ) + f ( u ) → min, subject to the constraint b ( u, q ) = 0 ∀ q ∈ Q. As for optimization problems with constraints, we apply a Lagrangian L ( u, q ) = J ( u ) + b ( u, q ) . Then we have that each solution u, p ∈ V × Q is a saddle-point of the Lagrangian and itholds L ( u, q ) ≤ L ( u, p ) ≤ L ( v, p ) ∀ ( v, q ) ∈ V × Q. Continuous Galerkin method for the Navier-Stokes equations
The pressure p ∈ Q may be interpret as Lagrangian multiplier associated with the incom-pressibility constraint ∇ · u = 0.In order to prove the existence and uniqueness for elliptic partial differential equations, wecould show the continuity and coercivity of the bilinear form and the lemma of Lax-Milgramwill guarantee the unique solvability. Unfortunately, this holds not for saddle-point prob-lems and an additional condition has to be fullfilled. Brezzi’s theorem for mixed methodsis used therefore.A general mixed variational form involves the two bilinear forms a ( · , · ) : V × V → R and b ( · , · ) : V × Q → R and two linear forms f ( · ) : V → R and g ( · ) : Q → R . Thus, with theseforms we define the following mixed problem a ( u, v ) + b ( v, p ) = f ( v ) ∀ v ∈ V, (5.1.3) b ( u, q ) = g ( q ) ∀ q ∈ Q. (5.1.4)We define the space V of the kernel of the bilinear form b ( · , · ) V = { v ∈ V : b ( v, q ) = 0 ∀ q ∈ Q } . (5.1.5) Theorem 1 ( Brezzi’s Theorem ) . Let a ( · , · ) : V × V → R and b ( · , · ) : V × Q → R betwo given bilinear forms, that fulfill the conditions:(i) The bilinear forms are continuous, a ( u, v ) ≤ α (cid:107) u (cid:107) V (cid:107) v (cid:107) V ∀ u, v ∈ V, (5.1.6) b ( u, q ) ≤ β (cid:107) u (cid:107) V (cid:107) q (cid:107) Q ∀ u ∈ V, q ∈ Q. (5.1.7) (ii) a ( · , · ) is coercive on the kernel, i.e. there exists an α > so that, a ( u, u ) ≥ α (cid:107) u (cid:107) V ∀ u ∈ V . (5.1.8) (iii) The Ladyshenskaja-Babuˇska-Brezzi (LBB) condition of the constraint b ( · , · ) is full-filled, i.e. there exists an β > such that, sup v ∈ V b ( v, q ) (cid:107) v (cid:107) V ≥ β (cid:107) q (cid:107) Q ∀ q ∈ Q. (5.1.9) Then the mixed method from Equation (5.1.3) and (5.1.4) is uniquely solvable and thesolution fullfills the stability estimate, (cid:107) v (cid:107) V + (cid:107) p (cid:107) Q ≤ C ( (cid:107) f (cid:107) V ∗ + (cid:107) g (cid:107) Q ∗ ) . (5.1.10)Here denotes V ∗ and Q ∗ the corresponding dual space to V and Q . For more details werefer to [BBF13]. The proof of Theorem 1 is given in [Sch18] and the analysis of the Stokesproblem are shown in [Led16].The weak formulation of the Navier-Stokes equations is already given in Problem 8.45 Continuous Galerkin method for the Navier-Stokes equations
In the following, we apply the continuous Galerkin discretization for the mixed problem ofform of Equation (5.1.3) and (5.1.4). We define the finite-dimensional subspaces V h ⊂ V , V hD = { v ∈ V h : v = v D on ∂ Ω D } and Q h ⊂ Q . The h refers to a discrete space. Thediscrete form of the variational problem is: Problem 12 ( Discrete formulation of the steady Stokes equations ) . Find u h ∈ V hD and p h ∈ Q h satisfying a ( u h , v h ) + b ( v h , p h ) = f ( v h ) ∀ v h ∈ V h , (5.2.1) b ( u h , q h ) = g ( q h ) ∀ q h ∈ Q h . (5.2.2)Discrete stability of Problem 12 is not inherited from the continuous problem. The conti-nuity of the bilinear forms follows from the infinite dimensional case since V h ⊂ V . Thediscrete kernel ellipticity is a ( u h , u h ) ≥ α h (cid:107) u h (cid:107) V ∀ u h ∈ V h , (5.2.3)where V h is defined as V h = { v h ∈ V h : b ( v h , q h ) = 0 ∀ q h ∈ Q h } . (5.2.4)The discrete LBB-condition readssup v h ∈ V h b ( v h , q h ) (cid:107) v h (cid:107) V ≥ β h (cid:107) q h (cid:107) Q ∀ q h ∈ Q h . (5.2.5)It follows for the Stokes problem the conditionsup v h ∈ V h (cid:82) Ω ( ∇ · v h ) q h dx (cid:107) v h (cid:107) V ≥ β h (cid:107) q h (cid:107) Q ∀ q h ∈ Q h , (5.2.6)which is the constraint that arises for the definitions of the finite spaces. Since the LBB-condition for Stokes equations holds in the continuous level, for a fixed pressure space thevelocity space can be enlarged to get discrete LBB-condition. The enlargement can be doneby increasing the polynomial order or refining the mesh. Therefore, the same polynomialdegree of order for the velocity and pressure space may lead to an unstable discretization.The discrete formulation of the unsteady Navier-Stokes equations is given below. Problem 13 ( Discrete formulation of the Navier-Stokes equations ) . Find u h ∈ V hD and p h ∈ Q h such that ( ∂u h ∂t , v h ) + a ( u h , v h ) + b ( v h , p h ) + c ( u h , u h , v h ) = f ( v h ) , (5.2.7) b ( u h , q h ) = 0 , (5.2.8) is satisfied for all v h ∈ V h and q h ∈ Q h . Continuous Galerkin method for the Navier-Stokes equations
The discretized weak formulation of Problem 12 can only be solved if the velocity u h fulfillsthe incompressibility constraint, (cid:90) Ω ( ∇ · u h ) q h dx = 0 ∀ q h ∈ Q h , which is called a discrete divergence-free property. Nevertheless, this does not generallyhold in a strong formalism.Assume we have the property that the space of the divergence of the test functions of thevelocity space is a subspace of the pressure space, thus ∇ · V h ⊂ Q h . (5.2.9)Then a discrete divergence-free velocity is also exactly divergence-free, namely (cid:90) Ω ( ∇ · u h ) q h dx = 0 ∀ q h ∈ Q h ⇒ ∇ · u h = 0 . This property has many advantages, especially u h leads to a better approximation of thevelocity and proper physical behavior.From [Leh10] and [Led16], the impact of exactly divergence-free velocity approximationmay be seen in the kinetic energy loss. Consider the momentum equation of the Navier-Stokes equations from Equation (2.4.2) with constant density ρ = 1 and no volume force f . Due to the friction of the fluid, the kinetic energy should decrease over time. The rateof change of energy with respect to time is then dEdt = (cid:90) Ω ∂ ( u · u ) ∂t dx = (cid:90) Ω u · ∂u∂t dx = (cid:90) Ω (cid:0) − ν ∇ u : ∇ u − u · (( u · ∇ ) u ) − ( ∇ · u ) p (cid:1) dx. The convective term can be rewritten as (cid:90) Ω u · (( u · ∇ ) u ) dx = 12 (cid:90) Ω ( u · ∇ )( u · u ) dx = − (cid:90) Ω ( u · u )( ∇ · u ) dx. Using the incompressibility constraint, it follows dEdt = (cid:90) Ω − ν ∇ u : ∇ u dx ≤ . This does not automatically hold in the approximate sense dE h dt = (cid:90) Ω − ν ∇ u h : ∇ u h + 12 ( u h · u h )( ∇ · u h ) dx. since ( u h · u h ) / ∈ Q h and u h is only discrete divergence-free.If Equation (5.2.9) is fullfilled, then we obtain the physically correct behavior dE h dt = (cid:90) Ω − ν ∇ u h : ∇ u h dx ≤ . Especially in turbulent flows, where conditions with relatively small viscosity (low moleculardiffusion) and very rapid mixing of the transport quantities (high turbulent diffusion), theproperty of exact divergence-free velocity is essential.47
Continuous Galerkin method for the Navier-Stokes equations
We have seen in the previous section that the used finite spaces have to fulfill the discreteLBB-condition in order to achieve an unique and stable solution.We assume a triangulation T which is regular and consists of elements T ∈ T with thecorresponding set of vertices V and edges/faces E / F .We introduce the most common finite element pairing for the Stokes problem, the Taylor-Hood element [HT73]. This discretization consists of standard H -conforming elements forthe velocity and pressure space, while the polynomial order of the pressure space is onedegree lower (at least linear polynomials for pressure space). The notation P p - P p − is usedfor the finite element pairing, in this case for the Taylor-Hood element.Let P p ( T ) be the space of all element-wise polynomials on T up to degree p . The finiteelement spaces are then chosen as V h = [ P p ( T )] ∩ [ C (Ω)] , (5.3.1) Q h = P p − ( T ) ∩ C (Ω) . (5.3.2)As usual, the test functions for the velocity live in the space V h = { v h ∈ V h : v h = 0 on ∂ Ω D } . (5.3.3)Due to the continuity of the velocity and pressure, the variational formulation of the Navier-Stokes (Problem 8) does not have to be changed for the discretization with Taylor-Hoodelements.The Taylor-Hood element satisfies the discrete LBB-condition for the Stokes Problem 5.2.1and 5.2.2, its proof may be seen in [Che14]. For P - P , it can be shown that these elementshave a quadratic convergence rate if the solution is smooth enough (cid:107) u − u h (cid:107) H (Ω) + (cid:107) p − p h (cid:107) L (Ω) ≤ h | u | H (Ω) + h | p | H (Ω) . (5.3.4)The big drawback of this choice of elements is that it only peserves discrete divergence-freeproperty since ∇ · V h (cid:54)⊂ Q h . After the spatial discretization of the unsteady Navier-Stokes equations (Problem 13), thereare still two aspects to consider to obtain a complete discretization. Firstly, an appropriatediscretization of the time derivative has to be discussed. Secondly, an approach to solvethe nonlinear convective term of the momentum equation. Therefore, an implicit-explicit(IMEX) splitting scheme will be discussed (refer to [ARS97]). For the sake of simplicity,the volume force term f will be neglected here.The approximated velocity u h ( x, t ) and pressure p ( x, t ) is given by u h ( x, t ) = N V (cid:88) i =1 u i ( t ) φ i ( x ) , p h ( x, t ) = N Q (cid:88) i =1 p i ( t ) ψ i ( x ) , Continuous Galerkin method for the Navier-Stokes equations while { φ i ( x ) } N V i =1 and { ψ i ( x ) } N Q i =1 are the basis functions of the finite element spaces V h and Q h . By that we define the matrices M ∈ R N V × N V M ij = (cid:90) Ω φ i ( x ) · φ j ( x ) dx ∀ i, j = 1 ... N V ,A ∈ R N V × N V A ij = (cid:90) Ω ν ∇ φ i ( x ) : ∇ φ j ( x ) dx ∀ i, j = 1 ... N V ,B ∈ R N V × N Q B ij = − (cid:90) Ω (cid:0) ∇ · φ i ( x ) (cid:1) ψ j ( x ) dx ∀ i = 1 ... N V ; j = 1 ... N Q ,C ( u h ) ∈ R N V C i = (cid:90) Ω (cid:0) ( u h · ∇ ) u h (cid:1) · φ i ( x ) dx ∀ i = 1 ... N V . For appropriate initial conditions, the problem is given as:
Problem 14 ( Matrix form of spatial discretized Navier-Stokes equations ) . Find u i ∈ R N V ∀ i = 1 ... N V and p i ∈ R N Q ∀ i = 1 ... N Q satisfying M ∂u i ( t ) ∂t + Au i ( t ) + Bp i ( t ) + C (cid:0) u i ( t ) (cid:1) u i ( t ) = 0 in [0 , T ] ,B T u i ( t ) = 0 in [0 , T ] ,u i ( t = 0) = u i, . For the time discretization, we use the first order IMEX scheme. The main idea is to handlethe convective term explicitly and use it as a kind of force term, while the diffusion andthe incompressibility constraint are processed implicitly.Generally, explicit methods are computationally cheap and can incoporate nonlinearitieswithout solving a nonlinear system. However, they are conditionally stable and may getunstable if the time step is not restricted.Implicit methods have the advantage that they are generally unconditionally stable, butare very expensive since at each time step a system of equations has to be solved.Such decomposition methods are called IMEX splitting schemes. We define the time step∆ t ≥ M u i ( t + ∆ t ) − u i ( t )∆ t + Au i ( t + ∆ t ) + Bp i ( t + ∆ t ) = − C (cid:0) u i ( t ) (cid:1) u i ( t ) ,B T u i ( t + ∆ t ) = 0 , which can be rewritten in matrix form (cid:18) M + ∆ t A ∆ t B ∆ t B T (cid:19) (cid:18) u i ( t + ∆ t ) p i ( t + ∆ t ) (cid:19) = ∆ t (cid:18) − C (cid:0) u i ( t ) (cid:1) u i ( t ) + τ M u i ( t )0 (cid:19) . In order to obtain the residual form, the equations are extended (cid:18) M + ∆ t A ∆ t B ∆ t B T (cid:19) (cid:18) u i ( t + ∆ t ) − u i ( t ) p i ( t + ∆ t ) − p i ( t ) (cid:19) = ∆ t (cid:18) − C (cid:0) u i ( t ) (cid:1) u i ( t ) − Au i ( t ) − Bp i ( t ) − B T u i ( t ) (cid:19) . We define M ∗ = (cid:18) M + ∆ t A ∆ t B ∆ t B T (cid:19) , Continuous Galerkin method for the Navier-Stokes equations and D = (cid:18) − C (cid:0) u i ( t ) (cid:1) − A − B − B T (cid:19) . We obtain the final system of equations (cid:18) u i ( t + ∆ t ) p i ( t + ∆ t ) (cid:19) = ( I + ∆ t M ∗− D ) (cid:18) u i ( t ) p i ( t ) (cid:19) . (5.4.1)After each time step, only the convective part has to be updated with the new velocity, ifthere is no change in the time step size.Unfortunately, for RANS and LES/VMS, the eddy-viscosity gets recalculated in each step,forcing an update of M ∗ , D and therefore M ∗− . The approximation of the RANS/SFNS equations (Problem 2 and 7) is very similar to thediscretization of the Navier-Stokes equations of the previous sections. The only differenceis the diffusion term. In turbulence modelling, the divergence of the Reynolds stress tensor(sub-grid scale tensor in LES/VMS) is approximated as a turbulent diffusion by the Boussi-nesq hypothesis. This new term gets incoporated by the molecular diffusion, resuming ina condensed diffusion term. The total viscosity consists of molecular- and eddy-viscosity ν total = ν + ν T . In the Navier-Stokes equations, ν = const and therefore the diffusion term could be sim-plified to ∇ · (cid:0) ν ( ∇ u + ∇ u T ) (cid:1) = ν ∆ u, with the incompressibility constraint ∇ · u = 0.In the case of the RANS/SFNS equations, this simplification does not hold since ν total isno constant anymore ∇ · (cid:0) ν total ( ∇ u + ∇ u T ) (cid:1) (cid:54) = ν total ∆ u. The bilinear form a ( u h , v h ) then changes to a turb ( u h , v h ) = (cid:90) Ω ν + ν T ) S ( u h ) : S ( v h ) dx, with S ( u h ) = ( ∇ u h + ∇ u Th ). The theory for the incompressible Navier-Stokes equationswith non-constant viscosity is given in [Kai14].The discretization of the two-equation eddy-viscosity models (Problem 3, 4 and 5) using thecontinuous Galerkin method with H -conforming elements is relatively straight forward.For ease of presentation we only consider the approximation of the K − (cid:15) model, thederivation of the other models is analogously.Firstly, as usual we start with the weak formulation. Define the function space for theturbulent kinetic energy and for the dissipation rate R D = H (Ω , R ) = { r ∈ H (Ω , R ) : r = r D on ∂ Ω D } ,R = H (Ω , R ) = { r ∈ H (Ω , R ) : r = 0 on ∂ Ω D } . Continuous Galerkin method for the Navier-Stokes equations
As for the Navier-Stokes equations, we assume Dirichlet and homogeneous Neumann bound-ary conditions ∂ Ω = ∂ Ω D ∪ ∂ Ω N , g K | ∂ Ω N , g (cid:15) | ∂ Ω N = 0 on a bounded domain Ω ⊂ R .We multiply Equation (4.2.5) with test functions k ∈ R and Equation (4.2.6) with testfunctions e ∈ R , integrate over the whole domain and integrate by parts. Then we obtain: Problem 15 ( Weak formulation of the K − (cid:15) equations ) . Find K ∈ R D and (cid:15) ∈ R D such that (cid:90) Ω ∂K∂t k dx + (cid:90) Ω ( (cid:104) u (cid:105) · ∇ K ) k dx = (cid:90) Ω (Π − (cid:15) ) k dx − (cid:90) Ω ( ν + ν T σ K ) ∇ K · ∇ k dx, (cid:90) Ω ∂(cid:15)∂t e dx + (cid:90) Ω ( (cid:104) u (cid:105) · ∇ (cid:15) ) e dx = (cid:90) Ω ( C (cid:15) Π (cid:15)K − C (cid:15) (cid:15) K ) e dx − (cid:90) Ω ( ν + ν T σ (cid:15) ) ∇ (cid:15) · ∇ e dx, is satisfied for all k ∈ R, and e ∈ R . This problem has generally the form of two scalar convection diffusion problems, where thesource and sink terms depend on each other.The weak variational formulation of the steady convection diffusion problem of form ofProblem 6 is seen below.
Problem 16 ( Weak formulation of the scalar convection diffusion problem ) . Find φ ∈ R D such that (cid:90) Ω ( (cid:104) u (cid:105) · ∇ φ ) ψ dx + (cid:90) Ω ( ν + ν T ) ∇ φ · ∇ ψ dx = (cid:90) Ω f ψ dx, (5.5.1) is satisfied for all ψ ∈ R . In order to prove existence and uniqueness of a solution of the weak Problem in Equa-tion (5.5.1), we have to check if the bilinear form is continuous and coercive. Let a CD ( · , · ) : R × R → R be a bilinear form defined as, a CD ( φ, ψ ) = (cid:90) Ω ( (cid:104) u (cid:105) · ∇ φ ) ψ dx + (cid:90) Ω ( ν + ν T ) ∇ φ · ∇ ψ dx. If (cid:104) u (cid:105) = 0 then we obtain a symmetric bilinear form.Let (cid:104) u (cid:105) ∈ L ∞ (Ω) then one gets with the Cauchy–Schwarz inequality, H¨older’s inequality,and the Poincar´e–Friedrichs inequality (refer to [BS02]), a CD ( φ, ψ ) ≤ ( ν + ν T ) (cid:107) ∇ φ (cid:107) L (Ω) (cid:107) ∇ ψ (cid:107) L (Ω) + (cid:107) (cid:104) u (cid:105) (cid:107) L ∞ (Ω) (cid:107) ∇ φ (cid:107) L (Ω) (cid:107) ψ (cid:107) L (Ω) ≤ ( ν + ν T ) (cid:107) ∇ φ (cid:107) L (Ω) (cid:107) ∇ ψ (cid:107) L (Ω) + α P F (cid:107) (cid:104) u (cid:105) (cid:107) L ∞ (Ω) (cid:107) ∇ φ (cid:107) L (Ω) (cid:107) ∇ ψ (cid:107) L (Ω) = α ∗ (cid:107) ∇ φ (cid:107) L (Ω) (cid:107) ∇ ψ (cid:107) L (Ω) = α ∗ | φ | H (Ω) | ψ | H (Ω) , for all φ, ψ ∈ R with α ∗ = ( ν + ν T + α P F (cid:107) (cid:104) u (cid:105) (cid:107) L ∞ (Ω) ). Hence, the bilinear form isbounded.The convective part may be rewritten as, (cid:90) Ω ( (cid:104) u (cid:105) · ∇ ψ ) ψ dx = − (cid:90) Ω ( ∇ · (cid:104) u (cid:105) ) ψ dx. Continuous Galerkin method for the Navier-Stokes equations
Inserting this relation into the bilinear form a CD ( ψ, ψ ) = (cid:90) Ω ( ν + ν T ) ∇ ψ · ∇ ψ −
12 ( ∇ · (cid:104) u (cid:105) ) ψ dx, then for all ψ ∈ R it follows if ( ∇ · (cid:104) u (cid:105) ) ≥ a CD ( ψ, ψ ) ≥ (cid:90) Ω (cid:0) ( ν + ν T ) ∇ ψ · ∇ ψ dx = α (cid:107) ψ (cid:107) R . Thus, a CD ( ψ, ψ ) is coercive. It must be mentioned that for convection-dominated flux, asit is usually the case for turbulent flows, coercivity loss may occur. Although the exactproblem is well-posed, instabilities are possible.The lemma of Lax-Milgram then states that for each bounded functional f ∈ R ∗ thereexists an unique solution φ ∈ R .For the discretization, the standard continuous Galerkin method is used, which just replacesthe space R by R h ⊂ R in the variational formulation. The space R hD incoporates theinhomogeneous Dirichlet boundary condition. Problem 17 ( Discrete formulation of the scalar convection diffusion problem ) . Find φ h ∈ R hD , such that a CD ( φ h , ψ h ) = (cid:90) Ω f ψ h dx, is satisfied for all ψ h ∈ R h . The existence of an unique solution of the discrete problem follows directly from the theoremof Lax–Milgram, since R h is a closed subspace of the Hilbert space R and the propertiesof the bilinear form carry over from R to R h . The use of standard H -conforming finiteelements is sufficient.As for the Navier-Stokes equations, a first order IMEX splitting scheme is used for the two-equation eddy-viscosity models. The diffusive term is handled implicitly and the convectiveterm as well as the source and sink terms are handled explicitly.For each time step, the turbulent quantities of the eddy-viscosity models are updated firstlywith respect to the values of the previous time step. After that, the RANS/SFNS equationsare calculated with the new updated eddy-viscosity.52 Hybrid discontinuous Galerkin method for the Navier-Stokes equations
In the previous chapter, we discussed the concept of the standard continuous Galerkin finiteelement method. It uses an approximation of the weak formulation of the PDE, which isachieved by replacing the infinite dimensional space in which the variational formulation isposed by a finite dimensional subspace. This finite element space normally uses piecewisepolynomials, which are continuous across element interfaces. As we already have observedin the previous chapter, its disadvantage is the conservation property and that no stabi-lization can be used for convection dominated flows.Discontinuous Galerkin method overcomes this problem by using a discretization, whichis continuous on each element but discontinuous across elements. While more degrees offreedom are required, it offers generally more flexibility. Its drawback is the high compu-tational effort, because of larger system of equations with less sparsity due to a lot morecouplings of unknowns.In this chapter we introduce a relatively new discretization method that was introduced byJoachim Sch¨oberl and Christoph Lehrenfeld in [Leh10] and [LS16], which manages to solvethose drawbacks. It uses a hybridized version of the divergence-conforming DG methodfrom [CKS05], called an H (div)-conforming hybrid discontinuous Galerkin (HDG) finiteelement method. Firstly, we want to use a discontinuous finite element approximation of the weak formulationof the Navier-Stokes equations. Therefore, we have to look at the space of element-piecewise H ( T ) functions, which form the broken Sobolev space H ( T ) = { v ∈ L (Ω) , v ∈ H ( T ) ∀ T ∈ T } . Due to the reason that functions are no longer continuous over Ω, applying integration byparts is no longer valid. Therefore, we are allowed to integrate by parts on each element T ∈ T . As the functions do not belong to H (Ω) anymore, we use the interior penaltymethod introduced in [Arn82] to weakly enforce continuity, which would lead to a DGformulation.However, in this method we use a semi-discontinuous approach called H (div)-conformingdiscretization. The basic idea is to decompose the velocity into continuous normal velocitycomponents and discontinuous tangential velocity components over edges of finite elements H (Ω , R ) = { u ∈ H ( T, R ) : (cid:74) u · n (cid:75) = 0 on E ∈ F }∩{ u ∈ H ( T, R ) : (cid:74) u × n (cid:75) = 0 on E ∈ F } , where E is an element of the triangulation skeleton F and (cid:74) · (cid:75) is the jump on a commonedge E of two neighbouring elements T and T (cid:74) u · n (cid:75) E = ( u | T − u | T ) · n E , Hybrid discontinuous Galerkin method for the Navier-Stokes equations and (cid:74) u × n (cid:75) E = ( u | T − u | T ) × n E . The sobolev space for H (div)-conforming functions is defined as H (div , Ω) = { v ∈ L (Ω , R ) : ∇ · v ∈ L (Ω , R ) } . (6.1.1)The compatibility condition for H (div)-conformity is the continuity of the normal compo-nent over element edges. We introduce the finite dimensional space W h = { v Th ∈ [ P p ( T )] : (cid:74) v Th · n (cid:75) = 0 on E ∈ F } , (6.1.2)for that holds W h ⊂ H (div , Ω). Note that here the superscript T should indicate theelement T ∈ T and should not be confused with transposing.For the pressure space, a discontinuous finite element space is used Q h = { q h ∈ P p − ( T ) } ⊂ L (Ω) . (6.1.3)It directly follows that ∇ · W h = Q h and thus the exact divergence-free property is fulfilled.Tangential continuity is not included in the velocity space so we have to account for itin another way. This is done in a hybrid DG formulation. Of course, a standard DGformulation for weakly enforcing continuity of the tangential component of the velocityis possible. Nevertheless, all degrees of freedom of neighbouring elements would coupledirectly. In order to reduce the coupling of the system matrix, we introduce an additionalfinite tangential facet space on the skeleton F F h = { v Fh ∈ [ P p ( F )] : v Fh · n = 0 on E ∈ F } . (6.1.4)The HDG scheme has a computational advantage. Although this comes with additionaldegrees of freedoms, the coupling is significantly reduced since element DOF do not couplewith each other and the same for facet DOF. By static condensation, the resulting linearsystem of equations only accounts the facet DOF. The size of the corresponding schurcomplement is then fairly reduced.Figure 6.1: Normal and tangential continuity of H (div)-conforming hybrid DG method.54 Hybrid discontinuous Galerkin method for the Navier-Stokes equations H (div) -conforming HDG method for theNavier-Stokes problem As for the previous derivations, the boundary of the domain Ω is divided into Dirichletand homogeneous Neumann boundaries. The new finite compound space for the velocityis V h = W h × F h . The following notation u h = ( u Th , u Fh ) ∈ V h for the solution and v h =( v Th , v Fh ) ∈ V h for the test functions is used. Dirichlet boundary conditions are posed onthe facet functions only. Thus, the discrete space is given as V h = { u h ∈ V h : u T,nh = 0 , u Fh = 0 on ∂ Ω D } , (6.2.1) V hD = { u h ∈ V h : u T,nh = u T,nh,D , u Fh = ( u Fh,D ) t on ∂ Ω D } . (6.2.2)The jump of the tangential component on the element to the facet is defined as (cid:74) u th (cid:75) = ( u T,th − u Fh ) , and (cid:74) v th (cid:75) = ( v T,th − v Fh ) . The superscript t and n indicates the tangential and normal velocity component respectively u nh = ( u h · n ) n, and u th = u h − u nh . The continuity of the normal component is automatically fulfilled by the definition of W h .For the viscous part, we define the bilinear form a HDG ( · , · ) : V h × V h → R as in [LS16].First we integrate by parts on each element − (cid:88) T ∈T (cid:90) T ν ∆ u Th · v Th dx = (cid:88) T ∈T (cid:18) (cid:90) T ν ∇ u Th : ∇ v Th dx − (cid:90) ∂T ν ( ∇ u Th n ) · v Th ds (cid:19) , while ( ∇ u Th n ) denotes the matrix vector product of the velocity vector gradient and elementboundary normal vector.We add a consistency term with the facet variables ˆ v h = v T,nh + v Fh (cid:88) T ∈T (cid:90) ∂T ν ( ∇ u Th n ) · ˆ v h ds = (cid:90) ∂ Ω N ν ( ∇ u Th n ) · ˆ v h ds, which is in our case zero (homogeneous Neumann boundary conditions) to the previousequation. Then we obtain (cid:88) T ∈T (cid:18) (cid:90) T ν ∇ u Th : ∇ v Th dx − (cid:90) ∂T ν ( ∇ u Th n ) · (cid:74) v th (cid:75) ds (cid:19) , since v Th − ˆ v h = v T,th − v Fh .We add two additional terms for symmetry and stability, which are all zero for the exact55 Hybrid discontinuous Galerkin method for the Navier-Stokes equations solution. The final bilinear form then reads a HDG ( u h , v h ) = (cid:88) T ∈T (cid:18) (cid:90) T ν ∇ u Th : ∇ v Th dx − (cid:90) ∂T ν ( ∇ u Th n ) · (cid:74) v th (cid:75) ds − (cid:90) ∂T ν ( ∇ v Th n ) · (cid:74) u th (cid:75) ds + (cid:90) ∂T ν αp h (cid:74) u th (cid:75) · (cid:74) v th (cid:75) ds (cid:19) , where α is the stabilization parameter.For the pressure we again integrate by parts on each element and define the bilinear form b HDG ( · , · ) : V h × Q h → R such as b HDG ( u h , q h ) = − (cid:88) T ∈T (cid:90) T ( ∇ · u Th ) q h dx. The final discrete H (div)-conforming HDG formulation of the Stokes problem is: Problem 18 ( Discrete HDG formulation of the steady Stokes equations ) . Find u h ∈ V hD and p h ∈ Q h such that a HDG ( u h , v h ) + b HDG ( v h , p h ) = f ( v h ) ∀ v h ∈ V h .b HDG ( u h , q h ) = 0 ∀ q h ∈ Q h . To show well-posedness of the method, we will make use of Brezzi’s theorem for saddle-point problems again. By using this theorem, the necessary conditions are the coercivityof the bilinear form a HDG and the LBB-condition of b HDG . For the continuous problemwe already showed that these conditions hold.In [Leh10] it is shown that for a sufficiently large stabilization parameter the coercivitycondition holds independently of h . Using the norm (cid:107) v h (cid:107) HDG = (cid:88) T ∈T (cid:0) (cid:107) ∇ v h (cid:107) L ( T ) + p h (cid:107) (cid:74) v th (cid:75) (cid:107) L ( ∂T ) (cid:1) , the discrete LBB-condition issup v h ∈ V h b HDG ( v h , q h ) (cid:107) v h (cid:107) HDG ≥ β (cid:107) q h (cid:107) L ∀ q h ∈ Q h , and the paramater β is independent of h .In [Led16], a p-version of the discrete LBB-condition is given.The discretization of the convective term of the Navier-Stokes equations uses an upwindstabilization technique. Due to the normal continuity of the velocity, it is only required tohave to treat the tangential part in the upwind fashion. The upwind function u uph is definedas u uph = u T,nh + (cid:40) u Fh , w · n < u T,th , w · n ≥ . Hybrid discontinuous Galerkin method for the Navier-Stokes equations (a) (b)
Figure 6.2: The HDG upwind scheme. The curved arrows are the wind and the blue linedtriangle is the reference element. (a): At the inflow edge (red) the facet variableis taken. (b): At the outflow edge (red) the tangential component of the elementvariable is taken.where w · n denotes the normal component of the wind of the convection. This means thatthe value which comes from the direction where the wind originates is chosen. On outflowedges of the element boundary, we choose the tangential element value, but on inflow edgesthe facet value is taken.The DG upwind scheme is derived by partial integration on each element and choosing theupwind value for the element boundary integral. The bilinear (trilinear) form c HDG ( w, · , · ) : V h × V h → R as c HDG ( w, u h , v h ) = (cid:88) T ∈T (cid:18) − (cid:90) T ( w ⊗ u Th ) : ∇ v Th dx + (cid:90) ∂T ( w · n ) u uph · v Th ds (cid:19) . Obviously, the unknowns of different elements do not couple at all, because facet unknownsjust couple with one neighbouring element, the downwind element. At the outflow of theelement boundary, an additional constraint is applied to overcome this problem. Thisconstraint glues the facet values on the trace of the upwind element (in a weak sense). Theadditive constraint reads (cid:88) T ∈T (cid:90) ∂T out ( w · n )( u Fh − u T,th ) · v Fh ds, while ∂T out denotes the outflow edges.Thus, the final HDG bilinear form for the convective part is c HDG ( w, u h , v h ) = (cid:88) T ∈T (cid:18) − (cid:90) T ( w ⊗ u Th ) : ∇ v Th dx + (cid:90) ∂T ( w · n ) u uph · v Th ds + (cid:90) ∂T out ( w · n )( u Fh − u T,th ) · v Fh ds (cid:19) . Finally we obtain the spatial discretization of the steady Navier-Stokes equations by settingthe wind as the velocity itself. 57
Hybrid discontinuous Galerkin method for the Navier-Stokes equations
Figure 6.3: The facet is glued to the outflow boundary.
Problem 19 ( Discrete HDG formulation of the Navier-Stokes equations ) . Find u h ∈ V hD and p h ∈ Q h such that, a HDG ( u h , v h ) + b HDG ( v h , p h ) + c HDG ( u h , u h , v h ) = f ( v h ) ∀ v h ∈ V h ,b HDG ( u h , q h ) = 0 ∀ q h ∈ Q h . The time discretization for the H (div)-conforming HDG method uses the same IMEXsplitting technique as already described in the previous chapter. H (div) -conforming finite elements The most common examples for H (div)-conforming finite elements are the Raviart-Thomas RT elements from [RT06] and the Brezzi-Douglas-Marini BDM elements, see [BDM85]. Aswe have seen before that H (div , Ω) vector fields have a continuous normal component anda discontinuous tangential component. For ease of use and simplicity, we assume Ω ⊂ R in this section. For the construction of higher order H (div)-conforming elements we referto [Zag06].The RT p element of order p ≥ RT p ( T ) = [ P p ( T )] + x ˆ P p ( T ) , with the number of degrees of freedom per element dim (cid:0) RT p ( T ) (cid:1) = ( p + 1)( p + 3) andˆ P p ( T ) denotes homogeneous polynomials.The lowest order RT element is defined by one degree of freedom per edge E i of the element(low order edge-based DOF), such that N E : ψ → (cid:90) E i ψ · n ds. By the definition of the barycentric coordinates, the shape functions associated with theedge is ψ E ( λ , λ ) = λ ∇ × (cid:18) λ λ (cid:19) − λ ∇ × (cid:18) λ λ (cid:19) . The global RT finite element space is then RT ( T ) = span { ψ E : ∀ E ∈ F } ⊂ H (div , Ω) . Hybrid discontinuous Galerkin method for the Navier-Stokes equations
The
BDM p element of order p ≥ BDM p ( T ) = (cid:8) v ∈ [ P p ( T )] : (cid:74) v · n (cid:75) = 0 on E ∈ F } . with the number of degrees of freedom per element dim (cid:0) BDM p ( T ) (cid:1) = ( p + 1)( p + 2).The linear BDM is defined by one additional degree of freedom per edge E i of the element(high order edge-based DOF), namely N E : ψ → (cid:90) E i ψ · nv ds ∀ v ∈ P ( E i ) . The corresponding shape functions are ψ E ( λ , λ ) = ∇ × (cid:18) λ λ λ λ (cid:19) . The global
BDM finite element space is then BDM ( T ) = span { ψ E , ψ E : ∀ E ∈ F } ⊂ H (div , Ω) . Indeed, there holds ∇ · ψ E = 0 , N E ( ψ E ) = 0 ∀ E ∈ ∂T. Figure 6.4: Degrees of freedom for the first order BDM element.Both elements fulfill the exact divergence-free condition ∇ · RT ( T ) = ∇ · BDM ( T ) = P ( T ) . But the approximation is better for
BDM as[ P ( T )] ⊂ RT ( T ) ⊂ BDM ( T ) = [ P ( T )] . However, the RT elements need less degrees of freedom to achieve the divergence-freecondition.For both finite elements, higher order versions exist. The construction of the higher order59 Hybrid discontinuous Galerkin method for the Navier-Stokes equations shape functions mimics the exact sequence property of the spaces H , H (curl), H (div) and L called de Rham Complex (see [Zag06]). As we have seen above, the BDM elementuses low and high order edge-based DOF. For p >
1, high order cell-based divergence-freeand non divergence-free DOF are used.One very interesting remark about the physical interpretation of the space separation wasmentioned in [Leh10]. As already known, in turbulent flow eddies of different sizes occur.In a discretized domain, the largest eddies may cover one or more vertices, while smallerwhirls are located on an edge. The smallest ones lie within one element. To resolve thelargest eddies, the lowest order RT element is sufficient. Higher order elements properlyrepresent the eddies located at one edge or within an element. The same issue as already explained in the Section 5.5 applies as well for the HDG scheme.Instead of using the Laplacian of the velocity ν ∆ u , the strain rate tensor is applied ∇ · (cid:0) ( ν + ν T )( ∇ u + ∇ u T ) (cid:1) . Therefore, only the bilinear form a HDG ( u h , v h ) changes to a HDG,turb ( u h , v h ) = (cid:88) T ∈T (cid:18) (cid:90) T ν + ν T ) S ( u Th ) : S ( v Th ) dx − (cid:90) ∂T ν + ν FT )( S ( u Th ) n ) · (cid:74) v th (cid:75) ds − (cid:90) ∂T ν + ν FT )( S ( v Th ) n ) · (cid:74) u th (cid:75) ds + (cid:90) ∂T ν + ν FT ) αp h (cid:74) u th (cid:75) · (cid:74) v th (cid:75) ds (cid:19) , with S ( u h ) = ( ∇ u h + ∇ u Th ). The variable ν FT corresponds to the eddy-viscosity calculatedby the facet variable of the transported turbulent quantities.The HDG version of the two-equation turbulence models for RANS simulation are veryanalogously to the HDG version of the scalar convection diffusion equation. The derivationof this scheme is very similar to Stokes and Navier-Stokes problem and is fairly detailedexplained in [Leh10].We define the compounded finite space for the transported quantity R h = { ( φ, φ F ) : φ ∈ P p ( T ) , φ F ∈ P p ( F ) } , and the two bilinear forms a CD,HDG : R h × R h → R and c CD,HDG : R h × R h → R a CD,HDG ( φ h , ψ h ) = (cid:88) T ∈T (cid:18) (cid:90) T ( ν + ν T ) ∇ φ h · ∇ ψ h dx − (cid:90) ∂T ( ν + ν FT )( ∇ φ h · n ) (cid:74) ψ h (cid:75) ds − (cid:90) ∂T ( ν + ν FT )( ∇ ψ h · n ) (cid:74) φ h (cid:75) ds + (cid:90) ∂T ( ν + ν FT ) αp h (cid:74) φ h (cid:75)(cid:74) ψ h (cid:75) ds (cid:19) , Hybrid discontinuous Galerkin method for the Navier-Stokes equations with (cid:74) φ h (cid:75) = φ h − φ Fh and (cid:74) ψ h (cid:75) = ψ h − ψ Fh and c CD,HDG ( φ h , ψ h ) = (cid:88) T ∈T (cid:18) − (cid:90) T (cid:104) u (cid:105) φ h ∇ ψ h dx + (cid:90) ∂T ( (cid:104) u (cid:105) · n ) φ uph ψ h ds + (cid:90) ∂T out ( (cid:104) u (cid:105) · n ) (cid:74) φ h (cid:75) ψ Fh ds (cid:19) , as the upwind value φ uph is defined as φ uph = (cid:40) φ Fh , (cid:104) u (cid:105) · n < φ h , (cid:104) u (cid:105) · n ≥ . Adding diffusion and convection together, we obtain the final equation:
Problem 20 ( Discrete HDG formulation of the steady scalar convection dif-fusion equations ) . Find φ h ∈ R hD such that a CD,HDG ( φ h , r h ) + c CD,HDG ( φ h , r h ) = f ( r h ) ∀ r h ∈ R h The two-equation models are adapted from Problem 20 and the time derivative and corre-sponding source and sink terms are appended.61
Numerical test case
This chapter is dedicated to the setup and results of the plane channel flow test case.This case demonstrates the accuracy and validity of the different modelling principles anddiscretization techniques for wall-bounded turbulent flows. Furthermore, its relatively lesscostly computational effort is also beneficial.One main focus of this thesis is on the comparison of the standard discretization methodand the HDG approach previously described in Chapter 5 and 6.
Firstly, the simulation of fully developed turbulent flow in a plane channel at frictionReynolds number Re τ = 395 is considered. The case consists of two infinite parallel platesbordering an equilibrium flow. In order to approximate this configuration, a finite sub-domain is taken and periodic boundaries are applied in streamwise x and spanwise x directions. In normal direction to the walls x , no-slip Dirichlet boundary conditions arechosen for the velocity. The overall dimensions of the computational domain are 4 δ × δ × δ ,allowing sufficiently large structures to envelope. For all computations, a channel half widthof δ = is taken and the domain was discretized in form of uniform hexahedral meshes.All results given in this thesis are also compared to the well-resolved DNS dataset per-formed by Moser [MKM99] with the same friction Reynolds number. The results of thisDNS benchmark may be considered as an exact solution of the Navier-Stokes equations forthe purpose of this thesis.All channel flow simulations have the following properties:Description Quantity Valuefriction Reynolds number Re τ U b Re b ν . × − friction velocity u τ . U b through the channel is adjusted tobe equal to the value in Table 7.1 by varying an imposed streamwise pressure gradient.The time-averaged value of this quantity is equivalent to the mean wall shear stress.As initial conditions for the velocity, a sightly pertubated laminar parabolic velocity profileis taken (except for RANS). Normally, the laminar flow takes many flow-through timesbefore small disturbances produced by numerical errors trigger the transition to the turbu-lent state. Random pertubations seem very ineffective, the pertubated flow field does not62 Numerical test case obey divergence-free velocity and has no structure. To accelerate the laminar-turbulenttransition, a method by [DV06] creates initial wavelike structures in the near-wall region.These structures have the statistical characteristics of the near-wall streaks described inChapter 3.6.2 and their interaction with the superposed laminar channel flow cause linearinstabilities.The transition to turbulence occurs relatively rapidly and is shown for different times andplanes in Figure 7.2. In the first few time units, the streaks form to a very regular pat-tern. After about three flow-through times, the flow is becoming chaotic and the near-wallstructures strongly effect the mean flow. The laminar flow regime completely breaks downand vortical structures and eddies are becoming widespread in the flow domain. At t ≈ y + -positions is shown. The three differentpositions are y + = 5 (viscous sublayer), y + = 40 (buffer layer) and y + = 100 (logarithmiclayer). The energy spectrum is calculated from the discrete Fourier transformation of thevelocity fluctuations in streamwise and spanwise direction.The following sections are divided into the four different simulation principles of turbulentflow. Figure 7.1: Fully turbulent channel flow.63 Numerical test case (a) t = 3(b) t = 4(c) t = 5(d) t = 8 Numerical test case (e) t = 14(f) t = 25(g) t = 50 Figure 7.2: Velocity magnitude at different times and planes. Left figures show x , x -planeat x = 1 and right figures show x , x -plane at x = 0 . Numerical test case
Three different two-equation eddy-viscosity models from Chapter 4 are used for calculatingthe RANS equations. The K − (cid:15) model has been left out here, since we did not obtainvery good results in this case with this model and the superiority of the K − ω models isunsurpassable in wall-bounded flows. Additionally, a second version of the K − ω model byPeng [PYF15] was introduced. This model incoporates a weakend form of the cross-diffusionterm and different damping functions. The following values for the respective constantsand the corresponding damping functions have been used throughout all computations. • K − ω version from 1998:Description Value/Function C µ . Re ∗ Re ∗ β ∗ . (cid:0) +( Re ∗ ) Re ∗ ) (cid:1) C ω . (cid:0) + Re ∗ . Re ∗ . (cid:1) C ω . Re ∗ Kνω σ K σ ω • K − ω version from Peng [PYF15]:Description Value/Function C µ .
025 + (1 − e − ( Re ∗ ) )(0 .
975 + 0 . e − ( Re ∗ ) ) /Re ∗ β ∗ . (cid:0) − . e − ( Re ∗ ) (cid:1) C ω . (cid:0) . e − ( Re ∗ . ) (cid:1) C ω . Re ∗ Kνω σ K . σ ω . ν T K ∇ K · ∇ ω is added to the right hand-side of the ω equa-tion. • K − ω SST: 66
Numerical test case
Description Value/Function a . β ∗ . min (Π , . Kω ) C ω χ ( , . C ω χ (0 . , . C ω . σ K χ (0 . , σ ω χ (0 . , . K − ω SST is not using any damping-functions and the constants are calculated viathe blending function Equation (4.2.16).The initial field for K and ω is given as K = 1 . T int ) ,ω = √ K l T , where the turbulent intensity is initially set to T int = 0 .
05 and the turbulent length scaleto l T = 0 . δ .The Dirichlet boundary conditions at the solid wall are K D = 0 ,ω D = 10 . The high value of ω at the wall demands fine meshing and small time stepping, especiallyin the first phase of the flow development. Strictly speaking, ω is infinity directly at thewall but this is actually not feasible.Figure 7.3: Mean velocity magnitude field.In the case of the channel flow with streamwise and spanwise periodic boundaries, theReynolds-averaged equations would even allow a one-dimensional consideration of the prob-lem. However, we choose a domain with dimensions 2 δ × δ in x and x direction, sinceit is more represantive for this setup. The use of stretched uniform mesh in x -direction67 Numerical test case allows a very high grid density in the vicinity of the walls with nearest-wall cell spacing of∆ y + = 1.The polynomial order of the Taylor-Hood elements is Q − Q . For the HDG version, BDM elements with tangential finite element space of the same order and the finite spaceof piecewise linear discontinuous functions for the pressure.Once a steady-state is achieved, the corresponding quantities are obtained and the resultsare seen below. In Figure 7.3, the magnitude of the mean velocity field is shown.The maintained mean streamwise velocity profiles compare relatively well with the DNSdata, as shown in Figure 7.4. The agreement in the viscous sublayer is very good. Although,all models significantly suffer in the buffer layer region, since this area is of greatest difficultyto correctly model. As expected, the K − ω SST model gives the best result. Furthermore,in Figure 7.5 the normalized turbulent kinetic energy K and the shear stress componentof the modeled RST are shown. Especially, the Reynolds stress term shows a suprisinglygood match with the DNS data for all model approaches. In the last figure, the viscosityratio is given. It can be seen, that the different curves of the models considerably divergeafter about y + ≈ K − ω SST even higherthan 50.One remarkable notice is that the results obtained from the standard CG discretizationwith Taylor-Hood elements (marked as TH in the figures) and from the HDG approachdo not differ at all. This is mainly due to two reasons. Firstly, the mean velocity field isstill very uniform and any form of vortical motions and eddies are missing. The absenceof such turbulent structures allows the standard CG method to perform quite well sincethe conservation error might be relatively small. Secondly, as the turbulent diffusion iscomparatively much higher than the viscous one, the modeled eddy-viscosity is more thanone order of magnitude higher than the kinematic viscosity as it can be seen in Figure 7.5c.Therefore, the highly increased total viscosity gives (locally) an appearance of a laminarflow.Overall, the outcome of the various RANS models has a reasonably good agreement withthe DNS data and all models work very well for both types of discretization techniques.68
Numerical test case (a)(b)(c)
Figure 7.4: Normalized mean streamwise velocity profiles calculated with (a): K − ω v.1998,(b): K − ω Peng [PYF15], (c): K − ω SST model and compared with DNS datafrom Moser [MKM99]. 69
Numerical test case (a)(b)(c)
Figure 7.5: (a): Turbulent kinetic energy, (b): shear stress component of the modeled RSTand (c): eddy-viscosity, normalized by the squared friction velocity or kinematicviscosity, calculated with RANS turbulence models and compared with DNSdata from Moser [MKM99]. 70
Numerical test case
Within the scope of this thesis, a fully-resolved DNS of the turbulent channel case for Re τ = 395 would have exceeded the computational effort. Nevertheless, a so called quasiDNS or implict LES of sufficiently well resolution gives surprisingly good results and there-fore the outcome is discussed in this work.As initially expected, a DNS calculation with the H -conforming method (without any ad-ditional stabilization technique) early becomes very unstable and obtaining a stable steady-state solution is not possible. The forming of highly three-dimensional vortical structuresleads to a strong mixing of the flow quantities. The property of only discrete divergence-free velocity is not sufficient in such highly turbulent diffusive processes. It is observed,that the divergence of the velocity is far from zero, reaching one order of magnitude atmaxmimum. From a physical point of view, it seems that the conservation of momentumis kind of ”lost” because of the insufficiently satisfied incompressibility constraint.However, the HDG discretization method works very well for DNS of the channel flow andits outcome is shown in the next figures. Again, a stretched uniform hexahedral mesh with BDM elements is used and the dimensions of the grid are 20 × ×
15. The wall-adjacentcell height is ∆ y + = 2.By the time it reaches a statistical steady-state of the fully turbulent flow, time and spatialaveraging is done and the following mean streamwise velocity, as seen in Figure 7.6a, isobtained. The velocity profile agrees very well with the DNS benchmark within all regionsof the turbulent boundary layer.In Figure 7.7, each component of the RST is given, which compares relatively well withthe reference data. The normal Reynolds stresses v (cid:48) v (cid:48) and w (cid:48) w (cid:48) give slightly overpredictedvalues in the whole y + range. This effect is consistent with published results for under-resolved meshes, and is primarily due to excessive resolved scale motion. Especially, theshear stress component shows good conformity. The averaged trace of the RST is seen viathe turbulent kinetic energy in Figure 7.6b.By looking at the energy spectrum of the resolved fluctuations in Figure 7.8, a better ideaof the effect of the resolution of the turbulent scales on the frequency range can be observed.The plots are a one-dimensional spectral representation of the turbulent kinetic energy atdifferent planes in the flow, corresponding to its y + values. The dark dotted line is the k − power curve, which corresponds to the gradient of the inertial range by Kolmogorovhypothesis, while the light dotted line represents the k − curve associated with the inverseenergy cascade.Generally, it is clear that most of the turbulent energy is expressed as lower frequencyeddies. Up to k < , the obtained energy spectrum from the quasi DNS is in good agree-ment with the benchmark spectra for almost all positions. None of the profiles exhibits awell developed inertial range, mostly because the region where isotropic turbulence dom-inates is fairly small in a channel flow. In the streamwise direction, however, all different y + positions tend toward the predicted k − slope. As seen in all plots, after about k > a significant drop in the energy spectrum is observed. The fact that the coarser mesh doesnot reproduce the higher frequency eddies very well is in line with the drop in the corre-sponding curve. In spanwise direction, the drop seems less pronounced as in the streamwisedirection, meaning that the turbulent structures are slightly better resolved there.71 Numerical test case (a)(b)(c)
Figure 7.6: (a): Mean streamwise velocity, (b): turbulent kinetic energy and (c): shearstress component of the RST, normalized by the squared friction velocity, cal-culated with DNS and compared with DNS data from Moser [MKM99].72
Numerical test case (a) (b)(c)
Figure 7.7: The different diagonal components of the RST normalized by the squared fric-tion velocity, calculated with DNS and compared with DNS data from Moser[MKM99]. 73
Numerical test case (a) y + = 5 (b) y + = 40(c) y + = 100(d) y + = 5 (e) y + = 40(f) y + = 100 Figure 7.8: (a), (b), (c): Total energy spectrum in streamwise direction, (d), (e), (f): to-tal energy spectrum in spanwise direction, normalized by the squared frictionvelocity and channel half width, calculated with DNS and compared with DNSdata from Moser [MKM99]. 74
Numerical test case
LES computations are conducted using the Smagorinksy model from Equation (4.3.14)combined with the Van Driest damping function. The rightful choice of the Smagorinksyconstant C S seemed to be not trivial and different values have been suggested in literature.By using the value of C S = 0 .
1, the filtered flow appears to be very diffusive and stronglydamping any fluctuations. Therefore, the parameter is decreased to C S = 0 .
05 and appliedto all LES and VMS simulations.Since we still resolve large eddies and coherent structures in LES, the use of the standardCG discretization method with Taylor-Hood elements again shows very bad results andnumerical instabilities are highly likely to occur. Thus, no steady-state in statistical sensewas achieved and any results are omitted here. The H (div)-conforming HDG method per-formes very well for LES.Compared to the DNS case from the previous section, the same type of grid and finiteelments of same order are used but the mesh resolution is reduced to 10 × ×
10 with anearest-wall cell spacing of ∆ y + = 4. This mesh resolution is fairly coarse even for a LES,but despite that the final results are reasonably fine.The mean velocity profile given in Figure 7.9a deviates slightly in the viscous and loga-rithmic layer region with the DNS benchmark. The reasons are mainly that the excessiveturbulent eddy-viscosity damps the near-wall eddies and the coarse near-wall resolution isincapable of carrying the fine turbulence producing features.The turbulent kinetic energy produced by LES (Figure 7.9b) overestimates in the regionof the viscous and buffer layer and underrates in higher regions. This outcome verifies thesmall deviations in the streamwise mean flow. The total Reynolds shear stress consists ofthe sum of the shear stress component of the resolved RST and averaged modeled SGStensor as given in Figure 7.9c. It can be clearly seen that the resolved part makes up thelargest percentage of the total Reynolds stress. Further on, each diagonal component of theresolved RST is shown in Figure 7.10. The low values of the v (cid:48) v (cid:48) and w (cid:48) w (cid:48) components aremainly due to the excessive damping and shortcomings of assuming isotropy turbulence ofthe SGS turbulence model.As it was expected, the energy spectrum of the coarser LES mesh resolves less vorticalstructures than the DNS from the previous section. The cut-off wave number is set to alower value therefore and the drop in the curve starts at smaller wave numbers. This isconsistent with the obtained spectral representation.75 Numerical test case (a)(b)(c)
Figure 7.9: (a): Mean streamwise velocity, (b): turbulent kinetic energy and (c): resolvedand modeled shear stress components, normalized by the squared friction veloc-ity, calculated with LES and compared with DNS data from Moser [MKM99].76
Numerical test case (a) (b)(c)
Figure 7.10: The different diagonal components of the RST normalized by the squared fric-tion velocity, calculated with LES and compared with DNS data from Moser[MKM99]. 77
Numerical test case (a) y + = 5 (b) y + = 40(c) y + = 100(d) y + = 5 (e) y + = 40(f) y + = 100 Figure 7.11: (a), (b), (c): Total energy spectrum in streamwise direction, (d), (e), (f): totalenergy spectrum in spanwise direction, normalized by the squared frictionvelocity and channel half width, calculated with LES and compared with DNSdata from Moser [MKM99]. 78
Numerical test case
For the VMS method, the same mesh and eddy-viscosity model for the unresolved scales asfor LES is used. The polynomial order of the L h is chosen as p = 1 and therefore assumingless turbulent activities in the whole computational domain. In this case as well, thestandard discretization method with the described finite element pairing from Section 5.3clearly failed to correctly predict a stable solution and the HDG approach is in all aspectsahead.For this choice of the polynomial order for L h , the VMS method performes significantlybetter than the LES with the same setup. The profile of the mean velocity is shown inFigure 7.12a and basically agrees pretty good with the DNS reference in all regions.The second order statistics also show better conformity with the benchmark than with thetraditional LES method. Still, in the vicinity of the wall the u (cid:48) u (cid:48) component of the RST issignificantly overestimated as seen in all turbulence simulation principles before. A betterresolution in the vicinity of the walls would improve this behavior. However, the normalstresses v (cid:48) v (cid:48) and w (cid:48) w (cid:48) shows astonishing good agreement with the DNS reference.As we know, in VMS the large and small scales are resolved and the impact of the modeledunresolved scales only influences the small scales. The choice of the large scale deformationtensor space determines the amount of the resolved small scales among all resolved scalesand therefore restricts the influence of the modeled unresolved structures. For p = 1, theeffect of the modeled scales can be observed in Figure 7.12c for the Reynolds shear stress.There it can be clearly seen, that on average the impact of the unresolved to small scales isdiminishing for such higher order. Although, except to the small peak in the buffer region,the total u (cid:48) v (cid:48) stress coincides with the reference solution.The streamwise and spanwise energy spectrum is given in Figure 7.14. As obvious, theVMS method adequately resolves less scales than the quasi DNS, but shows little bit betterresults than the LES. The reason for that is mainly the less excessive damping of the modeland therefore the better approximation of the statistics.A comparison of the normalized mean streamwise velocity of all simulation principles isgiven in Figure 7.15. 79 Numerical test case (a)(b)(c)
Figure 7.12: (a): Mean streamwise velocity, (b): turbulent kinetic energy and (c): re-solved and modeled shear stress components, normalized by the squared fric-tion velocity, calculated with VMS and compared with DNS data from Moser[MKM99]. 80
Numerical test case (a) (b)(c)
Figure 7.13: The different diagonal components of the RST normalized by the squaredfriction velocity, calculated with VMS and compared with DNS data fromMoser [MKM99]. 81
Numerical test case (a) y + = 5 (b) y + = 40(c) y + = 100(d) y + = 5 (e) y + = 40(f) y + = 100 Figure 7.14: (a), (b), (c): Total energy spectrum in streamwise direction, (d), (e), (f): totalenergy spectrum in spanwise direction, normalized by the squared frictionvelocity and channel half width, calculated with VMS and compared withDNS data from Moser [MKM99].82
Numerical test case
Figure 7.15: Comparison of the normalized mean streamwise velocity of all simulation prin-ciples and the DNS data from Moser [MKM99].83
Conclusion
In this thesis, simulations of the turbulent plane channel flow with the H (div)-conforminghybrid discontinuous Galerkin method have provided insight into the capabilities and ef-ficiency of this relatively new discretization method in order to predict wall-bounded tur-bulent flows. Basically, this was done by attempting to compare its result with a standard H -conforming method with the well-known Taylor-Hood pairing.These discretization techniques have been applied to the four main principles of simulat-ing incompressible turbulent flow. Numerical experiments clearly showed the supremacyof the HDG scheme in resolving turbulent coherent structures and noisy eddies. A goodcoincidence of both types has been observed for the RANS case. Herein, no improve-ment of performance with respect to the new method have been noticed and the standardmethod worked faultless. For LES/VMS and DNS, computations with the conventionalmethod with the Taylor-Hood elements have arised stability issues and the property ofonly discrete divergence-free velocity has been shown to be not sufficient. The quasi DNSsimulation has provided surprisingly good agreement with the reference data, even thoughthe smaller scales have not been adequately resolved. A separation of the turbulent scalesleads to the LES/VMS approach. Both modelling principles have proven their abilities ofproviding good approximations in this numerical test case. Within this comparison, theVMS slightly outperformed the traditional LES method through the different definiton ofthe scale separation.On this basis, we conclude that relatively to its computational effort, the H (div)-conformingHDG method produces qualitatively very good results compared to the given benchmarkcase for all principles. Since VMS methods are quite new in the field of turbulence simulation, future researchcould continue in investigating this method combined with the HDG discretization tech-nique. The used eddy-viscosity models for modelling the unresolved scales allows furtherresearch and modifications. As well as, to give more insight into the correlation betweenthe used model and the space of the strain rate tensor of the large scales.Further on, interesting questions for future research could examine more advanced and im-proved HDG methods for incompressible turbulent flow. In the work of Lederer [LGS18],a new formulation of the Navier-Stokes equations within a HDG scheme was posited. Ex-amination of the application of the new method to turbulence would possibly bring furtherimprovements. 84 ist of Figures k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.2 Schematic representation of the turbulent channel flow. . . . . . . . . . . . 173.3 Normalized viscous shear stress profile and turbulent shear stress profile. . . 183.4 Comparison of DNS data ( Re τ = 395) with the universal velocity profile. . . 193.5 Illustration of laminar-turbulent transition process. . . . . . . . . . . . . . . 203.6 Low-speed streaks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.7 (a): Normalized energy balance Π − (cid:15) . (b): Normalized energy flux Γ . . . . 244.1 Comparison of velocity signal u of different turbulence simulation principlesover time t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.2 Comparison unfiltered and filtered (top-hat, sharp spectral and Gaussianfilter kernel) velocity signal with filter width ∆ δ = . (a): Time domain.(b): Frequency domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.3 Comparison unfiltered and filtered energy spectrum. . . . . . . . . . . . . . 356.1 Normal and tangential continuity of H (div)-conforming hybrid DG method. 546.2 The HDG upwind scheme. The curved arrows are the wind and the bluelined triangle is the reference element. (a): At the inflow edge (red) the facetvariable is taken. (b): At the outflow edge (red) the tangential componentof the element variable is taken. . . . . . . . . . . . . . . . . . . . . . . . . . 576.3 The facet is glued to the outflow boundary. . . . . . . . . . . . . . . . . . . 586.4 Degrees of freedom for the first order BDM element. . . . . . . . . . . . . . 597.1 Fully turbulent channel flow. . . . . . . . . . . . . . . . . . . . . . . . . . . 637.2 Velocity magnitude at different times and planes. Left figures show x , x -plane at x = 1 and right figures show x , x -plane at x = 0 .
5. . . . . . . . 657.3 Mean velocity magnitude field. . . . . . . . . . . . . . . . . . . . . . . . . . 677.4 Normalized mean streamwise velocity profiles calculated with (a): K − ω v.1998, (b): K − ω Peng [PYF15], (c): K − ω SST model and comparedwith DNS data from Moser [MKM99]. . . . . . . . . . . . . . . . . . . . . . 697.5 (a): Turbulent kinetic energy, (b): shear stress component of the modeledRST and (c): eddy-viscosity, normalized by the squared friction velocity orkinematic viscosity, calculated with RANS turbulence models and comparedwith DNS data from Moser [MKM99]. . . . . . . . . . . . . . . . . . . . . . 707.6 (a): Mean streamwise velocity, (b): turbulent kinetic energy and (c): shearstress component of the RST, normalized by the squared friction velocity,calculated with DNS and compared with DNS data from Moser [MKM99]. . 727.7 The different diagonal components of the RST normalized by the squaredfriction velocity, calculated with DNS and compared with DNS data fromMoser [MKM99]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7385 ist of Figures ist of Tables eferences [ARJR15] N. Ahmed, T. Rebollo, V. John, and S. Rubino. A Review of Variational Multi-scale Methods for the Simulation of Turbulent Incompressible Flows.
Archivesof Computational Methods in Engineering , 24:115–164, 12 2015.[Arn82] D. N. Arnold. An interior finite element method with discontinuous elements.
SIAM Journal on Numerical Analysis , 19, 1982.[ARS97] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kuttamethods for time-dependent partial differential equations.
Applied NumericalMathematics , 25, 1997.[BBF13] D. Boffi, F. Brezzi, and M. Fortin.
Mixed Finite Element Methods and Applica-tions , volume 44. 01 2013.[BDM85] F. Brezzi, J. Douglas, and L. Marini. Two Families of Mixed Finite Elementsfor Second Order Elliptic Problems.
Numerische Mathematik , 47:217–235, 061985.[BIL06] L.C. Berselli, T. Iliescu, and W.J. Layton.
Mathematics of Large Eddy Simula-tion of Turbulent Flows . Scientific Computation. Springer, 2006.[BS02] S. Brenner and L.R. Scott.
The Mathematical Theory of Finite Element Methods .Texts in Applied Mathematics. Springer New York, 2002.[Che14] L. Chen. A simple construction of a Fortin operator for the two dimensionalTaylor–Hood element.
Computers and Mathematics with Applications , 68, 112014.[CKS05] B. Cockburn, G. Kanschat, and D. Sch¨otzau. A locally conservative LDGmethod for the incompressible Navier-Stokes equations. mathematics of Com-putation , 74, 2005.[DV06] E. De Villiers. The Potential of Large Eddy Simulation for the Modeling of WallBounded Flows.
Cycle , 01 2006.[Hin75] J.O. Hinze.
Turbulence . McGraw-Hill classic textbook reissue series. McGraw-Hill, 1975.[HMJ00] T.J.R. Hughes, L. Mazzei, and K.E. Jansen. Large eddy simulation and thevariational multiscale method.
Computing and Visualization in Science , 3(1):47–59, May 2000.[HN68] F.H. Harlow and P.I. Nakayana.
Transport of turbulence energy decay rate .LA-3854. 1968. 88 eferences [HT73] P. Hood and C. Taylor. A numerical solution of the Navier-Stokes equationsusing the finite element technique.
Computers and Fluids, Volume 1 , 1973.[JK10] V. John and A. Kindl. A Variational Multiscale Method for Turbulent FlowSimulation with Adaptive Large Scale Space.
J. Comput. Physics , 229:301–312,01 2010.[JKM05] V. John and S. Kaya Merdan. A Finite Element Variational Multiscale Methodfor the Navier–Stokes Equations.
SIAM J. Scientific Computing , 26:1485–1503,01 2005.[Joh04] V. John.
Large Eddy Simulation of Turbulent Incompressible Flows: Analyticaland Numerical Results for a Class of LES Models . Springer-Verlag, 2004.[Kai14] K. Kaiser. Finite Element Methods for the Incompressible Stokes Equationswith Non-Constant Viscosity. Master’s thesis, Freie Universit¨at Berlin, Berlin,2014.[Kuh14] H. Kuhlmann.
Str¨omungsmechanik: Eine kompakte Einf¨uhrung f¨ur Physikerund Ingenieure . Pearson Deutschland GmbH, 2014.[Led16] P. Lederer. Pressure Robust Discretization for Navier-Stokes Equations:Divergence-free Reconstruction for Taylor-Hood Elements and High Order Hy-brid Discontinuous Galerkin Methods. Master’s thesis, Technische Universit¨atWien, Wien, 2016. Diplomarbeit.[Leh10] C. Lehrenfeld. Hybrid discontinuous Galerkin methods for solving incompress-ible flow problems. Master’s thesis, Rheinisch Westfalischen TechnischenHochschule Aachen, Aachen, 2010.[LGS18] P. Lederer, J. Gopalakrishnan, and J. Schoeberl.
A mass conserving mixed stressformulation for the Stokes equations . PhD thesis, Technische Universit¨at Wien,06 2018.[LS16] C. Lehrenfeld and J. Sch¨oberl. High order exactly divergence-free Hybrid Discon-tinuous Galerkin Methods for unsteady incompressible flows.
Computer Methodsin Applied Mechanics and Engineering , 307:339 – 361, 2016.[Men94] F. R. Menter. Two-Equation Eddy-Viscosity Turbulence Models for EngineeringApplications.
AIAA Journal , 32(8):1598–1605, 1994.[MKM99] R. Moser, J. Kim, and N. Mansour. Direct Numerical Simulation of TurbulentChannel Flow up to Re=590.
Physics of Fluids - PHYS FLUIDS , 11:943–945,04 1999.[OBR15] H. Oertel, M. B¨ohle, and T. Reviol.
Str¨omungsmechanik: f¨ur Ingenieure undNaturwissenschaftler . Springer Fachmedien Wiesbaden, 2015.[Pop00] S.B. Pope.
Turbulent Flows . Cambridge University Press, 2000.[PYF15] B. Peng, H. Yan, and H. Fang. Modification of K − ω turbulence modelfor predicting airfoil aerodynamic performance. Journal of Thermal Science ,24:221–228, 2015. 89 eferences [RT06] P. Raviart and J. Thomas.
A Mixed Finite Element Method for Second OrderElliptic Problems , volume 606, pages 292–315. 11 2006.[Sch97] J. Sch¨oberl. NETGEN An advancing front 2D/3D-mesh generator based onabstract rules.
Computing and Visualization in Science , 1:41–52, 07 1997.[Sch07] F.G. Schmitt. About Boussinesq’s turbulent viscosity hypothesis: historicalremarks and a direct evaluation of its validity.
Comptes Rendus M´ecanique ,335(9-10):617–627, 10 2007.[Sch14] J. Sch¨oberl. C++11 Implementation of Finite Elements in NGSolve.
Institute ofAnalysis and Scientific Computing, Vienna University of Technology , 09 2014.[Sch18] J. Sch¨oberl. Numerical Methods for Partial Differential Equations.
Institute ofAnalysis and Scientific Computing, Vienna University of Technology , 10 2018.[Wil94] D.C. Wilcox.
Turbulence Modeling for CFD . DCW Industries, Incorporated,1994.[Zag06] S. Zaglmayr.