Space-time discontinuous Galerkin method for the numerical simulation of viscous compressible gas flow with the k-omega turbulence model

In this article we deal with the numerical simulation of the non-stationary compressible turbulent flow described by the Reynolds-Averaged Navier-Stokes (RANS) equations. This RANS system is equipped with two-equation k-omega turbulence model. The discretization of these two systems is carried out separately by the space-time discontinuous Galerkin method. This method is based on the piecewise polynomial discontinuous approximation of the sought solution in space and in time. We use the numerical experiments to demonstrate the applicability of the shown approach. All presented results were computed with the 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 the simulation of the Navier-Stokes equations [1][2][3][4][5][6].In the case of compressible turbulent flow the finite volume -space-time discontinuous Galerkin method was used [7][8][9], where the equations of the turbulence model were discretized by the implicit finite volume method.This article is devoted to the discretization of viscous compressible turbulent gas flow using ST-DG applied also for the equations of the turbulence model.The flow is described by the system of the RANS equations equipped with the system of k − ω equations.These two systems are solved separately and discretized by ST-DG.In this article Wilcox k −ω turbulence model was chosen [10].The flow around airfoil RAE2822 will be used for the numerical simulations.These results are compared with experiments.Computational results show that the ST-DG method is good approach for the solution of these types of problems.

Formulation of the k − ω turbulence model
We consider compressible turbulent flow in a bounded domain Ω ⊂ I R 2 .We assume that the boundary Ω consists of three disjoint parts ∂Ω = Γ I ∪ Γ O ∪ Γ W , where Γ I is the inlet, Γ O is the outlet and Γ W is impermeable wall.
The system of the RANS equations describing the viscous compressible turbulent flow can be written in the form ∂f p s (w) ∂x s a e-mail: jan.cessa@seznam.cz,cesenek@vzlu.cz where for s = 1, 2 we have We use the following notation: and is equipped with the initial condition w(x, 0) = w 0 (x), x ∈ Ω and the following boundary conditions a) 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 .Similary we can express f p s (w) = A p s (w)w, s = 1, 2. The viscous terms R s (w, ∇w) can be expressed in the form where K s,k (w) ∈ I R 4×4 are matrices depending on w.
The aforesaid system of the RANS equations is completed by the Wilcox's turbulence model (see [10]) for µ T .For the sake of the stability we introduce new variable ω = ln ω. (For details see [11]).Then the turbulence model reads where for s = 1, 2 we have Here ω is the turbulence dissipation, k is the turbulence kinetic energy and µ T is the eddy viscosity.We can write the production term as where Limited eddy viscosity ω lim is given by the term τrs τrs The cross-diffusion term C D is defined as The coefficients β, β * , σ k , σ ω , σ D , α ω , C lim , Pr T are chosen by [10]: 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 .Similary like in the RANS case we can express convect terms and diffusion terms where Ks (w) ∈ I R 2×2 are matrices depending on w.

Space discretization of the problem
By Ω h we denote polygonal 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.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 Γ .We shall look for an approximate solution of the problem in 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 Γ we denote the values of Φ on Γ considered from the interior and the exterior of K L Γ , respectively, and set which denotes the average and jump of Φ on Γ.
The discrete problem is derived in the following way: For arbitrary t ∈ [0, T ] we can multiply the system by a test function S p h , integrate over K ∈ T h , apply Green's 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 and then we can define the following forms.
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 (2) of f s .Let us define the matrix where n = (n 1 , n 2 ), n 2 1 + n 2 2 = 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 where the boundary state wR h is evaluated with the aid of the local linearized Riemann problem described in [1].For the discretization of the viscous terms we use the property (3) 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 the viscous terms.Further, we define the turbulent forms p h , k h , the interior and boundary penalty form J σ h and the right-hand side form l h in the following way: where σ is a parameter of the method and 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 forms, see [3].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 DG discretization
Let 0 = t 0 < t 1 < ... < t M = T be a partition of the interval [0, T ] and let us denote We define the space S p,q h,τ = (S p,q h,τ ) 4 , 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 (t 0 ))-projection of w 0 on S p h (t 0 ): 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) Further we define the turbulent form sh , the interior and boundary penalty form J σ h and the right-hand side form lh in the following way: Boundary state wB is defined on the basis of the Dirichlet boundary conditions and extrapolation.

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 [12].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 [12]) and second case is characterized by the Mach number M ∞ = 0.73 and the angle of the attack α = 3.19 o ('case 9' in [12]).For both cases we set free-stream turbulence intensity T u = 1%.

Conclusion
In this paper we dealt with the space-time discontinuous Galerkin method for the numerical solution of the viscous 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 viscous compressible turbulent flow.
show the pressure distribution, Mach number distribution and turbulent kinetic energy distribution.Comparison of the computed pressure coefficient c p with the experimental data is shown in Figures 4 and 8. Figure9shows detail of the mesh which consist of 64150 elements in total.

Fig. 1 .
Fig. 1.The distribution of the pressure for the case 3.

Fig. 2 .
Fig. 2. The distribution of the Mach number for the case 3.

Fig. 3 .
Fig. 3.The distribution of the turbulent kinetic energy for the case 3.

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

Fig. 5 .
Fig. 5.The distribution of the pressure for the case 9.

Fig. 6 .
Fig. 6.The distribution of the Mach number for the case 9.

Fig. 7 .
Fig. 7.The distribution of the turbulent kinetic energy for the case 9.

Fig. 8 .
Fig. 8.Comparison of the computed pressure coefficient c p with the experimental data for the case 9.