Finite volume-space-time discontinuous Galerkin method for the solution of compressible turbulent flow

In this article we deal with numerical simulation of the non-stationary compressible turbulent flow. Compressible turbulent flow is described by the Reynolds-Averaged Navier-Stokes (RANS) equations. This RANS system is equipped with two-equation k-omega turbulence model. These two systems of equations are solved separately. Discretization of the RANS system is carried out by the space-time discontinuous Galerkin method which is based on piecewise polynomial discontinuous approximation of the sought solution in space and in time. Discretization of the two-equation k-omega turbulence model is carried out by the implicit finite volume method, which is based on piecewise constant approximation of the sought solution. We present some numerical experiments to demonstrate the applicability of the method using own-developed code.


Introduction
During the last decade the space-time discontinuous Galerkin finite element method (ST-DG), which is based on piecewise polynomial discontinuous approximations of the sought solution, became very popular in the field of numerical simulation of the fluid flow.This method of higher order was successfully used for simulation of the Navier-Stokes equations ( [2], [4], [5]).Unfortunately this method became unstable, when it was used for the whole system of the RANS equations with k − ω equations.Thus we solve this two systems separately.We use the ST-DG for the discretization of the RANS system of equations and the finite volume method for the equations of the k − ω turbulence model.As numerical experiments show this is good approach for the solution of this problem.

Formulation of the k − ω turbulence model
We consider compressible turbulent flow in a bounded domain Ω ⊂ I R d , d = 2, 3. We assume that the boundary of Ω consists of three disjoint parts: ∂Ω = Γ I ∪ Γ O ∪ Γ W , where Γ I is the inlet, Γ O is the outlet and Γ W is the impermeable wall.
The system of the RANS equations can be written in the following form where for s = 1, ..., d we have w = (w 1 , ..., w d+2 ) T = (ρ, ρv 1 , ..., ρv d , E) T ∈ I R d+2 , a e-mail: jan.cessa@seznam.cz,cesenek@vzlu.cz where We use the following notation: u = (v 1 , ..., v d ) -velocity, ρ -density, p -pressure, θ -absolute temperature, E -total energy, γ -Poisson adiabatic constant, κ -heat conduction coefficient, c v -specific heat at constant volume, c p -specific heat at constant pressure, where c p = γc v , Pr is the laminar Prandtl number, which can be express in the form Pr = c p μ L κ , Pr T is the turbulent Prandtl constant number, μ T is the eddy-viscosity coefficient, μ L is the dynamic viscosity coefficient dependent on temperature via Sutherland's formula.The above system is completed by the thermodynamical relations with given data w 0 , ρ D , u D .It is possible to show that f s (αw) = α f s (w) for α > 0. This property implies that where are the Jacobi matrices of the mappings f s .The viscous terms R s (w, ∇w) can be expressed in the form where K s,k (w) ∈ I R (d+2)×(d+2) are matrices depending on w.
The aforesaid system of the RANS equations has to be completed by a turbulence model for μ T : where for s = 1, ..., d we have w = ( w1 , w2 ) T = (ρω, ρk) Here ω is the turbulence dissipation, k is the turbulence kinetic energy and μ T = ρk ω is the eddy viscosity.We can write the production terms as The cross-diffusion term C D is defined as The coefficients β, β * , σ k , σ ω , σ D , α ω are chosen by [9]: This system is also equipped with the initial condition w(x, 0) = w0 (x), x ∈ Ω and the following boundary conditions a) with given data w0 , ω D , k D , ω wall .
3 Space-time discretization of the RANS equations

Space discretization
By Ω h we denote a polygonal (d=2) or polyhedral (d=3) approximation of the domain Ω.Let T h be a partition of the domain Ω h into finite number of closed elements with mutually disjoint interiors such that Ω h = K∈T h K.In 2D problems, we usually choose K ∈ T h as triangles or quadrilaterals.In 3D, K ∈ T h can be, e.g., tetrahedrons, prisms or hexahedrons.By F h we denote the system of all faces of all elements K ∈ T h .Further, we introduce the set of boundary faces h the normal n Γ has the same orientation as the outer normal to ∂Ω h .For each Γ ∈ F I h there exist two neighbouring elements We use the convention that K R Γ lies in the direction of n Γ and K L Γ lies in the opposite direction to n Γ .If Γ ∈ F B h , then the element adjacent to Γ will be denoted by K L Γ .Because of the space discretization we define the space of piecewise polynomial functions where p > 0 is an integer and P p (K) denotes the space of all polynomials on K of degree ≤ p.A function Φ ∈ S p h is, in general, discontinuous on interfaces Γ ∈ F I h .By Φ L Γ and Φ R Γ we denote the values of Φ on Γ considered from the interior and the exterior of K L Γ , respectively, and set which denotes average and jump of Φ on Γ.
The discrete problem is derived in the following way: For arbitrary t ∈ [0, T ] we can multiply the system (1) by a test function EFM 2015 theorem, sum over all elements K ∈ T h , use the concept of the numerical flux and introduce suitable terms mutually vanishing for a regular exact solution.Moreover, we carry out a suitable partial linearization of nonlinear terms.
In order to evaluate the integrals over Γ ∈ F h in inviscid term we use the approximation where H is a numerical flux.For the construction of the numerical flux we use the properties (3) of f s .Let us define the matrix where n = (n 1 , ..., n d ), n 2 1 + ... + n 2 d = 1.Then we have It is possible to show that the matrix P is diagonalizable.It means that there exists a nonsingular matrix T = T(w, n) and a diagonal matrix Λ = Λ(w, n) such that where λ i = λ i (w, n) are eigenvalues of the matrix P. Then we can define the "positive" and "negative" parts of the matrix P by where λ + = max(λ, 0), λ − = min(λ, 0).Using this concept, we introduce the so-called Vijayasundaram numerical flux This numerical flux has suitable form for a linearization.Now we can define inviscid form in the following way: The boundary state wR h is evaluated with the aid of the local linearized Riemann problem described in [8].
For the discretization of the viscous terms we use the property (4) and get the viscous form a h ( wh , w h , Φ h ) := We set Θ = 1 or Θ = 0 or Θ = −1 and get the so-called symmetric version (SIPG) or incomplete version (IIPG) or nonsymetric version (NIPG), respectively, of the discretization of viscous terms.Further, we define the interior and boundary penalty form J σ h and the right-hand side form l h in the following way: where σ| Γ = C W Re d(Γ) , where d(Γ) is the diameter of Γ ∈ F h , Re is the Reynolds number and C W > 0 is a suitable sufficiently large constant.The boundary state w B is defined on the basis of the Dirichlet boundary conditions and extrapolation: In the vicinity of discontinuities or steep gradients nonphysical oscillations can appear in the approximate solution.In order to overcome this difficulty we employ artificial viscosity form, see [7].They are based on the discontinuity indicator where [ ρh ] is the jump of the function ρh (= the first component of the vector function wh ) on the boundary ∂K, d(K) denotes the diameter of K and |K| denotes the area of the element K. Then we define the discrete discontinuity indicator and the artificial viscosity forms βh ( wh , w h , Φ h ) := with constants ν 1 and ν 2 .All these forms are linear with respect to w h and nonlinear with respect to wh .
Finally, we set

Full space-time DGM discretization
Let 0 = t 0 < t 1 < ... < t M = T be a partition of the interval [0, T ] and denote We define the space S p,q h,τ = (S p,q h,τ ) d+2 , where with integers p, q ≥ 1. P q (I m ) denotes the space of all polynomials in t on I m of degree ≤ q.Moreover for Φ ∈ S p,q h,τ we introduce the following notation: Approximate solution w hτ of the problem will be sought in the space S p,q h,τ .Since the functions of this space are in general discontinuous in time, we ensure the connection between I m−1 and I m by the penalty term in time The initial state w hτ is included by the L 2 (Ω h )-projection of w 0 on S p h : Now we introduce a suitable linearization.We can use two possibilities.1) We put whτ (t) := w h (t − m−1 ) for t ∈ I m .2) We prolong the solution from the time interval I m−1 to the time interval I m .
We say that a function w hτ ∈ S p,q h,τ is the approximate solution of the problem (1) obtained by the ST-DG method, if it satisfies the following conditions 4 Finite volume discretization of the For an element K i , not containing any boundary side S j , we set γ(i) = ∅.Moreover we define S (i) = s(i) ∪ γ(i) and n i j = ((n i j ) 1 , ..., (n i j ) d ) as the unit outer normal to ∂K i on the side Γ i j .We again consider a partition 0 = t 0 < t 1 < ... < t M = T of the interval [0, T ] and denote I m = (t m−1 , t m ), τ m = t m − t m−1 for m = 1, ..., M. Now we can integrate the system (5) over set K i × (t m−1 , t m ), apply Green's theorem and approximate the integrals with the one-point rule and we get the following implicit finite volume scheme, Here wm Γ i j resp.∇ wm Γ i j denote wm resp.∇ wm at the center of the edge Γ i j at time instant t m .By wm i we shall denote approximate solution of w on the element K i at time instant t m For the numerical flux H we again use the concept of the Vijayasundaram numerical flux in the following way , where wB denotes boundary conditions.One possibility of linearization of Rs and s is via Taylor expansion which is used in our case.For more details see [6].

Numerical experiments
In order to demonstrate the applicability of the developed method we shall present two numerical simulations of compressible turbulent flow in 2D.For this purpose we chose two regimes for the profile RAE2822 and compared it with experimental data, which can be found in [1].First regime is characterized by the far-field Mach number M ∞ = 0.6 and the angle of the attack α = 2.57 o ('case 3' in [1]) and second case is characterized by the Mach number M ∞ = 0.73 and the angle of the attack α = 3.19 o ('case 9' in [1]).Figures 1 and 2 show the pressure distribution.Comparison of the computed pressure coefficient c p with the experimental data is shown in Figures 3 and 4.

Conclusion
In this paper we dealt with the finite-volume -space-time discontinuous Galerkin method for the numerical solution of compressible turbulent flow.The applicability of the proposed method was demonstrated on the examples, which show, that the presented method is an efficient numerical scheme for the solution of the 2D compressible turbulent flow.

Acknowledgment
This result originated with the support of Ministry of Industry and Trade of the Czech Republic for the long-term

/
k and is equipped with the initial condition w(x, 0) = w 0 (x), x ∈ Ω and the following boundary conditions a) ρ| Γ I = ρ D , DOI: 10.1051/ C Owned by the authors, published by EDP Sciences, 201

Figure 1 .
Figure 1.The pressure distribution for the case 3.

Figure 2 .
Figure 2. The pressure distribution for the case 9.

Figure 3 .
Figure 3.Comparison of the computed pressure coefficient c p with the experimental data for the case 3.

Figure 4 .
Figure 4. Comparison of the computed pressure coefficient c p with the experimental data for the case 9.
two-equation k − ω turbulence modelSimilarly as in Section 3.1 let T h = {K i } i∈I be a partition of the domain Ω h into finite number of closed elements with mutually disjoint interiors such that Ω h = i∈I K i , where I ⊂ Z + = {0, 1, 2, 3, ...} is a suitable index set.If two elements have a common face, then we call them neighbours and put Γ i j = Γ ji = ∂K i ∩ ∂K j .For each i ∈ I we define the index set s(i) = j ∈ I; K j is a neighbour of K i .The boundary ∂Ω h is formed by a finite number of sides of elements K i adjacent to ∂Ω h .We denote all these boundary sides by S j , where j ∈ I b ⊂ Z − = {−1, −2, −3, ...} and set γ