Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0135 Acta Polytechnica 61(SI):135–147, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague ON THE USE OF A FLUX-SPLITTING SCHEME IN THE NUMERICAL FLUTTER ANALYSIS OF A LOW-PRESSURE TURBINE STAGE Marek Pátýa,∗, Jan Halamab a Czech Technical University in Prague, Center of Advanced Aerospace Technology, Technická 4, 16607, Prague 6, Czech Republic b Czech Technical University in Prague, Faculty of Mechanical Engineering, Department of Technical Mathematics, Technická 4, 16607, Prague 6, Czech Republic ∗ corresponding author: marek.paty@seznam.cz Abstract. The endeavour to increase the power output of steam turbines results in the design of low-pressure stages with large diameters. Such designs, featuring long and thin blades, are increasingly susceptible to unfavourable aeroelastic effects. The interaction of structure and flow may induce blade vibrations, known as flutter, which act detrimentally on the operational life of the machine. The present work employs a time-marching numerical simulation to investigate the flutter behaviour of a low-pressure transonic turbine cascade. Its blades are subject to a harmonic motion based on the results of a structural analysis and its susceptibility to flutter is evaluated via the energy method. The computations are performed with an in-house Finite Volume Method code. The flow model is based on 2D Euler equations in Arbitrary Lagrangian-Eulerian formulation with the AUSM+-up scheme for inviscid flux discretization. A higher-order spatial accuracy is achieved by using a MUSCL approach, for which both the gradient reconstruction and the slope limiting are given a careful examination-by comparing the convergence and accuracy of multiple methods. The computational model is validated by experimental data on the Fourth Standard Configuration turbine cascade. Keywords: Turbine flutter, CFD, time-marching simulation, gradient reconstruction, limiters. 1. Introduction The interaction of flow and structure leading to vi- brations has long been recognised as one of the most essential problems faced by turbine blade designers. Uncontrolled oscillations may induce high-cycle fa- tigue or even a complete failure of the machine, caus- ing extensive damage. The demand for the higher power of industrial steam turbines leads to large ex- haust areas of the last low-pressure (LP) stages with blades typically over 1 meter in length [1]. The aerody- namic requirements, especially in the near-tip region with highly supersonic flow, dictate that the thickness of the blades is reduced. As a result, the last-stage blades become increasingly susceptible to aeroelastic effects. Two types of aeroelastic phenomena are usually investigated in turbomachinery: forced response and flutter ([2, 3]). The former is a result of multiple blade row interaction, where the inlet of a downstream blade row is subjected to a circumferentially distorted flow field coming from the upstream blade row. The present aerodynamic or mechanical damping is sufficient to keep the vibrations stable and the danger of a struc- tural failure is not imminent. However, the oscillations may lead to a low or high-cycle fatigue and reduce considerably the operational life of the turbine. In contrast to a forced response, the flutter phenomenon investigated in this work represents self-excited un- stable oscillations whose magnitude grows with every cycle. The fluid-structure interaction mechanism is similar to the infamous wing flutter encountered in aeronautics, but with an added degree of complexity due to the interaction of other blades within the blade row. Aside from the danger of a fatal structural fail- ure, flutter has been reported to cause problems such as the cracking of the blade root [4]. The modelling of flutter is complicated by the num- ber of effects which need to be considered, such as the non-linear behaviour of both flow and elastic structure, change of domain and mesh in time, flow viscosity and turbulence [5]. Efforts to restrict the model to include only the most essential effects have therefore been made, such as the assumption of a potential flow and zero wing thickness in the analytical model of Theodorsen [6]. The advent of Computational Fluid Dynamics (CFD) and the rise in computational power allows for treating aeroelasticity problems with numer- ical approaches. In order to keep the computational demands reasonable, flutter was first modelled by su- perimposing linear perturbations to a steady-state non-linear solution and casting the equations to the frequency domain. However, the consideration of non- linear effects, such as large amplitudes of blade motion or strong unsteady shocks, may be crucial for the ac- curacy of the solution ([7, 8]), which was addressed by the non-linear harmonics (NLH) method introduced by Ning and He [9]. 135 https://doi.org/10.14311/AP.2021.61.0135 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Marek Pátý, Jan Halama Acta Polytechnica The inherent limitation of the frequency-domain based methods is the inclusion of only a few select harmonics in the unsteady solution. The unsteady fluid-structure interaction may be resolved fully by time-marching methods, employing often the Arbi- trary Lagrangian-Eulerian (ALE) formulation to allow a deformation of the physical domain independent of the fluid particle motion ([10–12]). The solution of the complete fluid-structure problem is complicated by the incompatibility of approaches used traditionally for the treatment of isolated structure and fluid dy- namics problems. Typically, CFD is solved with Finite Volume Methods (FVM) and Computational Struc- tural Mechanics (CSM) with Finite Element Meth- ods (FEM), differing in solution and discretization techniques [13]. Therefore, it is often simplified by imposing only a one-way coupling between the fluid and structural dynamics. It is assumed that the effect of aerodynamic forces on the change of the structural dynamic properties can be neglected in turbomachin- ery thanks to the large density ratio between fluid and structure. The blade eigenmodes and eigenfre- quencies are determined in advance in a structural solver and they are used to prescribe its motion in an unsteady fluid dynamics solver ([2, 7, 14, 15]). The motion of the neighbouring blades is not necessarily in phase and usually a range of inter-blade phase angles (IBPAs) is considered. For each IBPA, the unsteady aerodynamic load on the blade surface during one cycle is extracted and expressed by the aerodynamic damping coefficient. The susceptibility to flutter is then evaluated by comparing the present aerodynamic and mechanical damping. This approach is commonly called the energy method [3]. Although several efforts have been made to include the full two-way coupling ([16–18]), the energy method is a well-proved tool for flutter analysis [15] and it is employed in this work. The transonic or supersonic flow-field in the last stages of steam turbines is often characterised by a complex shock pattern. This requires that the nu- merical scheme, which is employed for the inviscid flux approximation in the FVM, is capable of captur- ing the shockwaves sharply, such as the flux-splitting scheme AUSM devised by Liou and Steffen [19]. In order to achieve a sufficient resolution of the flow features, many FVM strategies include a gradient reconstruction step providing a higher spatial order of accuracy. The techniques are often based on the Monotone Upstream-Centered Schemes for Conser- vation Laws (MUSCL) approach [20], comprising of two stages: firstly, the cell-wise constant values from neighbouring finite volumes are used to reconstruct the local gradient in each cell and secondly, the slope is limited to prevent spurious oscillations at the points of local discontinuities [21]. There exists a plethora of techniques for implement- ing the MUSCL-type reconstruction, and although this approach has become a routine practise in CFD, there is no general agreement regarding which one is the most suitable. In the first step, the cell-average values from the standard FVM representation are employed to reconstruct gradients in each cell. This is most commonly performed using either a Green- Gauss (GG) method, or a least squares minimiza- tion (LSQ) [22]. The GG techniques further differ by the method used to interpolate cell-face values, while the least squares employ various stencils and they may be either weighted or unweighted. In order to achieve a second order of accuracy overall, the gra- dient approximation needs to be at least first-order accurate [23]. Mavriplis [24] found that unweighted least squares severely underpredict normal gradients on some meshes, which could be remedied by the use of inverse distance weighting or by resorting to a Green-Gauss method. Sozer et al. [22] compared three GG-based approaches with a LSQ minimization and a curvilinear gradient method and demonstrated that the commonly used GG technique with either a simple or an inverse-distance weighted averaging is in- consistent and 0th order accurate for irregular meshes. They further found gradient operators with compact stencils to exhibit generally lower errors. The unsatis- factory performance of the GG method with simple averaging was confirmed by Mishriky and Walsh [25], who derived analytically that it has a 0th order of accuracy on a mesh with arbitrary spacing. They also showed that the GG method with a more suit- able face-value interpolation achieves a 1st order of accuracy, similarly to the unweighted LSQ method. Syrakos at al. [26] also concluded that a 1st order of accuracy can be achieved with a GG method via a suitable choice of face-value interpolation. They found the inverse distance weighted LSQ approach to be at least 1st accurate and, with a particular choice of the weighting exponent (-3/2), up to 2nd order accurate on some types of meshes. The second step of MUSCLE-type procedure im- poses a limiter on the reconstructed gradients in order to prevent spurious oscillations. Such limiters are for one-dimensional problems usually derived from a total-variation diminishing criterion (TVD). However, this approach is not directly extendable to multidi- mensional problems on unstructured grids, due to the difficulty of defining the TVD criterion in such cases [27]. One of the earliest attempts to apply multidi- mensional limiting on unstructured grids was devised by Barth and Jespersen [28], invoking the monotonic- ity principle to ensure that the reconstructed values within the cell do not at any point exceed the mini- mum/maximum of neighbouring centroid values. Ap- plied to a computation of flow around airfoils, their technique produced significantly more accurate results in comparison with a first-order method. However, the discontinuous nature of the limiter resulted in prob- lems with steady-state convergence. To address this problem, Venkatakrishnan [29] modified the method by using a smooth limiter function. He succeeded in achieving a better convergence behaviour, but his lim- 136 vol. 61 Special Issue/2021 Numerical flutter prediction iter used a constant which needed to be tuned for each case. Other researchers tried to exploit the multidi- mensionality of the problem and instead of multiplying the gradient by a scalar slope limiter, they used its projection onto a boundary of a maximum principle region ([27], [23]). Kim and Kim [30] developed a multidimensional limiting process for structured grids, which was extended by Park et al. to unstructured grids [31] and applied successfully to a transonic flow computation around a NACA0012 airfoil. Recently, a novel edge-based limiting procedure was devised by Delis and Nikolos ([32], [33]). Instead of defining the limited gradient as a constant value over the entire cell, the limiting procedure is performed individually for each cell-edge. The inevitable loss of accuracy accompanying the activation of the limiter can thus be minimized. Moreover, the edge-based limiting pro- cedure employs the well-known 1D limiters, allowing the user to chose the one that is most suitable for their case. The aim of this paper is to present a computational model for flutter analysis of a low-pressure turbine blade-row. Although the model employs relatively well-established techniques such as the AUSM+-up scheme of Liou [34] for inviscid flux approximation with a MUSCL-type reconstruction, the implemen- tation of each step is carefully revisited in order to achieve high fidelity. The gradients are approximated by the weighted least-squares method, for which a study analysing the suitable choice of the weighting coefficient and of the stencil extent is presented. A comparison of gradient limiting procedures is also shown, comprising of the techniques of Barth and Jaspersen [28], Venkatakrishnan [29], Park and Kim [31] and Delis and Nikolos [32]. The numerical results are validated by comparison with experimental data on the Fourth Standard Configuration (STCF4) of Bölcs and Fransson [35]. Results are also shown for a low-pressure turbine blade-row with a highly subsonic inlet and a coupled mode of vibration. 2. Computational Methods 2.1. Flow Model We denote Ωt the computational domain occupied by fluid at time t ∈ 〈0, T 〉. We aim to find the fluid density ρ, velocity u = u(x, t) and static pressure p for x ∈ Ωt, where u = [ux, uy]T has two components in the directions of Cartesian axes x, y. Further we denote the total energy et = e+ 1 2 (u2 x +u2 y), using the notation e for internal energy. We consider the domain at time t = 0 as the ref- erence domain Ωref = Ω0 and introduce the ALE mapping [5]: At : Ωref → Ωt, X → x(X, t) = At(X), (1) describing the time-dependent position x ∈ Ωt of a point from reference domain X ∈ Ωref . The ALE Figure 1. Computational domain. velocity is defined s = ∂At/∂t. We consider the set of Euler equations [36]: ∂ ∂t ∫ Ωt W dΩ + ∮ Γ(t) F (W ,n, s)dΓ = 0, (2) where n = [nx, ny]T is the surface outward unit nor- mal. The vector of characteristic variablesW and the flux vector F are defined as: W =  ρ ρu ρet  , F = (u−s)·n  ρ ρu ρet + p −  0 −pn −ps · n  . (3) The system of equations is closed by the ideal gas law: p = (κ− 1) [ ρet − ρ u2 x + u2 y 2 ] (4) We distinguish four domain locations with a differ- ent boundary condition prescribed: the blade wall Γw, the domain inlet Γin, the domain outlet Γout and the periodic boundaries Γper,l, Γper,u (Fig. 1). The blade wall is treated with a free-slip boundary condition, requiring that the normal component of flow velocity is equal to the normal component of wall movement velocity (Eq. 5). [(u− s) · n]Γw = 0 (5) Uniform stagnation conditions are prescribed at the inlet and the local characteristic variables are calculated via isentropic relations, using the Mach number extrapolated from the inside of the domain. At the outlet, the circumferentially average static pressure is prescribed. The number of blades included in the domain is chosen in function of the IBPA, such that a spatial periodicity with W |Γper,l = W |Γper,u is achieved. 2.2. Numerical Scheme The numerical solution employs a Finite Volume Method for the discretization of Euler equations in ALE formulation. The cell-face fluxes are calculated 137 Marek Pátý, Jan Halama Acta Polytechnica by the AUSM+-up scheme, introduced by Liou [34] as a modification of the original AUSM+ scheme [19]. The reader is referred to the original paper for details, here we will only expand on the alterations arising due to the flux formulation within the ALE reference frame. The underlying idea of the AUSM scheme family is to split the flux into convective and pressure parts and to treat them separately: F̃1/2 = F̃ c 1/2 + F̃ p 1/2, (6) where the 1/2 subscript indicates evaluation at the cell interface. The definition of the pressure term can be found in the paper of Darracq et al. [37], who derived the ALE formulation for an M-split AUSM+ scheme. In contrast to the original Eulerian formulation, a new term ps · n appears: F̃ p 1/2 =  0 pn ps · n  1/2 (7) The convective flux for the AUSM+-up method is defined as [34]: F̃ c 1/2 = ṁ1/2  1 u et + p/ρ  L/R (8) The mass-flow calculation needs to be updated for the ALE formulation, such that it uses a relative Mach number Mr = (u − s) · n/a instead of an absolute one: ṁ1/2 = a1/2M r 1/2 { ρL if Mr 1/2 > 0, ρR otherwise (9) Special attention needs to be paid to the integration in time for the ALE formulation. The stability and accuracy of the time-stepping method is conditioned by its compliance with the Geometric Conservation Law (GCL) which ensures an exact reproduction of a constant solution. In case of the explicit Euler scheme used for time integration in the present work, this requires that the fluxes are evaluated on a midpoint grid in between the present and the next timestep [36]. 2.3. Gradient Reconstruction and Limiting The approximation of local gradients for the MUSCL- type reconstruction is performed by a weighted least squares method. For each cell a stencil is selected, including either only its neighbours with a common face, or also those sharing only a vertex. A plane is fitted to this data, resulting into an overdetermined linear system of equations [38]:  w1(x1 − xA) w2(y1 − yA) w2(x2 − xA) w2(y2 − yA) ... ... wm(xm − xA) wm(ym − yA) ∇qA=  w1(q1 − qA) w2(q2 − qA) ... wm(qm − qA)  (10) Here the subscript A denotes the concerned cell, the subscripts 1...m its neighbours included in the stencil, x and y are the cell centroid coordinates, q is a cell-averaged quantity. Each equation is multiplied by an inverse-distance weighting factor wi = [√ (xi − xA)2 + (yi − yA)2 ]−k (11) with the exponent k ∈ {0, 1, 3 2 , 2}. Solving the system via the normal equation approach, i.e. by mul- tiplying the system of equations Ax = b with AT , can lead to a compromised accuracy on highly stretched grids [39]. The solution is therefore realized via the Gram-Schmidt process with QR decomposition. In the second step of gradient reconstruction, a limiter to avoid spurious oscillations is applied. The performance of the following limiters is analysed in the present work: • Barth & Jespersen limiter [28]: The technique limits the gradient obtained by the LSQ method in a scalar manner, defining the reconstructed quantity within the element A as: q(x, y)|A = qA + ξ∇qA ·∆r (12) First the minimum qminA and the maximum qmaxA of the cell-average values on a stencil composed of the cell A and its face-neighbours is calculated and then the limiter ξ is set such that qminneighbour ≤ q(x, y)|A ≤ qmaxneighbour (13) For a linear reconstruction the extrema of q(x, y)|A occur at the vertices of the cell, implying that the above condition needs to be applied only at these points. Sometimes a less restrictive procedure is applied, requiring that the condition is satisfied only for face midpoints. The limiter of Barth & Jespersen is also sometimes referred to as the Limited Central Difference (LCD) approach [40]. • Projected LCD [27]: The author defines a maxi- mum principle (MP) region to constrain the limited gradient ∇qA: min(qk − qA, 0) ≤ rAk · ∇qA ≤ max(qk − qA, 0) (14) for k = 1...m, where rAk is the vector from the centroid of cellA to the midpoint of the face between cells A and k. The LCD procedure simply scales the gradient such that it lies within the MP region. The procedure proposed by Hubbard, on the other hand, projects the gradient onto the boundary of the MP region. The limiter is therefore less compressive 138 vol. 61 Special Issue/2021 Numerical flutter prediction and aims at a better accuracy. It should be noted that in our implementation the gradient is obtained by the weighted least squares method, which differs from the approach employed by Hubbard. • Multidimensional Limiting Process (MLP) [31]: The MLP bears resemblance to the limiter of Barth & Jespersen in that it satisfies the maximum principle, but it differs in the extent of the stencil involved. The procedure applies to each vertex of a given cell the condition qminneighbour ≤ qvtx ≤ qmaxneighbour (15) where qvtx is the reconstructed vertex value, ob- tained by Eq. 12, and qminneighbour, q max neighbour are the minimum and maximum of the cell-average values among the cells sharing the vertex. From this con- dition the constraints for the limiting function ξ are derived. Park et al. propose two definitions of ξ: one is denoted as MLP-u1 and it represents the steepest slope in the allowable region. The limiter is non-differentiable, which can cause con- vergence problems for steady computations. The latter, MLP-u2, is a differentiable function derived from the limiter of Venkatakrishnan [29]. It also employs a constant K that needs to be tuned to balance between accuracy and convergence of the limiter. The present computations were performed with K = 1. • Face-based limiter of Delis and Nikolos. [33]: Whereas the other procedures result in a limited gradient that is constant over the whole cell, the technique of Delis et al. applies the limiter indi- vidually to each face of the cell. The procedure is performed in two steps: first, the value at the inter- section of the cell-face with a line connecting the centroids of the two adjacent cells is reconstructed, and then another reconstruction is performed for the face midpoint value. Each of these steps employs a 1D limiter, working with a virtual upwind node. The motivation behind this approach is that scalar limiters are stricter a necessary because they affect all components of the gradient uniformly. More- over, this technique employs 1D limiters, whose behaviour has been thoroughly researched in the past and whose wide array allows the user to choose the most suitable one for their case. Here we present results for the min-mod limiter. 3. Test-case Description 3.1. Fourth Standard Configuration The Fourth Standard Configuration (STCF4) belongs to a collection of turbomachinery flutter test-cases formed in the 1980’s at École Polytechnique Fédérale De Lausanne and it is presented here for the purpose of solver validation. The geometries and experimental data were made public in [35] and they are available in a digitized form at [41]. From the wide range Parameter STCF4 M8 Chord, c 72 mm 273.3 mm Blade spacing, s 56.25 mm 229.6 mm Stagger from axial, γ −56.6° −72.3° Inlet total temperature, T01 330 K 348.5 K Inlet total pressure, p01 217 100 Pa 35 637 Pa Inlet flow angle from ax., β1 −28° −73.9° Outlet isen. Mach n., M2,is 1.39 (1.42 uns.) 1.56 Vibration mode translation coupled Maximum displacement, h/c 0.00406 0.0201 Reduced frequency, λ 0.164 0.4074 Table 1. Test-case parameters. Figure 2. Computational grids for the STCF4 (top) and the M8 (bottom) cascades with insets showing refinement in the LE and TE regions. of testing conditions with publicly available data we choose the Test Case 628, characterized by transonic flow with strong compressibility effects. As the value of the inlet total temperature is not specified in the original data, we use T01 = 330 K in accordance with Waite [42]. It has been discussed that the originally reported inlet flow angle β1 = −56.6° was measured inaccurately due to the proximity of the probe to the blade leading edge [43]. A value of β1 = −28° has been demonstrated to give a better match between numerical and experimental results [42], which is also what we employ here. The mode of blade oscillation is a harmonic translation in a direction of 63° from the 139 Marek Pátý, Jan Halama Acta Polytechnica Figure 3. Test-grids for gradient reconstruction. axial. An overview of the turbine cascade parameters and testing conditions is given in the Tab. 1. The computational domain was discretised in AN- SYS ICEM with a mixed-element tri/quad grid, con- taining 13824 elements per blade passage (Fig. 2-top). The movement of the blades requires that the mesh is updated accordingly in each step. The displacement of the grid vertices forming the surface of each blade is interpolated linearly onto the domain. The result- ing displacement of each grid point is thus a sum of individual contributions from all blades. 3.2. Low-Pressure Turbine M8 The second test-case represents a cylindrical cut of a low-pressure steam turbine M8 at 65% of span. Com- pared to the Test Case 628 with an inlet Mach num- ber of 0.21, the inlet conditions are high subsonic at M1 = 0.76. Another substantial difference arises in the mode of vibration which is coupled. The elas- tic axis (EA) of the blade travels along an elliptical trajectory while the blade is subject to a pitching mo- tion about the EA. The horizontal (xEA) and vertical (yEA) displacement of the EA and the pitching angle (α) are defined as: xiEA(t) = Axcos(2πft+ φx + i · IBPA), yiEA(t) = Aycos(2πft+ φy + i · IBPA), (16) αi(t) = Aαcos(2πft+ φα + i · IBPA), i = 1...n Here the amplitudes Ax, Ay and Aα, the phase- shifts φx, φy and φα and the frequency f are obtained via an a-priori structural analysis performed in AN- SYS. The index i denotes the position of the blade in the cascade, ranging from 1 to the overall number of blades n, and the IBPA signifies the phase shift be- tween two adjacent blades. Basic test-case parameters are listed in the Tab. 1. The mixed-element tri/quad grid contains 13.107 elements per blade passage (Fig. 2-bottom) and its movement is realized similarly to the STCF4 test case. 4. Results and Discussion 4.1. Gradient Reconstruction Reconstruction of gradients from the cell-wise constant representation is performed within the MUSCL frame- work to obtain a higher-order numerical method. The present technique employs the weighted least square approximation described at the beginning of Sec. 2.3. Grid A (tri) L2-norm L∞-norm k face vertex face vertex 0 0.99 1.07 0.95 0.93 1 0.99 1.14 0.94 0.98 1.5 1.00 1.68 0.93 1.01 2 1.00 1.10 0.92 0.90 Grid B (quad) L2-norm L∞-norm k face vertex face vertex 0 1.62 1.55 0.01 0.02 1 1.64 1.72 0.01 0.15 1.5 1.64 1.76 0.00 0.23 2 1.65 1.76 0.03 0.23 Grid C (tri/quad interface) L2-norm L∞-norm k face vertex face vertex 0 1.13 1.24 0.96 0.91 1 1.11 1.36 0.88 0.98 1.5 1.11 1.76 0.86 1.00 2 1.11 1.33 0.85 0.97 Table 2. Order of accuracy of gradient approximation by least squares on three types of meshes. Results for different values of the distance-weighting exponent k and for stencils including either f ace- or vertex- adjacent cells. In order to evaluate the suitable choice for the stencil size and for the value of the weighting exponent k (Eq. 11), the accuracy of the method was tested on three types of grids. The test-grids are presented in Fig. 3: grid (A) is composed of triangular cells, grid (B) of quadrilateral cells and grid (C) comprises two sections featuring triangular and quadrilateral cells respectively. A testing function f(x, y) = tanh(jx)tanh(jy) (17) was prescribed on the domain and the error of the numerical gradient approximation was evaluated for 7 levels of refinement on each of the three grids, halving the cell-edge length each time. The order of convergence with grid refinement is quantified in Tab. 2 for the overall error (L2-norm) and for the maximum error (L∞-norm). The results are presented for different values of the weighting exponent k ∈ {0, 1, 3 2 , 2} and for two sizes of stencil, including either cells with a common face only, or also cells sharing a vertex. The performance of the least squares method on the grids (A) and (C) was similar, achieving an order of accuracy at least 0.99 in the L2-norm and over 0.85 in the L∞-norm. The quad-element grid exhibited a sub- stantially worse convergence in the L∞-norm, with the face-stencil based approximation failing to converge 140 vol. 61 Special Issue/2021 Numerical flutter prediction Figure 4. Pseudo-schlieren of a steady-state flow in the STCF4 cascade. entirely and the vertex-based approximation achiev- ing at best an order of accuracy 0.23. The L2-norm error of the quad grid, on the contrary, converged the fastest of all grid-types. The failure to converge in the L∞-norm was caused by a small region of extremely skewed cells near the lower boundary of the domain, becoming even more distorted with each refinement. Such low-quality cells would normally be avoided at the grid generation stage. The accuracy of the approximation was generally higher for the larger stencil. Only the unweighted least squares (k = 0) in some cases achieved slightly better results with the faced-based stencil. The best convergence was obtained with the weighting exponent set to k = 3 2 , in agreement with the observations of Syrakos et al. [26]. The achieved order of accuracy for k = 3 2 combined with the vertex-based stencil was over 1.68 in the L2-norm on all grids and at least 1.0 in the L∞-norm everywhere except on the quad- type grid. It can be concluded that the weighted least squares method with these settings achieves an order of accuracy of at least 1 in both types of norms, providing that extremely low-quality cells are avoided. Thus it is fully satisfactory for the construction of a second-order accurate numerical scheme and it is used for the computations presented further in this paper. 4.2. Fourth Standard Configuration Turbine Steady-state results. The topology of the steady- state flow is illustrated in the Fig. 5 by means of density gradient contours, imitating a schlieren im- age. The flow with an average inlet Mach number 0.21 is accelerated in the blade passage up to approxi- mately M = 1.6, until it encounters a right-running shockwave formed at the TE of the upper blade. The shockwave impinges on the suction side (SS) of the Figure 5. Pressure coefficient distribution (top) and convergence history (bottom) for STCF4 computa- tions with different limiters. The top plot includes experimental data [41] marked by x. lower blade, reflects from the surface at 55% chord and propagates further downstream at an angle of −30° from the machine axis. A left-running oblique shock also forms at the TE, emanating at an angle of −18° from axial direction. The blade surface pressure distribution is quantified in Fig. 5-top in terms of a non-dimensional pressure coefficient Cp. The plots show results of computations using the set of limiters described in Sec. 2.3, together with experimental data from [41] for validation. The first-order method without MUSCL-type reconstruc- tion underpredicts the flow acceleration in the first 30% of the suction side (SS), resulting in a Cp higher by about 0.05 than in the other numerical cases. More- over, it fails to capture the local pressure minimum induced by the impingement of the right running TE shockwave at 55% SS. This region also highlights the differences between the individual MUSCL-type re- constructions which are negligible elsewhere on the blade. Although the location of the peak is captured 141 Marek Pátý, Jan Halama Acta Polytechnica consistently, its magnitude varies depending on the limiter employed. The sharpest resolution is obtained by the limiter of Barth & Jespersen, whereas the u2 variation of the MLP method is shown to be the most dissipative. The other three limiters capture the pres- sure jump with similar sharpness, although with an offset in the distributions. The maximum discrepancy between all limiters amounts to ∆Cp = 0.02, while the first-order computation differs by up to ∆Cp = 0.08. The results of computations employing the MUSCL- type reconstruction agree closely with the experimen- tal data over the whole pressure side (PS) and on the section of SS upstream of the shockwave impingement. The location of the pressure minimum induced by the shock matches the experimental data well, but its mag- nitude is underpredicted by about ∆Cp = 0.11. This discrepancy is accounted for by the viscous phenom- ena in the shockwave-boundary layer interaction and to the three-dimensional effects, neither of which are modelled in the present inviscid planar computations. The impact of both viscosity and three-dimensionality on the pressure distribution in this region was demon- strated by McBean et al. [3], although it has to be noted that their computations failed to capture the shock location accurately in either case. Together with the sharpness of extrema resolution, an equally important property of a limiter is its conver- gence behaviour. The bottom part of Fig. 5 displays the evolution of density residuals in L2-norm. The convergence of the MLP-u1 and Barth & Jespersen methods stops after 5× 104 iterations and the resid- uals begin to rise again nearly to their original level. The projected-LCD and MLP-u2 limiters perform more favourably and the convergence of the latter could be further improved by raising the constant K, but this would increase the its dissipativeness, already the highest of all tested limiters. The best conver- gence was obtained by the Delis and Nikolos limiter, although the overall drop of residuals of three orders of magnitude is still not comparable with the fast convergence of the first-order computation. Neverthe- less, this limiter proves to offer the most favourable combination of extrema resolution and convergence behaviour and it is employed for the computations presented further. Unsteady results. The blades are subject to har- monic vibrations in translatory mode as described in the Sec. 3.1. Fluctuations in blade-surface pressure are analysed by performing a Fast-Fourier Transform in time at each point and the magnitude and phase of the first harmonic are displayed in the top and lower plots of Fig. 6 respectively for the IBPA 180°. The first pressure harmonic is typically the largest [3] and thus most essential for the determination of aerodynamic damping. The pressure-side distribution matches well the experiments in terms of both phase and magnitude. Satisfactory agreement is achieved also in the front 50% of the suction side, although Figure 6. Amplitude (top) and phase (bottom) of the STCF4 first harmonic pressure coefficient. with a higher discrepancy than on the PS. The fail- ure to capture viscous and three-dimensional effects, discussed already for the steady results, influences the pressure oscillations downstream of the shockwave impingement at the 55% chord. The peak at 60% chord is underpredicted by ∆|Cp(1)| = 0.6, but more importantly, the phase at the 86% chord differs by 160°. The pressure fluctuations in the aft part of the blade SS are thus captured with a nearly opposite phase, meaning that the contribution to the integral aerodynamic damping has an opposite sign. The work of aerodynamic forces during one cycle of oscillation is characterised by the aerodynamic damp- ing coefficient [44]. We use the difference between the inlet total pressure p01 and the outlet static pressure ps2 for non-dimensionalization in order to keep con- sistency with the definition of the pressure coefficient Cp: Ξ = −Waero πbh2(p01 − ps2) (18) 142 vol. 61 Special Issue/2021 Numerical flutter prediction Figure 7. STCF4 cascade global damping coefficient. Here b is the blade span and h is the maximum displacement of all points on the blade surface during the oscillation cycle. A positive value of Ξ indicates a stable flutter behaviour. The aerodynamic damping coefficient is plotted in Fig. 7 for computations with IBPAs −90°, 0°, 90° and 180°. The plot also includes experimental results for the whole IBPA range with a discretization of 18°, except for the value at 0°, unavailable in the measure- ments [41]. The sign-convention for IBPA is chosen such that a positive value marks a forward-travelling wave in the direction of blade-row rotation. In spite of the discussed shortcomings of the computational model, the prediction of aerodynamic damping shows a reasonable agreement with experiments. The trends are predicted consistently, although the numerical re- sults are shifted by 0.15-0.26 towards the region of stability. It has to be noted that this discrepancy may also arise from the method of aerodynamic damping calculation from experimental data. The imaginary part of the unsteady pressure coefficient is integrated along the blade surface and given the limited number of transducers on the blade, this procedure is bound to be rather inaccurate. 4.3. Transonic Turbine M8 Steady-state results. The steady-state blade Mach number distribution is plotted in Fig. 9 and compared to the STCF4 turbine. As the M8 test-case represents a section of low-pressure turbine in the upper portion of span, the inlet Mach number is considerably higher at 0.76 compared to M1 = 0.21 of the STCF4. The flow is further accelerated in the channel and reaches a value of 1.72 at 55% SS when it encounters the right-running TE shockwave of the neighbouring blade. The subsequent deceleration is milder and the Mach number remains nearly constant at the remainder of the SS. The PS Mach number is high subsonic with a Figure 8. Pseudo-schlieren of a steady-state flow in the M8 cascade. Figure 9. Steady-state blade Mach number distribu- tions of the M8 and STCF4 turbine cascades. supersonic portion in the aft 11% chord induced by the flow acceleration in the inter-blade channel. The flow topology is illustrated by contours of the density gradient in Fig. 8. The shockwave pattern is similar to the STCF4, but the compression induced by the pair of left- and right-running TE shockwaves is weaker. The density gradients marking the wake region are also notably less pronounced. Unsteady results. The global aerodynamic damp- ing coefficient is plotted in Fig. 10 to characterise the aeroelastic stability of the turbine cascade. The computations were performed for the whole range of IBPA with a stepping of 30°. The damping coefficient is negative in all investigated cases, indicating an un- stable behaviour. Its values seem to lie on a nearly sinusoidal curve with the exception of the IBPA −30°. 143 Marek Pátý, Jan Halama Acta Polytechnica Figure 10. M8 cascade global damping coefficient. Flat-plate cascade model predictions of acoustic reso- nance inter-blade phase angles highlighted. The aerodynamic excitation is significantly less intense in this case with the damping coefficient standing out by 0.25 from the distribution. Such a sharp local drop or peak in aerodynamic damping is commonly referred to as a "flutter bite" and it is often a manifestation of acoustic resonance [42]. This condition describes a limiting state when acoustic waves induced by the blade movement are on the verge of propagating throughout the domain. Lane and Friedman [45] studied propagation of acoustic waves in a 2D flat plate cascade and derived an analytical relation for resonant IBPAs. We use an equivalent formula given by Donini [14]: σr = ( Mη ± √ 1−M2 ξ ) Mλs (1−M2)c , (19) whereMξ andMη are the axial and tangential Mach numbers respectively, λ is the reduced frequency, s the blade-to-blade spacing and c the chord. Using ei- ther the inlet or the outlet Mach number and reduced frequency, we obtain two pairs of resonant IBPAs. In case of the M8 turbine cascade, the resonant an- gles are σ1 r,in = −17.1°, σ2 r,in = 119.6° for the inlet and σ1 r,out = −52.3°, σ2 r,out = −12.8° for the outlet conditions. Given the simplifications of Lane’s and Friedman’s analytical model and the IBPA resolution of 30° of investigated cases, it cannot be expected that the resonant angles would match exactly the results of the computation. However, the location of the IBPA −30° in the neighbourhood of the first inlet and of both outlet resonant angles supports the conlusion that the local peak in damping is caused by acoustic resonance. 5. Conclusions This paper presents the results of a numerical flut- ter investigation for a low-pressure turbine blade row. The in-house solver is based on Euler equations for planar flow, discretized by the Finite Volume Method. The AUSM+-up scheme is used for flux approxima- tion and a MUSCL-type reconstruction is employed to construct a higher-order scheme. A weighted least- squares gradient approximation was tested for several types of grids and it was shown that the best accuracy is achieved with a vertex-based stencil and a weight- ing exponent − 3 2 . Such approximation achieved at least 1st order accuracy on all grids and it is thus considered sufficient for the construction of a 2nd or- der MUSCL reconstruction. Several gradient-limiting techniques were tested on a steady-state computation of the Fourth Standard Configuration turbine cascade. Although the method of Barth & Jespersen managed to capture shockwaves with the sharpest resolution, it suffered from a rather poor convergence behaviour. The face-based limiting technique of Delis and Nikolos offered the most favourable combination of improved convergence properties and accurate extrema resolu- tion. Both steady-state and unsteady computations of the STCF4 turbine showed a good overall agreement with experimental data. There were observed dis- crepancies in terms of the local pressure magnitude during suction-side shockwave impingement and in the unsteady pressure phase-shift further downstream along the suction side. Both these discrepancies were ascribed to the influence of viscous and three- dimensional effects which are not included in the nu- merical model. The global aerodynamic damping pre- diction for different inter-blade phase angles matched the measurements in terms of the overall trends, but the absolute values showed an offset. Apart from the computational model simplifications, this disagree- ment can also be ascribed to the limited accuracy of damping coefficient evaluation from experimental data. The aeroelastic behaviour was further analysed for a low-pressure turbine M8, operating at high subsonic inlet conditions and excited in a coupled mode of vibration. The global damping coefficient followed a sinusoidal curve over a range of inter-blade phase angles with a strong deviation observed for the IBPA −30°. An analytical model for the propagation of acoustic waves in a 2D flat plate cascade showed that the probable cause of this damping peak is acoustic resonance. 6. Acknowledgements Authors gratefully acknowledge support from the ESIF, EU Operational Programme Re- search, Development and Education, and from the Center of Advanced Aerospace Technology (CZ.02.1.01/0.0/0.0/16_019/0000826), Faculty of Mechanical Engineering, Czech Technical University in Prague, and from the Grant Agency of CTU in Prague, grant No. SGS19/154/OHK2/3T/12. 144 vol. 61 Special Issue/2021 Numerical flutter prediction List of symbols A amplitude ALE Arbitrary Lagrangien-Eulerian b blade span [m] c chord [m] Cp pressure coefficient, (ps − ps2)/(p01 − ps2) Cp(n) n-th harmonic of unsteady pressure coefficient, ps(n)/[h/c(p01 − ps2)] CFD Computational Fluid Dynamics CSM Computational Structural Mechanics h maximum blade displacement [m] e energy [J kg−1] EA elastic axis f frequency [Hz] F flux vector FEM Finite Element Methods FVM Finite Volume Methods GCL Geometric Conservation Law GG Green-Gauss IBPA inter-blade phase angle k weighting exponent LCD Limited Central Difference LE leading edge LP low-pressure LSQ least squares M Mach number ṁ mass-flow [kgm−2 s−1] MLP Multidimensional Limiting Process MP Maximum Principle MUSCL Monotone Upstream-Centered Schemes for Conservation Laws n surface normal [m] NLH non-linear harmonics p pressure [Pa] s ALE velocity [ms−1] s blade-to-blade spacing [m] t time [s] TE trailing edge TV D total-variation diminishing u velocity [ms−1] W vector of characteristic variables w weighting factor Waero work of aerodynamic forces in one cycle [J ] x, y coordinates [m] β flow angle from axial [deg] γ stagger angle from axial [deg] Γ domain boundary κ Poisson constant λ reduced frequency, 2πfc/uref ξ limiter function Ξ aerodynamic damping coefficient ρ density [kgm−3] σr resonant inter-blade phase angle [deg] φ phase-shift [rad] Ω computational domain 0 stagnation 1 inlet 2 outlet is isentropic t total ref reference w wall l lower u upper per periodic r relative c convective p pressure 1/2 cell interface L left R right vtx vertex s static η tangential ξ axial References [1] P. Petrie-Repar, V. Makhnov, N. Shabrov, et al. Advanced Flutter Analysis of a Long Shrouded Steam Turbine Blade. vol. Volume 7B: Structures and Dynamics of Turbo Expo: Power for Land, Sea, and Air. 2014. V07BT35A022, doi:10.1115/GT2014-26874. [2] W. Höhn. Numerical Investigation of Blade Flutter at or Near Stall in Axial Turbomachines. Ph.D. thesis, Royal Institute of Technology, Stockholm, Sweden, 2012. [3] I. McBean, K. Hourigan, M. Thompson, F. Liu. Prediction of Flutter of Turbine Blades in a Transonic Annular Cascade. Journal of Fluids Engineering 127(6):1053–1058, 2005. doi:10.1115/1.2060731. [4] P.-A. Masserey, I. McBean, H. Lorini. Analysis and improvement of vibrational behaviour on the nd37 a last stage blade. VGB Powertech Journal 92:42–48, 2012. [5] P. Sváček, M. Feistauer, J. Horáček. Numerical simulation of flow induced airfoil vibrations with large amplitudes. Journal of Fluids and Structures 23(3):391 – 411, 2007. doi:10.1016/j.jfluidstructs.2006.10.005. [6] B. Perry. Report no. 496: General theory of aerodynamic instability and the mechanism of flutter. Technical Report, NASA. Langley Research Center, Hampton, Virginia, United States, 2015. [7] H. Doi, J. J. Alonso. Fluid/Structure Coupled Aeroelastic Computations for Transonic Flows in Turbomachinery. vol. 4: Turbo Expo 2002, Parts A and B of Turbo Expo: Power for Land, Sea, and Air, pp. 787–794. 2002. doi:10.1115/GT2002-30313. [8] R. Kamakoti, W. Shyy. Fluid-structure interaction for aeroelastic applications. Progress in Aerospace Sciences 40:535–558, 2005. doi:10.1016/j.paerosci.2005.01.001. [9] W. Ning, L. He. Computation of Unsteady Flows Around Oscillating Blades Using Linear and Nonlinear Harmonic Euler Methods. Journal of Turbomachinery 120(3):508–514, 1998. doi:10.1115/1.2841747. 145 http://dx.doi.org/10.1115/GT2014-26874 http://dx.doi.org/10.1115/1.2060731 http://dx.doi.org/10.1016/j.jfluidstructs.2006.10.005 http://dx.doi.org/10.1115/GT2002-30313 http://dx.doi.org/10.1016/j.paerosci.2005.01.001 http://dx.doi.org/10.1115/1.2841747 Marek Pátý, Jan Halama Acta Polytechnica [10] J. Donea, S. Giuliani, J. Halleux. An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions. Computer Methods in Applied Mechanics and Engineering 33(1):689 – 723, 1982. doi:10.1016/0045-7825(82)90128-1. [11] D. Boffi, L. Gastaldi. Stability and geometric conservation laws for ale formulations. Computer Methods in Applied Mechanics and Engineering 193(42):4717 – 4739, 2004. doi:10.1016/j.cma.2004.02.020. [12] P. H. Saksono, W. G. Dettmer, D. Perić. An adaptive remeshing strategy for flows with moving boundaries and fluid-structure interaction. International Journal for Numerical Methods in Engineering 71(9):1009–1050, 2007. doi:10.1002/nme.1971. [13] A. Slone, C. Bailey, M. Cross. Dynamic solid mechanics using finite volume methods. Applied Mathematical Modelling 27(2):69 – 87, 2003. doi:10.1016/S0307-904X(02)00060-4. [14] N. Donini. Aeroelasticity of turbomachines linearized flutter analysis. Ph.D. thesis, Politecnico di Milano, Facoltá di Ingegneria Industriale, Milano, Italy, 2012. [15] M. May, Y. Mauffrey, F. Sicot. Numerical flutter analysis of turbomachinery bladings based on time-linearized, time-spectral and time-accurate simulations. In IFASD 2011 - 15th International Forum on Aeroelasticity and Structural Dynamics. 2011. [16] A. Slone, K. Pericleous, C. Bailey, M. Cross. Dynamic fluid-structure interaction using finite volume unstructured mesh procedures. Computers & Structures 80(5):371 – 390, 2002. doi:10.1016/S0045-7949(01)00177-8. [17] R. Rzadkowski, V. Gnesin. 3-d inviscid self-excited vibrations of a blade row in the last stage turbine. Journal of Fluids and Structures 23(6):858 – 873, 2007. doi:10.1016/j.jfluidstructs.2006.12.003. [18] V. Gnesin, L. Kolodyazhnaya, R. Rzadkowski. A numerical modelling of stator-rotor interaction in a turbine stage with oscillating blades. Journal of Fluids and Structures 19(8):1141 – 1153, 2004. doi:10.1016/j.jfluidstructs.2004.07.001. [19] M.-S. Liou, C. J. Steffen. A new flux splitting scheme. Journal of Computational Physics 107(1):23 – 39, 1993. doi:10.1006/jcph.1993.1122. [20] B. van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of Computational Physics 32(1):101 – 136, 1979. doi:10.1016/0021-9991(79)90145-1. [21] J. Blazek. Computational Fluid Dynamics: Principles and Applications. Elsevier, 2001. [22] E. Sozer, C. Brehm, C. C. Kiris. Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered CFD solvers. In 52nd Aerospace Sciences Meeting. 2014. doi:10.2514/6.2014-1440. [23] M. Berger, M. Aftosmis, S. Muman. Analysis of slope limiters on irregular grids. In 43rd AIAA Aerospace Sciences Meeting and Exhibit. 2012. doi:10.2514/6.2005-490. [24] D. Mavriplis. Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes. In 16th AIAA Computational Fluid Dynamics Conference. 2003. doi:10.2514/6.2003-3986. [25] F. Mishriky, P. Walsh. Towards understanding the influence of gradient reconstruction methods on unstructured flowsimulations. Transactions of the Canadian Society for Mechanical Engineering 41(2):169–179, 2017. [26] A. Syrakos, S. Varchanis, Y. Dimakopoulos, et al. A critical analysis of some popular methods for the discretisation of the gradient operator in finite volume methods. Physics of Fluids 29(12):127103, 2017. doi:10.1063/1.4997682. [27] M. Hubbard. Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids. Journal of Computational Physics 155(1):54 – 74, 1999. doi:10.1006/jcph.1999.6329. [28] T. Barth, D. Jespersen. The design and application of upwind schemes on unstructured meshes. In 27th Aerospace Sciences Meeting. 1989. doi:10.2514/6.1989-366. [29] V. Venkatakrishnan. On the accuracy of limiters and convergence to steady state solutions. In 31st Aerospace Sciences Meeting. 1993. doi:10.2514/6.1993-880. [30] K. H. Kim, C. Kim. Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part ii: Multi-dimensional limiting process. Journal of Computational Physics 208(2):570 – 615, 2005. doi:10.1016/j.jcp.2005.02.022. [31] J. S. Park, S.-H. Yoon, C. Kim. Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids. Journal of Computational Physics 229(3):788 – 812, 2010. doi:10.1016/j.jcp.2009.10.011. [32] A. I. Delis, I. K. Nikolos. A novel multidimensional solution reconstruction and edge-based limiting procedure for unstructured cell-centered finite volumes with application to shallow water dynamics. International Journal for Numerical Methods in Fluids 71(5):584–633, 2013. doi:10.1002/fld.3674. [33] A. Delis, I. K. Nikolos. On a solution reconstruction and limiting procedure for unstructured finite volumes. American Institute of Mathematical Sciences (AIMS) Journal 8:491–499, 2014. [34] M.-S. Liou. A sequel to AUSM, part II: AUSM+-up for all speeds. Journal of Computational Physics 214(1):137 – 170, 2006. doi:10.1016/j.jcp.2005.09.020. [35] A. Bölcs, T. Fransson. Aeroelasticity in Turbomachines: Comparison of Theoretical and Experimental Cascade Results + Appendix A5. EPFL Lausanne, 1986. [36] R. W. Smith. AUSM(ALE): A geometrically conservative arbitrary Lagrangian-Eulerian flux splitting scheme. Journal of Computational Physics 150(1):268 – 286, 1999. doi:10.1006/jcph.1998.6180. [37] D. Darracq, S. Champagneux, A. Corjon. Computation of unsteady turbulent airfoil flows with an aeroelastic AUSM+ implicit solver. In 16th AIAA Applied Aerodynamics Conference. 1998. doi:10.2514/6.1998-2411. 146 http://dx.doi.org/10.1016/0045-7825(82)90128-1 http://dx.doi.org/10.1016/j.cma.2004.02.020 http://dx.doi.org/10.1002/nme.1971 http://dx.doi.org/10.1016/S0307-904X(02)00060-4 http://dx.doi.org/10.1016/S0045-7949(01)00177-8 http://dx.doi.org/10.1016/j.jfluidstructs.2006.12.003 http://dx.doi.org/10.1016/j.jfluidstructs.2004.07.001 http://dx.doi.org/10.1006/jcph.1993.1122 http://dx.doi.org/10.1016/0021-9991(79)90145-1 http://dx.doi.org/10.2514/6.2014-1440 http://dx.doi.org/10.2514/6.2005-490 http://dx.doi.org/10.2514/6.2003-3986 http://dx.doi.org/10.1063/1.4997682 http://dx.doi.org/10.1006/jcph.1999.6329 http://dx.doi.org/10.2514/6.1989-366 http://dx.doi.org/10.2514/6.1993-880 http://dx.doi.org/10.1016/j.jcp.2005.02.022 http://dx.doi.org/10.1016/j.jcp.2009.10.011 http://dx.doi.org/10.1002/fld.3674 http://dx.doi.org/10.1016/j.jcp.2005.09.020 http://dx.doi.org/10.1006/jcph.1998.6180 http://dx.doi.org/10.2514/6.1998-2411 vol. 61 Special Issue/2021 Numerical flutter prediction [38] A. C. Haselbacher. A grid-transparent numerical method for compressible viscous flows on mixed unstructured grids, 1999. [39] W. Anderson, D. L. Bonhaus. An implicit upwind algorithm for computing turbulent flows on unstructured grids. Computers & Fluids 23(1):1 – 21, 1994. doi:10.1016/0045-7930(94)90023-X. [40] P. Batten, C. Lambert, D. M. Causon. Positively conservative high-resolution convection schemes for unstructured elements. International Journal for Numerical Methods in Engineering 39(11):1821–1838, 1996. doi:10.1002/(SICI)1097- 0207(19960615)39:11<1821::AID-NME929>3.0.CO;2- E. [41] Fourth standard configuration experimental data. https://www.rpmturbo.com/testcases/STCF/STCF4_ update/stcf4_update.htm. Accessed: 2020-05-27. [42] J. J. Waite. Physical Insights, Steady Aerodynamic Effects, and a Design Tool for Low-Pressure Turbine Flutter. Ph.D. thesis, Duke University, Durham, North Carolina, 2016. [43] T. H. Fransson, J. M. Verdon. Updated report on standard configurations for the determination of unsteady flow through vibrating axial-flow turbomachine- cascades. Proceedings of the 6th International Conference on Aeroelasticity in Turbomachines 1991. [44] D. Schlüß, C. Frey. Time domain flutter simulations of a steam turbine stage using spectral 2D non-reflecting boundary conditions. Proceedings of the 15th International Symposium on Unsteady Aerodynamics, Aeroacoustics & Aeroelasticity of Turbomachines, ISUAAAT15. 2018. [45] F. Lane, M. Friedman. Technical note 4136: Theoretical investigation of subsonic oscillatory blade-row aerodynamics. Technical Report, New York University, New York, United States, 1958. 147 http://dx.doi.org/10.1016/0045-7930(94)90023-X http://dx.doi.org/10.1002/(SICI)1097-0207(19960615)39:11<1821::AID-NME929>3.0.CO;2-E http://dx.doi.org/10.1002/(SICI)1097-0207(19960615)39:11<1821::AID-NME929>3.0.CO;2-E http://dx.doi.org/10.1002/(SICI)1097-0207(19960615)39:11<1821::AID-NME929>3.0.CO;2-E https://www.rpmturbo.com/testcases/STCF/STCF4_update/stcf4_update.htm https://www.rpmturbo.com/testcases/STCF/STCF4_update/stcf4_update.htm Acta Polytechnica 61(SI):135–147, 2021 1 Introduction 2 Computational Methods 2.1 Flow Model 2.2 Numerical Scheme 2.3 Gradient Reconstruction and Limiting 3 Test-case Description 3.1 Fourth Standard Configuration 3.2 Low-Pressure Turbine M8 4 Results and Discussion 4.1 Gradient Reconstruction 4.2 Fourth Standard Configuration Turbine 4.3 Transonic Turbine M8 5 Conclusions 6 Acknowledgements List of symbols References