Finite volume-space-time discontinuous Galerkin method for the numerical simulation of compressible turbulent flow in time dependent domains

The article is concerned with the numerical simulation of the compressible turbulent flow in time dependent domains. The mathematical model of flow is represented by the system of non-stationary ReynoldsAveraged Navier-Stokes (RANS) equations. The motion of the domain occupied by the fluid is taken into account with the aid of the ALE (Arbitrary Lagrangian-Eulerian) formulation of the RANS equations. This RANS system is equipped with two-equation k − ω 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−ω 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.Due to the time dependent domain the RANS equations are transformed in the ALE (Arbitrary Lagrangian-Eulerian) formulation.In this article Wilcox k − ω turbulence model was choosen ( [10]).For numerical experiments will be used flow around moving airfoil.A solid airfoil with two degrees of freedom performs rotation around elastic axis and oscillations in the vertical direction.The discrete flow problem is coupled with the system of ordinary differential equations describing airfoil vibrations.Results of numerical simulation are presented.

Formulation of the k − ω turbulence model
We consider compressible turbulent flow in a bounded domain Ω(t) ⊂ I R 2 , depending on time t ∈ [0, T ].We assume that the boundary Ω(t) consists of three disjoint parts a e-mail: jan.cessa@seznam.cz,cesenek@vzlu.cz∂Ω(t) = Γ I ∪Γ O ∪Γ W (t), where Γ I is the inlet, Γ O is the outlet and Γ W (t) is impermeable wall, whose parts may move.
The time dependence of the domain is taken into account with the aid of a regular one−to−one ALE mapping A t : Ω(0) → Ω(t).
We define the ALE velocity z by the relations z(X, t) = ∂ ∂t A t (X), t ∈ [0, T ], X ∈ Ω(0), and the ALE derivative of a vector function w = w(x, t) defined for x ∈ Ω(t) and t ∈ [0, T ]: Then the system of the RANS equations describing the compressible turbulent flow can be written in the ALE form where for s = 1, 2 we have where 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 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 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 , z 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 : 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 terms as τ T sr ∂v s ∂x r , where Limited eddy viscosity ω 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 .

Space discretization of the flow problem
By Ω h (t) we denote polygonal approximation of the domain Ω(t).Let T h (t) be a partition of the domain Ω h (t) into finite number of closed elements with mutually disjoint interiors such that Ω h (t) = K∈T h (t) K.In 2D problems, we usually choose K ∈ T h (t) as triangles or quadrilaterals.By F h (t) we denote the system of all faces of all elements K ∈ T h (t).Further, we introduce the set of boundary faces 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 (t), 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 4 , S p h (t) = {v; v| K ∈ P p (K), ∀K ∈ T h (t)}, 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 (t) is, in general, discontinuous on interfaces Γ ∈ F I h (t).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 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 (t), integrate over K ∈ T h (t), apply Green's theorem, sum over all elements K ∈ T h (t), 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 (t) 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 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: The boundary state wR h is evaluated with the aid of the local linearized Riemann problem described in [9].
For the discretization of the viscous terms we use the property (4) 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 σ| Γ = C W Re d(Γ) , where d(Γ) is the diameter of Γ ∈ F h (t), 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 forms, see [8].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 , t) := with constants ν 1 and ν 2 .All these forms are linear with respect to w h and nonlinear with respect to wh .

Full space-time DGM 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 ):

Finite volume discretization of the k − ω turbulence model
The finite volume method is applied on the same mesh T h (t) as ST-DG method.So, similarly as in Section 3.1, let T h (t) = {K i (t)} i∈I be a partition of the domain Ω h (t) into finite number of closed elements with mutually disjoint interiors such that Ω h (t) = i∈I K i (t), 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 The boundary ∂Ω h (t) is formed by a finite number of sides of elements K i (t) adjacent to ∂Ω h (t).We denote all these boundary sides by S j (t), where For an element K i (t), not containing any boundary side S j (t), we set γ(i) = ∅.Moreover we define S (i) = s(i) ∪ γ(i) and n i j (t) = ((n i j (t)) 1 , (n i j (t)) 2 ) as the unit outer normal to ∂K i (t) on the side Γ i j (t).We again consider a partition 0 = t 0 < t 1 < ... < t M = T of the interval [0, T ] and denote Now we can integrate the system (5), apply Green's theorem and approximate the integrals with the one-point rule and we get the following implicit finite volume scheme The scheme of the vibrating airfoil.
+ j∈S (i) Rs ( wm 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 w(x, t m ) dx.
For the numerical flux H we again use the concept of the Vijayasundaram numerical flux in the following way: • n i j ) − wB for j ∈ γ(i), 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].

Equations for the moving airfoil
We shall simulate the motion of a profile in 2D with two degrees of freedom: H -displacement of the profile in the vertical direction (positively oriented downwards) and αthe rotation of the profile around the so-called elastic axis (positively oriented clockwise), see Figure 1.The motion of the profile is described by the system of ordinary differential equations where we use the following notation: m -mass of the airfoil, L(t) -aerodynamic lift force, M(t) -aerodynamic torsional moment, S α -static moment of the airfoil around the elastic axis (EA), k HH -bending stiffness, k αα -torsional stiffness.We define L(t) and M(t) by the terms where l is the depth of the airfoil, x EA s are the coordinates of the elastic axis and τ sr := −pδ sr +τ V sr are the components of the stress tensor.For the derivation of the system (9), see e.g.[11].
The system ( 9) is transformed to a first-order system and solved by the fourth-order Runge-Kutta method together with the discrete flow problem.

Numerical experiments
We performed numerical simulations in 2D for the profile NACA0012 with the following data and initial conditions: m = 0.
In the case of space-time discontinuous Galerkin method we employ linear polynomials in space and in time.Far field-velocity was chosen in the range 10 m s −1 to 40 m s −1 .Figures 2-13 show the computed displacement and rotation of the profile.Figures 14 -17 show different quantities for far-field velocity 37.5 m s −1 at the time instant 0.34 s. Figure 18 shows detail of the mesh which consist of 69766 elements in total.

Conclusion
In this paper we dealt with the finite-volume -space-time discontinuous Galerkin method for the numerical solution

Figure 14 .
Figure 14.The distribution of the pressure.

Figure 15 .Figure 16 .
Figure 15.The distribution of the velocity.

Figure 17 .
Figure 17.The distribution of the Mach number.

Figure 18 .
Figure 18.Detail of the mesh.