Space-time discontinuous Galerkin method for the numerical simulation of the compressible turbulent gas flow through the porous media

Abstract. The article is concerned with the numerical simulation of the compressible turbulent gas flow through the porous media using space-time discontinuous Galerkin method.The mathematical model of flow is represented by the system of non-stationary Reynolds-Averaged Navier-Stokes (RANS) equations. The flow through the porous media is characterized by the loss of momentum. 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 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 the simulation of the Navier-Stokes equations [1][2][3][4][5][6] or compressible turbulent flow ( [7][8][9][10]). This article is devoted to the discretization of viscous compressible turbulent gas flow through the porous media using ST-DG. 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 [11]. The flow through the porous media is characterized by the loss of momentum in the RANS equations. For the numerical simlations will be used the flow through the porous media. These results will be compared with the free stream flow through the dense fence. Computational results show that the ST-DG method is good approach for the solution of these types of problems. The simplified simulation of the porous media was also shown in [12][13][14].
2 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 * e-mail: jan.cessa@seznam.cz * * e-mail: cesenek@vzlu.cz 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 where for s = 1, 2 we have We use the following notation: u = (v 1 , v 2 ) -velocity, ρ -density, p -pressure, θ -absolute temperature, E -total energy, γ -Poisson adiabatic constant, κ -heat conduction coefficient, c v -specific heat at constant volume, c pspecific 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, C 0 is the parameter of the porous media. The above system is completed by the thermodynamical relations 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 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 [11]) for µ T . For the sake of the stability we introduce new variablẽ ω = ln ω. (For details see [7]). Then the turbulence model reads where for s = 1, 2 we havẽ Here ω is the turbulence dissipation, k is the turbulence kinetic energy and µ T is the eddy viscosity. We can write the production term as Limited eddy viscosity ω lim is given by the term The cross-diffusion term C D is defined as The coefficients β, β * , σ k , σ ω , σ D , α ω , C lim , Pr T are chosen by [11]: This system is also equipped with the initial conditioñ w(x, 0) =w 0 (x), x ∈ Ω and the following boundary conditions with given dataw 0 ,ω D , k D ,ω wall . Similary like in the RANS case we can express convect terms and diffusion terms whereK s (w) ∈ I R 2×2 are matrices depending onw.

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 F B h = {Γ ∈ F h ; Γ ⊂ ∂Ω h } and the set of inner 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 and Φ R Γ 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)  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 statew R h is evaluated with the aid of the local linearized Riemann problem described in [5]. For the discretization of the viscous terms we use the property (3) and get the viscous form 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 [6]. They are based on the discontinuity indicator with constants ν 1 and ν 2 . All these forms are linear with respect to w h and nonlinear with respect tow h .
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 I m = (t m−1 , t m ), τ m = t m − t m−1 for m = 1, ..., M. 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 ): ∀Φ hτ ∈ S p,q h,τ . Now we introduce a suitable linearization. We can use two possibilities.
1) We putw hτ (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) 4 Discretization of the k − ω turbulence model

Space discretization of the problem
We apply discontinuous Galerkin method in the similar way as in the previous section. Discretization is carried out on the same mesh for simplicity. An approximate solution of the problem is looked for in the spaceSp h = (Sp h ) 2 and a functionΦ ∈Sp h is, in general, discontinuous on interfaces Γ ∈ F I h .
If we defineP ± = ( ρv · n Γ ) ± I, where I = diag(1, 1), and use the property (5) then we can define convect form where the statew R h is based on the boundary conditions. For the discretization of the diffusion term we use the property (6) and then we havẽ Further we define the turbulent forms h , the interior and boundary penalty formJσ h and the right-hand side forml h in the following way: Boundary statew B is defined on the basis of the Dirichlet boundary conditions and extrapolation.

Full space-time DG discretization
For the full space-time discontinuous Galerkin discretization we use the same partition of the interval [0, T ] as in the previous section. We say that a functionw hτ ∈Sp ,q h,τ is the approximate solution of the problem (4)

Numerical experiments
In this section we will demonstrate the applicability of the developed method and present numerical simulations of the compressible turbulent gas flow through the porous media. Such approach can be used for the simulation of various permeable obstacles. In our case it will be the dense fence (figure 1). We will simulate the fence 0.1m high which will be placed in the free stream. We consider input flow with the total temperature θ 0 = 293K, the total pressure p 0 = 101325Pa and the velocity 20ms −1 and 30ms −1 . The area of the porous media is situated in the same place as is occupied by the fence and we set the parameter  Figure 2 shows detail of the mesh which consist of 346034 elements in total. From the proposed simulation we can guess that porous media is good approximation of the dense fence.
The main advantage of this approach is that we do not have to use so much fine mesh around the real fence.

Conclusion
In this paper we dealt with the space-time discontinuous Galerkin method for the numerical simulation of the compressible turbulent gas flow through the porous media. 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 these types of problems.
This result originated with the support of Ministry of Industry and Trade of the Czech Republic for the long-term strategic