Numerical Investigation of Aeroelastic Flutter in Two-Dimensional Cascade of Compressor Blades

Following contribution presents numerical study of aeroelastic flutter in two-dimensional section of flat wing cascade in wind tunnel. The investigation is conducted as a parametric study of varying pitch angle of one (middle) blade in the cascade with each computational case performed on fixed computational grid. This approach can be viewed as an approximation of fluid-structure interaction realized on moving mesh. Numerical predictions were carried by means of CFD open-source codes OpenFOAM® and Nektar++. The particular aim was focused on assessment of numerical performance and accuracy of the numerical solvers as well as several turbulence models.


Introduction
There are two major areas of applied aerodynamics research where flutter occurrence plays important role. The first one is aircraft industry with strict requirements of safety resulting in careful wing design preventing any kind of dynamic instability. The second one is gas turbine engines where jet-engine companies keep pushing the limits of efficiency and performance. However, internal high-speed flow through usually complex geometry in gas turbines makes it complicated to create aerodynamic functions used for aircraft flutter prediction. Therefore mathematical models of fluid-structure-interaction consisting of strongly coupled systems of equations describing time-dependent dynamics of both flow field and deformable solid structure are often chosen. The aim of this work is to asses performance and accuracy of numerical solvers as well as turbulent models with usage of simplified physical assumption of that time-dependent interaction of air flow and moving blade can be approximated by series of steady-state cases with varying position of fluttering blade as a parameter.

Mathematical model
As numerical results were obtained by means of two open-source codes OpenFOAM® [1] and Nektar++ [2] brief description of mathematical and turbulence models used in both codes is provided in this section.

OpenFOAM
The mathematical model of compressible flow is formed by the Favre-averaged Navier-Stokes equations. The system of equations is completed by constitutional relations for the ideal gas and two-equation SST turbulence model [3] with heat transfer modelled as constant Prandtl number. Two variants of SST model were used for comparison: standard model, kω-SST, and three equation model with bypass transition, γ-SST [4]. Furthermore, two numerical solvers were used for computation. The native OpenFOAM solver, rhoSimpleFoam, and coupled matrix-free LU-SGS solver developed by Fürst [5].

Nektar++
Simulation based on solution of Incompressible Navier-Stokes equations was the second approach. In this case, the problem was considered as evolutionary and without any turbulence model.
The IncNavierStokesSolver from Nektar++ library [2] was used in this case. It implements velocity-correction scheme [6] for discretization in time and allows high degree polynomial approximations in space (spectral/hp element method). The spectral/hp approach was chosen as it is promising in sense of fast/spectral convergence in comparison to finite volumes, see e.g. [7]. High order approximation in space lowers the numerical diffusion and dispersion. The spectral/hp method is used with intention to achieve higher accuracy. However, due to relatively high speed of the flow and limited computational resources we left ambition to achieve full spectral convergence which indicate high-precision solutions. Stabilization techniques of Spectral Vanishing Viscosity were used instead [8], [9].
As the problem was considered as evolutionary and incompressible it stays in contrast to approaches mentioned in previous paragraph.

Configuration and flow conditions
All the numerical simulations were executed using the same geometric configuration see Figure 1. The computational domain consists of 2D channel with blade cascade of 5 blades. Each blade is hereafter denoted as blade_0x with x varying from 1 to 5 as going from the bottom up. The blade profile is designed as flat plate of chord length c = 120 mm and thickness t = 5 mm with circular leading and trailing edge. The slope of cascade is determined by β = 31.5°. Width and height of computational domain is set as h ~ 195mm and w = 1000 mm. The pitch of the cascade is given by s = 74.52 mm. Parametric study was realized by changing the incidence angle α of blade_03 from α = -5° to α = +5° with an increment of 1°. In the case of Nektar++ the fluid was considered as incompressible and modeled as unsteady. The flow regime was set to keep similarity to the compressible model in sense of the Reynolds number. The setting followed from above values and the reference values, Tref = 274.43 K, Uref = 166.05 m•s -1 , ρref = 1.0702 kg•m -3 . Boundary conditions for velocity were set as "no-slip" on blade and channel walls and constant profile on inlet. On outflow, the stabilized boundary condition of Dirichlet type for pressure and Neumann type for velocity was prescribed with default setting provided from Nektar++. High order pressure boundary condition proposed in Karniadakis [6] was used on the no-slip walls and inlet.

Computational meshes
In the case of OpenFOAM the computational mesh contains about 435 000 cells of rectangular and triangular shape in free stream and rectangular shape in boundary layer. The nearest node is in the distance from the wall about y+ ≈ 5 in average. The case of incompressible flow solution was discretized in spatial coordinates on quadrilateral mesh of approximatively 19000 elements and uniform polynomial degree 6 for velocity components and degree 5 for pressure. Curved boundaries were represented by high-order deformed elements. Elements neighboring to boundaries with no-slip condition were narrowed in direction perpendicular to the wall.

OpenFOAM
Regarding to the rhoSimpleFoam, SIMPLE algorithm [10] based, steady-state solver for compressible fluids following numerical setting was used. The native limited central-upwind differencing scheme of 2 nd order was used for discretisation of both inviscid and viscous fluxes. The iterative process employed by SIMPLE was used to reach steady-state.
In the matter of the lower-upper symmetric Gauss-Seidel (LU-SGS) matrix-free compressible solver the inviscid fluxes were discretized with use of limited piece-wise linear reconstructions with the HLLC numerical flux [11] and viscous terms are discretized with the OpenFOAM build-in central differencing schemes. The local pseudo-time marching method based on 1 st order backward difference was utilized to reach steady-state.

Nektar++
Discretization in time employed the widely used implicit-explicit scheme and applied order was 2. Timestep 10 -5 was safely under limit for stability while using stabilization techniques of Spectral Vanishing Viscosity (SVV). The SVV technique was used accordingly to relatively high Reynolds number. Presently, Nektar++ offers three types of SVV kernel -Exponential, Power and DG. All these three versions introduce internal parameters, which affect its impact to accuracy of the solution and solver stability. Full resolution was beyond the scope of present study. Testing simulations for various kernels (and cut-off parameter) were performed to compare influence of SVV to results. Data arising from simulations with Exponential kernel and default parameters are presented in the following section.

Numerical results
This section presents some of the numerical results obtained by means of OpenFOAM and Nektar++ solvers. The attention is paid to the assessment of the skin friction coefficient which well documents differences between tested models of turbulence with regards to the prediction of flow separation and reattachment on blade profiles. Furthermore, the aerodynamic forces in a form of dimensionless coefficients are discussed. The comparison of those quantities is then made among finite volume approach as well as spectral element method as a contrast.
In the case of OpenFOAM it has been found that the LU-SGS solver provides much more robust approach than rhoSimpleFoam when dealing with higher Mach numbers with respect to the setup of the computational case. This is especially apparent when tuning up relaxation factors, initial conditions and bounding constraints in SIMPLE algorithm in order to find the parameters leading to the convergence of the case. This is well documented also in [5].
Nektar++ was set to provide both the instant and time-averaged fields of velocity and pressure. These were used for evaluation of friction coefficient distribution, see Figure 4 Other traced quantities, aerodynamic forces and moment, were evaluated at every 10 th timestep from the instant fields. The time-dependent data were statistically processed to obtain the final average values and dispersion, what is in contrast with data obtained using the finite volume method with turbulence models. Figure 3 illustrates the difference between an instant velocity field which exhibits small structures and averaged field which shows separation and reattachment. The computation was restarted at time t = 2, to avoid influence of the averaged quantities by the initial evolution of the flow field. Presented data come from averaging at time 2 ≤ t ≤ 3.
The flow separation and reattachment are represented by the skin friction coefficient, defined by the relaxation Cf = τw / pd,1 where pd,1 is the inlet dynamic pressure given by pd,1 = ½ Uref 2 ρref . The distribution of the Cf along blade_03 is depicted in Figure 4 which shows distribution of Cf various angles of attack (α = -5°, +0° and +5°) and different numerical settings. Figure 4a and 4f show graphs for pressure side of rotated blade for α = -5° and α = +5°, respectively. It can be seen that both solvers LU-SGS and rhoSimpleFoam give the same results when using k-ω-SST turbulence model. On the contrary γ-SST model predicts lower values of Cf along about 0.05 ≤ x/c ≤ 0.4. This documents the presence of bypass laminarturbulent transition caused by combination of leading edge curvature and non-zero angle of attack. The spectral element approach gives rather low values of Cf in comparison with finite volume solvers. This can be ascribed to the absence of turbulence model and also incompressibility assumption. Figure 4b and 4e show graphs for suction side of rotated blade for α = -5° and α = +5°, respectively. One can observe in Figure 4b that prediction of flow reattachment varies considerably among different solvers and also different turbulence models of OpenFOAM. Calculation using spectral element method results in separation located further downstream the leading edge and reattachment point is not clearly determined. On the other a hand distribution of Cf and, subsequently, the reattachment point in Figure 4e is influenced mainly by γ-SST turbulence model in case of finite volume approach.   In the case of Nektar++ the plotted results show qualitatively reasonable results when compared with OpenFOAM, nevertheless, in regards to the blade_01 larger discrepancy is visible in the case of all the coefficients CD, CL and especially CM. As mentioned previously this can be ascribed to the lack of turbulence modelling, incompressible fluid assumption, but also to flow field, which developed to structures insufficiently approximated on used spatial discretisation, what was discovered from insufficient convergence of expansion coefficients in elements between trailing edge of blade_01 and bottom wall.
Data from evolutionary solver contain additional information about fluctuations of the traced quantities. It can be seen in Figure 5, that these fluctuations vary across blades. Particularly, fluctuations of CM are quite stronger on blade_01 and blade_05 in comparison with the rotated blade_03 for α = -5°. Figs. 6 and 7 illustrate the difference between different approaches used in this study, namely, time-averaged and instant Mach number filed.   Qualitative similarities are observable from comparison of Mach number distribution among steady finite volume simulations and averaged fields of incompressible result, Figure 6 and 8, namely, a higher flow rate in channels between blades 01-03 and wake structures behind blades.
Small vortical structures are observable over downstream area in Figure 7. Detail physical interpretation of these structures is problematic due to restriction to 2D and incompressibility, but sum of these instant fields surprisingly form result (time-averaged field) with relatively good coincidence as compared to those from the finite volume approach. One can also observe from comparison of CM in Figure 5 and instant field in Figure 7, that stronger vortical structures in instant flow field wake do not imply stronger fluctuations of CM.

Conclusions
Numerical simulations of aeroelastic flutter in two-dimensional flat wing cascade were carried out with aim to assess performance and accuracy different numerical solvers and turbulence models. The flutter study was realized on fixed computational meshes as an approximative approach of full fluid-structure interaction with the possibility of future extension.
The finite volume approach was realized with help of two OpenFOAM numerical solvers, namely, rhoSimpleFoam and LU-SGS. Two models of turbulence k-ω-SST and γ-SST were also compared. Numerical results of finite volume solvers show good agreement in evaluated aerodynamic coefficients. The essential aspect of chosen models of turbulence has been portrayed on skin-friction coefficient where γ-SST model differed from k-ω-SST model at extremal angles of attack.
Spectral/hp element approach, as used in time-dependent simulation of incompressible flow, has shown good stabilisation properties of the spectral vanishing viscosity technique. Processing the model and data from evolutionary solver was in strong contrast to steady state OpenFOAM simulations, but an extra information about fluctuations of traced quantities was obtained. The spectral/hp approach have shown possibilities for achievement of more accurate results, but also limits of incompressible model in this case.