The discontinuous Galerkin method for the numerical simulation of compressible viscous flow

In this paper we deal with numerical simulation of the compressible viscous flow. The mathematical model of flow is represented by the system of non-stationary compressible Navier-Stokes equations. This system of equations is discretized by the discontinuous Galerkin finite element method in space and in time using piecewise polynomial discontinuous approximations. We present some numerical experiments to demonstrate the applicability of the method using own-developed code.


Introduction
During the last decade the discontinuous Galerkin finite element method (DGM), 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.In this paper we present the DGM for numerical simulation of compressible viscous flow around a moving airfoil.A solid airfoil with two degrees of freedom performs rotation around an elastic axis and oscillations in the vertical direction.Compressible viscous flow is described by the system of the Navier-Stokes equations.Due to the time dependent domain the Navier-Stokes equations are transformed in the ALE (Arbitrary Lagrangian-Eulerian) formulation.We consider two approaches to the discretization using the DGM, which differ in the time discretization.First we apply the backward difference formula (BDF) to aproximate the ALE derivative.Second, the full space-time discontinuous Galerkin method (ST-DG) is employed.The discrete flow problem is coupled with the system of ordinary differential equations describing airfoil vibrations.Results of numerical simulation are presented.

ALE formulation of the Navier-Stokes equations
We consider compressible viscous flow in a bounded domain Ω(t) ⊂ I R d , d = 2, 3 , depending on time t ∈ [0, T ].We assume that the boundary Ω(t) consists of three disjoint parts ∂Ω(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 We define the ALE velocity z by the relations a e-mail: cesenek@vzlu.czand the ALE derivative of a vector function w = w(x, t) defined for x ∈ Ω(t) and t ∈ [0, T ]: Then the system of the Navier-Stokes equations describing the compressible viscous flow can be written in the ALE form where for i, j = 1, ..., d we have We use the following notation: u = (v 1 , ..., v d ) -velocity, ρ -density, p -pressure, θ -absolute temperature, E -total energy, γ -Poisson adiabatic constant, c v -specific heat at constant volume, μ > 0, λ = −2μ/3 -viscosity coefficients, k -heat conduction coefficient.The above system is completed by the thermodynamical relations and is equipped with the initial condition w(x, 0) = w 0 (x), x ∈ Ω(0), 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 mapping f s .Expression of the matrices A s can be found in [4].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.Expression of the matricies K s,k (w) can be found in [3].

Space discretization of the flow problem
By Ω h (t) we denote polygonal (d=2) or polyhedral (d=3) 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.In 3D, K ∈ T h (t) can be, e.g., tetrahedrons, prisms or hexahedrons.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 is associated with a unit normal vector n Γ .For Γ ∈ F B h (t) the normal n Γ has the same orientation as the outer normal to ∂Ω h (t).For each Γ ∈ F I h (t) 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 (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 where p > 0 is an integer and P p (K) denotes the space of all polynomials on K of degree ≤ p.
Γ 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 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 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 [4].
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 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 , where d(Γ) is the diameter of Γ ∈ F h (t) 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 solu-tion.In order to overcome this difficulty we employ artificial viscosity forms, see [4].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 .

Time discretization by the BDF method of the flow problem
Let 0 = t 0 < t 1 < ... < t M = T be a partition of the interval [0, T ], and define the time step τ m = t m −t m−1 , m = 1, ..., M. For a time instant t m we use the approximation w(t m ) ≈ w m h , defined in Ω h (t m ).Moreover we set ŵl h = w l h (A t l (A −1 t m (x))), x ∈ Ω h (t m ), l < m.Then we approximate the ALE derivative at time t m by the backward finite difference formula (BDF) of order q where the coefficients α l , l = 0, .., q, depend on τ m−l , l = 0, .., q − 1.In the case m < q we set q := m.For the nonlinear parts of the forms we employ the extrapolation wm where coefficients β l , l = 1, .., q, depend on τ m−l , l = 0, .., q− 1.In the case m < q we set q := m.Expression of the coefficients α l , β l of order q ≤ 3 can be found in [1] or [2].Finally, we set We say that the function w m h ∈ S p h (t m ) is the approximate solution of the problem (1) obtained by the BDF method, EFM 2013 where w 0 h is S p h (t 0 )-approximation of w 0 .

Full space-time DGM discretization
Another way how to construct a method of high-order accuracy both in space and time is the full space-time discontinuous Galerkin (ST-DG) method.We again consider a partition 0 = t 0 < t 1 < ... < t M = T 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 (t 0 ))-projection of w 0 on S p h (t 0 ): Similary as in Section 3.2 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 Fig. 1.The scheme of the vibrating airfoil.

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 i are the coordinates of the elastic axis and τ i j := −pδ i j + τ V i j are the components of the stress tensor.For the derivation of the system (4), see e.g.[5].
The system (4) 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.086622 kg, S a = −0.000779673kg m, EPJ Web of Conferences 02015-p.4 For all numerical simulations we chose symmetric version of the viscous term (SIPG), the parameter C W = 500 for all Γ ∈ F h (t) except for Γ ⊂ Γ W (t), where we set C W = 5000.The constants in the artificial viscosity forms were set ν 1 = ν 2 = 0.1.The time step was choosen τ = 0.003299 c/ |u ∞ | , where |u ∞ | is the far-field velocity.We employed quadratic polynomials (p=2) for the space discretization.In the case of the BDF we used the second order approximation for the time discretization (BDF-p2q2).In the case of the ST-DG we used linear polynomials in time (ST-DG-p2q1).Far field-velocity was choose 20 ms −1 and 40 ms −1 .In this cases Reynolds number is between 10 5 and 10 6 .Figures 2-5 show the computed displacement and rotation of the profile obtained by both methods.At the end we present an example of high-speed flow for far-field velocity 680 ms −1 and Reynolds number 10 7 .In this case we left the data and initial conditions the same except bending stiffness and torsional stiffness which were chosen 1000 times higher.Figures 6-7 show the displacement and rotation of the profile obtained only by ST-DG, because for this regimes the BDF method became unstable.Figures 8 and 9 display distribution of the Mach number and pressure at time instant t = 0.00087s.

Conclusion
In this paper we dealt with the discontinuous Galerkin method for the numerical solution of compressible viscous flow.The applicability of the proposed method was demonstrated on the example of interaction of compressible viscous flow and a moving airfoil.Subject for the further work is including a turbulence model into the method.
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 P(w, n) := d s=1 (A s (w) − z s I)n s , where n = (n 1 , ..., n d ), n 2 1 + ... + n 2 d = 1.Then we have P(w, n)w := d s=1 g s (w)n s .