Numerical simulation of fluid flow through simplified blade cascade with prescribed harmonic motion using discontinuous Galerkin method

This paper deals with a numerical simulation of compressible viscous fluid flow around three flat plates with prescribed harmonic motion. This arrangement presents a simplified blade cascade with forward wave motion. The aim of this simulation is to determine the aerodynamic forces acting on the flat plates. The mathematical model describing this problem is formed by Favre-averaged system of Navier-Stokes equations in arbitrary Lagrangian-Eulerian (ALE) formulation completed by one-equation Spalart-Allmaras turbulence model. The simulation was performed using the developed in-house CFD software based on discontinuous Galerkin method, which offers high order of accuracy.


Introduction
In the forties of the 20th century, first papers focused on the aero-elasticity of blades in axial turbines appeared [1]. Empirical criterion on arising of flutter is introduced in the paper [2]. This criterion is based on a reduced blade frequency for bending blade vibration. If the reduced blade frequency is less than 0.33, the flutter is not excited. Nevertheless, other influences such as geometric configuration of blade cascade and angle of attack can also lead to flutter.
The problem of aero-elasticity is still subject of topical interest because modern turbines are designed for higher efficiencies and higher power under higher operational temperatures and flow rates. Higher operational, safety and economical demands force the designers to be more precise during the process of determination of safe operational condition laying out of the area with loss of stability. Originally, the term flutter was established for flow attached to the vibration of blade when phase displacement between blade motion and induced aerodynamic forces occurs. In dependence on time shift between the moving structure and unsteady aerodynamic forces acting on the structure, following cases can occur: • the energy of the flow is absorbed into the blade, • kinetic energy of the blade is transmitted to the flow, • the energy of the flow stays unchanged during particular cycles of blade vibration.
If the blades add the energy to the flow, we talk about damped vibration from the aero-elasticity point of view. If the flow adds the energy to blades and their amplitudes increase, we talk about unstable excited vibration. In the work [3] it has been shown that the critical flutter point, For purpose of this study, we choose the travelling wave mode approach for the numerical modelling of the unsteady flow and flutter analysis in the simplified blade cascade formed by three flat plates. Thus, all flat plates in the cascade vibrate and their mutual movement prescribes the travelling waves in the cascade. The aim of this numerical study is to determine aerodynamic forces acting on the flat plates in cascade, when the harmonic motion with phase delay is prescribed. The numerical simulation is performed using developed in-house CFD software based on the discontinuous Galerkin finite element method. In section 4 we present numerical results of unsteady fluid flow for the selected configuration of the simplified blade cascade described in [3]. For the numerical simulation one flow speed and forward travelling wave mode in the cascade were considered.

Mathematical model
The two-dimensional problem of compressible viscous fluid flow through simplified blade cascade with prescribed harmonic motion is described by Favre-averaged system of Navier-Stokes equations, in dimensionless arbitrary Lagrangian-Eulerian (ALE) formulation written as where i, j = 1, 2,¯ andp are dimensionless time-averaged values of density and pressure,ũ i andẽ are dimensionless mass-averaged velocity components and energy, Re is Reynolds number, Pr is Prandtl number, Pr t = 0.89 is turbulent Prandtl number for flow around flat plate, µ and µ t are dynamic and turbulent viscosity. U j are components of mesh velocity. Mass-averaged stress tensor and Reynolds stress tensor are given bỹ The equation of state for ideal gas is considered as The symbol D A Dt denotes the ALE derivative [5]. To include the influence of turbulent fluctuations on the mean flow a one-equation turbulence model of Spalart and Allmaras [4] is used. It is given by transport equation for eddy viscosityν written as and completed by the following relations where D is the distance to the nearest wall andΩ is the vorticity magnitude.
The system of Favre-averaged Navier-Stokes equations (1) and the transport equation for eddy viscosity (2) are numerically solved together, the whole mathematical model can be written in a dimensionless compact flux vector form as where w is a vector of flow variables, f j (w) are the inviscid fluxes, f v j (w, ∇w) are the viscid fluxes and p(w, ∇w) is a production term. The vectors are then written as
and P q (Ω k ) is the space of polynomials of q-th order at the most.
Multiplying the system of equations (3) by a test function v ∈ S h , integrating it over the element Ω k and using the Green's theorem, the following integral identity is obtained The key part of the method is an approximation of curve integrals. They are evaluated in the same way as it is done in the case of finite volume method. The inviscid fluxes are evaluated as where H is Lax-Friedrichs flux, [9], n = [n 1 , n 2 ] T is vector of outer normal, values w − resp. w + are the left resp. right limits of the vector w on the edge Γ int ∪ Γ b and w b is a boundary value of w computed from boundary conditions. Similarly, the value of viscous fluxes are evaluated using the appropriate numerical fluxes. In present paper, the approximation known as interior penalty (IP) method was chosen, [8].
Let w k be a part of solution w corresponding to control element Ω k and ϕ k,i be i-th basis function of the space P q (Ω k ). The m components w m k , m = 1, 2, ..., 4, of vector w k can be expressed as a linear combination Inserting it into the integral identity and evaluating the volume and curve integrals using the Gauss integration rules of appropriate order, the following system of ordinary differential equations is obtained where vector W consists of coefficients of linear combination (4) and M is a mass matrix. In this study, linear Lagrange polynomials are chosen as the basis functions of the space P q (Ω k ). Time integration of equation (5) is performed using the implicit backward Euler method, which is solved by the Newton's method The system of linear equations (6) is solved iteratively using the GMRES method with block diagonal Jacobi preconditioner.
For the computation of deformed mesh points coordinates, the blending function approach was used [10], [11]. For our purpose, the blending function was generalized for three independently moving blades. Let x 0 denotes the vector of initial coordinates of each mesh node. For each blade (i = 1, 2, ..., n, in our study n = 3), the vector of initial coordinates x 0 is transformed into a new vector of coordinates y i using the following formula: where p i = p i (t) is the time dependent translation vector and T(α) i = T(α(t)) i is the time dependent rotation matrix of i-th blade. The coordinates x = x(t) of the new mesh points are computed by applying the blending functions as The key part of this algorithm is to determine the blending functions b i . The blending function b i have to satisfy two basic properties:

Numerical results
The developed in-house CFD software based on the discontinuous Galerkin method was applied to the numerical simulation of compressible viscous fluid flow through a simplified blade cascade with prescribed harmonic motion. The geometric configuration of the simplified blade cascade and the angle of attack of the flow were adopted from the test case no. 3, which was experimentally studied in [3]. The only difference is that the thin airfoils of test cascade in [3] were replaced by 3 millimetres thick flat plates. The geometric configuration of the simplified test cascade is shown in Figure 1. The considered two-dimensional computational domain and corresponding unstructured triangular mesh with 11932 elements is shown in Figure 2. The detailed view of the mesh in the vicinity of the three flat plates is shown in Figure 3. Each plate performs kinematic harmonic motion in vertical direction y = A sin(2π f t + φ) with an interblade phase angle φ = π/2 for the upper plate, φ = 0 for the middle plate and φ = −π/2 for the bottom plate. The amplitude and the frequency of the harmonic motion are A = 0.003 m and f = 20 Hz, respectively. At the inlet of the computational domain, see Figure 2, the stagnation pressure p 01 = 101325 Pa, the stagnation temperature T 01 = 293.15 K and angle of attack α = 11.4 were prescribed. At the outlet, the static pressure p 2 = 95520 Pa was prescribed. The value of p 2 was calculated out of the free stream velocity 100 ms −1 under the consideration of standard sea level air conditions.
The pressure contours are shown at three selected times t 1 = 0.0433 s, t 2 = 0.1109 s and t 3 = 0.1345 s in Figures 4-6. The times t 1 , t 2 and t 3 correspond respectively to the zero, maximal and minimal vertical displacement of the middle flat plate. The time evolution of y-component of aerodynamic forces acting on upper, middle and bottom plates is shown in Figure 7. The results produced by the developed in-house solver were compared to the results obtained by the commercial CFD software package ANSYS Fluent. The y-components of aerodynamic forces acting on upper, middle and bottom plates obtained by both the in-house and commercial solvers are shown in Figure 8. Finally, vertical displacement of the middle plate along with y-component of the aerodynamic force acting on the middle plate are plotted in Figure 9.

Conclusion
The in-house CFD solver based on the discontinuous Galerkin method was developed for the numerical solution of compressible viscous flow problems in domains with moving boundaries. This solver was applied to the numerical simulation of compressible viscous fluid flow through the simplified blade cascade formed by three flat plates with prescribed harmonic motion.
The paper presents numerical results in the form of pressure distribution and time evolution of aerodynamic forces acting on each flat plate for the free stream velocity 100 ms −1 and for the forward travelling wave mode, i.e. the inter-blade phase angle φ = π/2. The y-components of aerodynamic forces acting on the flat plates obtained by the developed in-house CFD solver show qualitative agreement with those from the commercial CFD software package ANSYS Fluent. Comparing the time evolution of the middle plate vertical displacement with the time evolution of the y-component of aerodynamic force acting on the middle plate, one can observe a small time delay of the force. Such delay is typical for the case of forward wave motion through cascade when kinetic energy of the plate is transmitted to the flow.
The objective of our future study will be aimed at numerical computation of the y-components of aerodynamic forces acting on blades and their time delays to the blade's vertical displacements for different travelling waves through the blade cascade in order to assess the flutter stability of the cascade in the whole range of the inter-blade phase angle φ from 0 to 2π. Obtained numerical results will be experimentally validated in the Institute of Thermomechanics of the Czech Academy of Sciences.