Analysis of classical flutter in steam turbine blades using reduced order aeroelastic model

In the present paper classical flutter phenomena in LP steam turbine rotor is studied . A reduced order aeroelastic model (ROAM) with acceptable accuracy and fast execution is developed for this purpose. The aerodynamics damping (AD) is estimated using Traveling wave mode (TWM) method. Flow field is modeled using Panel method. For the structural part ROM non-linear beam element method (BEM) based FEM structural solver is used. Partitioned based ( loose) coupling approach is adopted to perform aeroelastic (flutter) cosimulation. Both 2D cascade flow and 3D cascade are modeled.The estimated stability parameters are compared with experimental data. Moreover, present ROAM shows significant reduction in computational time.


Introduction
In the quest to generate more power in the nuclear power plants giant steam turbines are being installed.The large low pressure steam turbine blades are more prone to mechanical vibration due to longer size and less stiffness, which leads to high-cycle fatigue failure and in worse case blade loss and damaging the whole system.The aeroelasticity of low-pressure turbine blade has been the subject of much attention in the past and recent years [1,2].Aerodynamic damping is considered as an crucial parameter to quantify the classical flutter or aeroelastic stability of the system.Especially for turbine blisks and centrifugal compressors where there is no friction between blades and disk and small material damping, the damping from the surrounding fluid is the main damping contributor.Therefore, large numbers of design iterations are required to optimize the blade for aeroelastic stability in short period of time.There are well developed high fidelity CFD-CSD available and are being used for product development purpose as an aeroelastic simulation tool.These methods provide very accurate representation of aeroelastic system if modeled accurately, however, they are computationally very costly.Therefore,fast ROAM model are best suited for the preliminary stage for LP steam turbine design and development.
The present research work is devoted to development of medium fidelity such ROAM tools for the aeroelastic simulation of steam turbine rotor.To model the flow field computationally less expensive potential flow based engineering model are used to estimate the aeroelastic parameters, e.g.aerodynamics damping for the blade cascade.The panel method (PM) is one of such methods and first proposed by Hess and Smith [3] to model the lifting and non-lifting potential flow around slender bodies.These methods are good compromise of MATEC Web of Conferences 211, 15001 (2018) https://doi.org/10.1051/matecconf/201821115001VETOMAC XIV speed and accuracy and can be used for complex geometry unless the flow fulfill the criteria of potential flow and without separation.These methods are widely adopted for aeroelastic modeling of wind turbines, helicopter rotors, and aircraft aeroelasticity, a short review is presented by [4,5].The more details about theoretical and numerical implementation can be found in Katz [6].However, instead of great potential of the PM to be used in LP turbine aeroelasticity very few researchers have used it in this field.A improved version of PM is used by [7,8] for cascade design and flow modeling, in his work only pressure distribution is estimated, whereas [9] used the similar technique to estimate the aeroelastic stability parameters, but in this work only camber surface of the blade is modeled not the actual geometry.For structural part Euler-Bernoulli beam theory based finite element model (FEM) is used to model the long steam turbine blades, the beam models are wildly used in wind turbine, helicopter rotor and propeller rotor structural modeling [10].The LP turbine blades are modeled as rotating cantilever flexible beams with 2 degree of freedom.Euler-Bernoulli beam model is good compromise of 3D complete FEM model, in terms of accuracy and speed.For the fluid structure coupling partitioned based loose coupling strategy is used where flow solver (PM) is loosely coupled with Structural solver (beam FEM) and both solver exchange date at predeceased time step.
To estimate the aerodynamic damping most common method include traveling wave mode (TWM) method and less known aerodynamic influence coefficient (AIC) method [11,12].The TWM method is used by many research to estimate the aerodynamics damping for vibrating cascade.In the present work TWM method is used for the AD calculation, a detail explanation of both the methods can be found in Fransson [11].Both 2D and 3D cascade aeroelastic modeling is carried out and results are compared against experimental results.

Unsteady panel method for flow modeling
Panel methods are the family of boundary element methods, and can be applied to simulate the aerodynamic flow over lifting and non lifting bodies.The panel method gives the solution of Laplace's potential flow equation Eq. 1, therefore, the flow field assume to be inviscid and irrotational.
where Φ(x, y, z) is the scalar velocity potential and it is function of (x, y, z) .There are many analytical solution which satisfy the Laplace's equation, these are called singularities solution e.g.vortex, source, doublet or mix of either of them.The detail formulation of unsteady PM can be found in Katz [6].

Beam element FEM model of bladed disk
The structural part of the blades are modeled with Euler-Bernoulli beam FEM .The beam based structural models are adopted by many researchers for aeroelastic modeling of turbomachinary rotors [13,14], however, according to authors knowledge it is not used for LP turbine blades modeling.Therefore, beam based FEM is used for structural modeling here.The 1D/2D beam elements are faster than 3D FEM models.An they can simulate the main structural dynamic parameters part in aeroelastic simulation.In the present work rotor blades are modeled as rotating cantilever beam with 1D beam element with 2 DOF namely , bending and torsion.More details about Euler-Bernoulli beam formulation can be found in any stander FEM text book.A shamanistic picture of rotating cantilever beam , beam element Fig. 1. and Structural model of Bladed disk is given in the Fig. 2. https://doi.org/10.1051/matecconf/201821115001VETOMAC XIV

Fluid structure interaction (FSI) coupling and co-simulation
As mentioned in the introduction section the FSI is based on partitioned approach, where both the FEM and flow solvers are loosely coupled and execute interdependently.Both solver are coupled using MATLAB/SIMULINK S-function environment.A schematic diagram of FSI is presented in the Fig. 3.More details about aerodynamic and structural grid mapping and co-simulation is presented in the [15].

Aeroelastic simulation of 2D cascade using ROAM
The unsteady PM described in sec. 2 is adopted to estimate the AD in the 2D cascade flow.The TWM method is used for the analysis.The cascade numbering, TWM and INFC method schematically presented in the top left of Fig. 3 more details about the method is given [11].
For the 2D test case NACA65 series airfoil is selected and flow condition are similar to Carta [16] which have constant inlet velocity of 61 m/set (200 ft/sec), mean camber line incidence angles (α 0 = 6 • ), pitching amplitudes ( ᾱ = 0.5 • ), reduced frequency k = 0.122 based on semi chord, and 8 IBPA (ϕ) = 0, ±45 • , ±90 • ,±135 • , ±180 • .More details about https://doi.org/10.1051/matecconf/201821115001VETOMAC XIV the experimental setup and data acquisition is presented in [16].This kind of superimposed velocity potential can be obtained for any airfoil in the cascade.The oscillatory motion of the cascade blade can be given in complex form by Eq. 2 where n is blade number Fig. 3 and the ω is the frequency of oscillation in rad/sec and ω = 2.k.U ∞ chord .Hence the effective angle of attack (α e f f ) at any time instance is given by Eq. 3 Comments: The simulated aerodynamic damping (Ξ) for α 0 = 6 • and ᾱ = 0.5 • at reduce frequency k = 0.122 and wind speed U ∞ = 61.2,Re = 4.0e + 06, Mach = 0.18 is compared with the experimental result as presented by Carta [16].In the simulation IBPA is −180 • to +180 • with increment of 10 • is used.In the Fig. 5 simulated aerodynamic damping shows good agreement with experimental results, and simulated results capture the stable and unstable zone with great accuracy and the stability limit (Ξ = 0) are well predicted, however, in the unstable zone the PM based method over estimate the aerodynamics damping magnitude over > +60 • IBPA.The main cause of this overshoot is not clear yet, but it might be caused by disagreement in phase lead or acoustic resonance at those IBPA, however it is to be investigated in more details.

Aeroelastic simulation for 3D annular cascade using ROAM
After successful implementation of the method for 2D cascade, the ROAM based solver is adopted to estimate the AD in 3D blade row.For the test case experimental setup presented by [17] is selected.The simulated AD for different IBPA is compared against experimental results for nominal inflow angle (α)= −26 • and Mach=0.28 and k=0.1 [18].The test setup and corresponding PM model (50 chorwise x15 span panels) is given in the Fig. 6 and Fig. 7.
First the steady flow span-wise pressure distribution over the reference blade is estimated at nominal flow condition and presented in the Fig. 8.The AD for the same flow condition for torsional mode oscillation at different IBPA is compared in Fig. 9 The simulated AD using 3D PM flow model at different IBPA in torsional vibration in TWM agrees well with experimental results in Fig. 9 and PM based solve accurately catches the S − curve running from unstable to stable zone.However, the ROAM model overestimates the magnitude of  the AD at higher IBPA in the stable region.This can be caused by acoustic resonance or flow separation from the neighboring oscillating blades.Since there is no flow separation or acoustic resonance model in the present ROAM, therefore, it over estimates the magnitude.
In the calculation of AD only the effects of ∓1 and ∓2 on the reference blade is considered, this is to be more consistence with experimental setup.It is worth mentioning that the stability boundaries are based on the -sign and + sign of AD value, the -sign of AD value leads to unstable system and probable cause of flutter,whereas + sign of AD value keeps the system stable and no chance of flutter.The stability curve vs IBPA in Fig. 5 and Fig. 9, follows asymptotic behavior in stable and unstable zone and can be noticed by s − curve in the figures.

Conclusion
A medium fidelity ROAM has been successful develop and implemented to simulate the flow in LP power turbine.The aerodynamic damping for low pressure turbine both 2D and 3D cascade is estimated.The present method is good compromise of speed and accuracy.The simulated results for aerodynamic damping has good agreement with experimental results.Therefore, the present method can be applied to estimate the aeroelastic stability parameters for low frequency and low amplitude oscillation cases.Furthermore, in the ROAM reduced order structural model based on beam theory coupled with flow solver (PM) using partitioned based FSI approach.This type of model will significantly reduce MATEC Web of Conferences 211, 15001 (2018) https://doi.org/10.1051/matecconf/201821115001VETOMAC XIV the product development time especially for turbomachinery field, where virtual prototyping is a key aspect during R& D. It will also provide engineers and researchers active in the aeroelastic field more flexibility and ease to design more complex and aeroelastically stable system in short time span.The proposed model will also facilitate the both time and frequency domain aeroelastic simulation.Moreover, the present ROAM prove to be computationally more than 10 times cheaper compared to traditional CFD-CSD models.