Efficient prediction of classical flutter stability of turbomachinery blade using the boundary element type numerical method

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Efficient prediction of classical flutter stability of turbomachinery blade using the boundary element type numerical method Chandra Shekhar Prasad∗, Luděk Pešek Institute of Thermomechanics, Prague, Czech Republic

a r t i c l e

i n f o

Keywords: Turbomachinery Aeroelasticity Panel method (PM) Aerodynamic damping (AD) Traveling wave mode (TWM) Classical flutter Boundary element method (BEM)

a b s t r a c t In this paper classical flutter phenomena in power turbine rotor is studied. A medium fidelity 2D flow solver based on boundary element method is developed for this purpose. The classical flutter parameters in turbomachinery cascades such as aerodynamic damping at different inter-blade phase angle are estimated using flutter stability analysis flow solver. The flow solver is developed using 2D unsteady potential flow based panel method. The unsteady aerodynamic loading on the vibrating cascade is estimated using a traveling-wave mode oscillation. The simulated aerodynamic damping and pressure coefficient using boundary element method are compared against both experimental data and the computational fluid dynamics model’s results at different flow conditions. The boundary element based method results demonstrate good agreement with experimental data. The boundary element based flow solver shows significant reduction in computational time compared to computational fluid dynamics model.

1. Introduction The worldwide increasing need for electrical energy and the same time government’s policies around the world to reduce the dependency on fossil fuels to generate electricity force the modern nuclear power plants to consider installing more efficient and even larger steam turbines. An example of such gigantic steam turbine is Arabelle steam turbine by GE power, which has low pressure (LP) turbine blade length of 75 inches to capture the massive volume of steam produced by nuclear power plants. The Arabelle turbines are now used in one-third of the world’s nuclear power stations. These large LP stage steam turbine blades are more prone to mechanical vibration due to longer size and thinner dimension, thus less stiffer. Therefore, aeroelastic stability of low pressure (LP) turbine has been the subject of much attention in the past and recent years [1,2] to ensure uninterrupted and safe power production in nuclear power plants. Flutter is one of such common aeroelastic instability phenomena occurs in turbomachines e.g. power turbines, jet engines. There are different types of flutter depending upon the flow condition inside turbomachines e.g. dynamic stall flutter, choke flutter, shock-induced flutter, supersonic non-stalled flutter, classical (or low amplitude vibration) flutter, Fransson [3]. However, in this research work only the classical flutter is analyzed and studied which is well suited to underlying assumptions in potential flow. The classical flutter is usually defined as a type of flutter when the flow is fully attached to the blade surface in all time and it involves directly the

∗

phase lag between the blade oscillation and the induced unsteady timedependent aerodynamic forces on the blade cascade. Therefore, aerodynamic damping is considered as an important 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, thus the damping from surrounding fluid is the main damping contributor. The aerodynamic damping is nothing but work done by the fluid on the vibrating structure per oscillation cycle. Although in the most cases the problem caused by unsteady aerodynamic forces are of non-negligible amplitude, it is much easier to only study the low amplitude classical flutter phenomena, this is because of the tendency towards flutter must be ignored and if no self excitation will appear at small amplitudes then it will not grow further. This straight forward or linearized way of looking at the problem of aeroelastic instability is only justified as long as only flutter stability limits are studied. However, in some cases, nonlinear effects may have to be taken into account to accurately estimate the flutter stability limits during design phase [3]. Therefore, engineers and designers in this field need to invest a fair amount of time and resources to insure the aeroelastic stability issue well within the optimal operating condition, thus large numbers of design iterations are requited within short period of time to optimize the blade for aeroelastic stability. There are well developed computational fluid dynamics (CFD) and computational structural dynamics (CSD) based aeroelastic models available for this purpose e.g. ANSYS/Fluent or

Corresponding author. E-mail addresses: [email protected] (C.S. Prasad), [email protected] (L. Pešek).

https://doi.org/10.1016/j.enganabound.2020.01.013 Received 30 January 2019; Received in revised form 21 January 2020; Accepted 23 January 2020 0955-7997/© 2020 Elsevier Ltd. All rights reserved.

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

CFX, NUMECA, Start CCM+ etc, however, theses CDF-CSD based models are computationally very expensive, time consuming and required highly skilled manpower and supercomputers to solve. In CFD-CSD based aeroelastic models flow computation is responsible for majority of computational time. Therefore, CFD based methods can not be the best choice for industrial use where product development time is vital. Therefore, in this research paper, application of computationally less expensive potential flow based boundary element method (BEM) flow model is studied to model the flow in low pressure (LP) turbine blade cascade and to estimate the aeroelastic stability (classical flutter) parameters e.g. aerodynamic damping for the 2D blade cascade. The panel method (PM) is one of such methods and first proposed by Hess and Smith [4] to model the lifting and non-lifting potential flow around slender bodies. These methods are good compromise of speed and accuracy and can be used for complex geometry until the flow fulfills the criteria of potential flow and free from flow separation which is exactly the case of classical flutter flow condition. However, modified PM can be used for the separated flow case [5]. These methods are widely adopted for aeroelastic modeling of wind turbines, helicopter rotors, and aircraft aeroelasticity problem [6,7]. The more details about theoretical and numerical implementation of PM can be found in Katz [8]. Instead of having great potential of the PM based BEM to be used in LP turbine classical flutter study, a very less attention have been paid to it in the past. Most of the researchers or engineers have adopted Euler or Navier Stokes (NS) based CFD model for unsteady flow modeling in turbomachinery application and very few researchers have adopted PM based BEM for unsteady flow modeling blade cascade to study aeoelastic stability more specifically classical flutter. An improved surface singularity method similar to PM is used by McFariand [9], Lei and Liang [10] for cascade blade design and flow modeling for compressible inviscid flow in 2D cascade. Both researchers in their work have simulated only pressure distribution and the Mach number over airfoil surface in steady condition. Furthermore, Barbarossa et al. [11] used the similar technique to estimate the aeroelastic stability parameters in LP turbine, however, in his work only camber surface of the blade is modeled using lifting line method not the actual geometry is considered. This makes the Barbarossa et al. [11] method little restrictive to be used for complex/thick airfoil shapes, moreover, only pressure difference over airfoil’s reference surface is estimated not the pressure distribution over actual surface of the airfoil. Therefore, in the present work 2D unsteady surface panel method with Dirichlet boundary condition (BC) is adopted to discretized the actual blade cascade airfoil surface with distribution of constant strength source and doublet singularity elements for flow field construction. The key reason for it, that the lifting line theory based model used by Barbarossa et al. [11] or similar vortex lattice method (VLM) based model ignores the thickness of the airfoil and the boundary condition are applied on the mean reference surface and the pressure difference (Δp) distribution is estimated over mean surface, instead of actual pressure distribution over lower and upper surface, whereas PM has no such restrictions, this makes the PM based BEM flow solver more versatile without any further restriction compared to VLM or lifting line model. Furthermore, to estimate the aerodynamic damping most common method are traveling wave mode (TWM) method and less known aerodynamic influence coefficient (AIC) method [12,13]. The TWM method is used by many researchers to estimate the aerodynamic damping in vibrating cascade. In the TWM all the blades vibrate with identical mode shape, frequency and amplitude, but maintains a constant inter-blade phase angle (IBPA) between the oscillating neighboring blades. In contrast to TWM in the AIC method only a single blade (reference blade) in the cascade vibrates and influence of this on its non-oscillating blade are superimposed to estimate the aerodynamics damping, a detail explanation of both the methods can be found in Fransson [12]. The estimated aerodynamic damping using either of these two methods are theoretically identical [12,14]. In this research work TWM method is adopted because it is more representative to actual vibration pattern in the power turbine rotor. In reality all the blades in the turbine rotor

Fig. 1. Cross-section of blade, wake, tangent, normal and flow domain.

oscillate simultaneously with respective eigenmode and a fixed frequency and multiple nodal diameters (ND). Due to appearance of these nodal diameters of the disk the neighboring blades in the blade row experience a phase shift, termed as inter-blade phase angle (IBPA). Depending on the IBPA a traveling wave appears to be traveling with (forward traveling) or against (backward traveling) the direction of rotation [14]. Unlike TWM oscillation it is not very common to have only single oscillating blade in a rotating bladed disk. Moreover, the TWM oscillation provides most accurate data for the controlled flutter experiments [15]. Therefore, TWM technique is used for estimation of aerodynamic damping in the oscillating cascade. The estimated results for NACA65 series airfoil profile 2D cascade are compared against the experimental results presented by Carta [16] and same has been compared with CFD based simulated results. For the CFD simulation ANSYS/Fluent 17.2 is used with double precision accuracy setting. 2. Unsteady panel method for flow modeling As it is discussed earlier PMs are the family of boundary element methods used in simulation of aerodynamic flow over lifting and nonlifting bodies. The panel method gives the solution of Laplace’s potential flow equation Eq. (1), therefore, the flow field assumed to be incompressible, inviscid and irrotational. Laplace’s equation is given as ∇2 Φ(𝑥, 𝑦, 𝑧) =

𝜕2 Φ 𝜕2 Φ 𝜕2 Φ + + =0 𝜕2 𝑥 𝜕2 𝑦 𝜕2 𝑧

(1)

where Φ(x, y, z) is the scalar velocity potential and it is function of (x, y, z) or (x, y) for 2D. In the unsteady PM to construct the flow field and to find the solution of the flow field variables, singularity element panels are placed over the upper and lower surface of the blade and on the time varying wakes. There are number of well known singularity elements available for this purpose e.g. vortex singularity elements, source singularity elements or doublet singularity elements. To construct the flow filed, these elements can be used either individually or in a combination of two which depends upon the user’s choice and the problem. In this research work constant strength source (𝜎) and doublet (𝜇) singularity elements are used. Therefore, each surface panel includes combination of a source and a doublet distribution while the wake panels feature only doublet panels, since it dose not carry any aerodynamic load. More details about choice of singularity elements can be found in Katz [8]. The boundary conditions are imposed on control points lying on the mid-point of each panel, Figs. 1 and 2 represents a 2D crosssection of 3D blade with source and doublet elements distribution. The blade’s airfoil surface is denoted by Sb , that of the wake by Sw and the far-field boundary by S∞ . The velocity potential for the onset flow (outer flow region) is denoted by Φ∞ and inner velocity potential is Φi (inside the bounded region of airfoil). The flow region of interest is the one lying outside the airfoil and the one inside the airfoil is a fictitious flow 329

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 2. The velocity potential near a solid boundary Sb .

which is single equation with three unknowns, 𝜇b , 𝜎 b and 𝜇 w . The doublet and source strengths are defined as

field. The Fig. 1 also depicts the free stream Q∞ and an infinitesimal area dS on the airfoil’s surface whose influence on point P(x, y, z) is sought. The vector n is a unit vector normal to the surface at dS. The source (𝜎) and doublet (𝜇) singularity element are the know analytical solutions to Laplace’s equation, therefore, they automatically satisfy the Eq. (1). Subsequently, the panels have been arranged so as to represent the real blade geometry. However, the strength of the panels will determine the flow field and this has not been calculated yet. In order to perform this calculation three boundary conditions (BC) must be applied. These BC are given as

− 𝜇 = Φ − Φ𝑖

𝜕Φ 𝜕Φ𝑖 − = ∇(Φ − Φ𝑖 ) ⋅ 𝐧 (8) 𝜕𝑛 𝜕𝑛 The Eq. (6) itself is an integral and cannot be solved for a general blade geometry. The solution is to discretize the geometry, as discussed earlier, and replace the integrals by sums, and apply Eq. (6) to the control point of each panel, such that is can be expressed as set of linear algebraic equation −𝜎 =

1. Impermeability: The blade’s surface is a solid boundary, therefore, all flow velocity components normal to it must vanish, i.e. ∇Φ ⋅ 𝐧 = 0

𝐀𝝁𝑏 + 𝐁𝝈𝒃 + 𝐂𝝁𝑤 = 𝟎

(2)

lim ∇(Φ − Φ∞ ) = 0

(3)

𝐀(𝑁𝑏 𝑥𝑁𝑏 )

⎛ 𝑎1,1 ⎜𝑎 = ⎜ 2,1 ⎜ ⋮ ⎜𝑎 ⎝ 𝑁𝑏 ,1

𝑎1,2 𝑎2,2 ⋮ 𝑎𝑁𝑏 ,2

… … ⋱ …

𝑎1,𝑁𝑏 ⎞ 𝑎2,𝑁𝑏 ⎟⎟ ⋮ ⎟ 𝑎𝑁𝑏 ,𝑁𝑏 ⎟⎠

𝐁(𝑁𝑏 𝑥𝑁𝑏 )

⎛ 𝑏1,1 ⎜𝑏 = ⎜ 2,1 ⎜ ⋮ ⎜𝑏 ⎝ 𝑁𝑏 ,1

𝑏1,2 𝑏2,2 ⋮ 𝑏𝑁𝑏 ,2

… … ⋱ …

𝑏1,𝑁𝑏 ⎞ 𝑏2,𝑁𝑏 ⎟⎟ ⋮ ⎟ 𝑏𝑁𝑏 ,𝑁𝑏 ⎟⎠

𝐂(𝑁𝑏 𝑥𝑁𝑤 )

⎛ 𝑐1,1 ⎜𝑐 = ⎜ 2,1 ⎜ ⋮ ⎜𝑐 ⎝ 𝑁𝑏 ,1

𝑐1,2 𝑐2,2 ⋮ 𝑐𝑁𝑏 ,2

… … ⋱ …

𝑐1,𝑁𝑤 ⎞ 𝑐2,𝑁𝑤 ⎟⎟ ⋮ ⎟ 𝑐𝑁𝑏 ,𝑁𝑤 ⎟⎠

where r is the distance from point P(x, y, z) to the body and Φ∞ = 𝐐 ∞ ⋅ 𝐱

(4)

for 𝐱 = 𝑥𝐢 + 𝑦𝐣 + 𝑧𝐤, is the potential induced by the free stream. This boundary condition is satisfied automatically because source and doublet panels are singularities and their effects become negligible at great distances. This statement also applies to vortex rings. The difference Φ′ = Φ − Φ∞ is known as the perturbation potential [8]. 3. Kutta condition: While modeling flow over slender body e.g. airfoil shape, the flow must leave the trailing edge smoothly which is one of the criteria of potential flow. Kutta condition is applied at the trailing edge of the airfoil to ensure this criteria. This condition is enforced by appropriate choice of the strengths of the first row of wake panels whose strength is equal to the difference upper and lower trailing edge surface panels.

𝝁𝑏 = [𝜇𝑏1 … 𝜇𝑏𝑁 ]𝑇 is the vector of the unknown doublet strengths of 𝑏

the body surface panels, 𝝈 𝑏 = [𝜎1 … 𝜎𝑁𝑏 ]𝑇 is the vector of the unknown source strengths of the body surface panels, 𝝁𝑤 = [𝜇𝑤1 … 𝜇𝑤𝑁 ]𝑇 is the 𝑤 vector of the unknown doublet strengths of the wake panels, Nb is the number of body panels, Nw is the number of wake panels. Each elements of influence coefficient matrices A i.e. ai,j represents the induced potential on 𝑖th collocation point by 𝑗 th constant strength body doublet panel, similarly element of matrices B i.e. bi,j represents the induced potential on 𝑖th collocation point by 𝑗 th constant strength body source panel and element of matrices C i.e. ci,j represents the induced potential on 𝑖th collocation point by 𝑗 th wake doublet panel. For the isolated airfoil/blade the matrices A and B are evaluated once at beginning of time step, but for the oscillating cascade they need to be evaluated at each time step to include the effect of vibrating neighboring blades. More details about solution of Eq. (9) can be found in [8]. The solution of Eq. (9) for each collocation points give the values of three unknowns, 𝜇b , 𝜎 b and 𝜇 w . Once the perturbation potential or unknown 𝝁b is known, the local tangential velocity components at each collocation point can be calculated by differentiating the perturbation potential along the tangential direction 𝜇 ( 𝑗 , 𝑡 ) − 𝜇𝑏 ( 𝑗 + 1 , 𝑡 ) 𝑞𝑖𝑛𝑑 (𝑗 , 𝑡) = 𝑏 (10) Δ𝑙(𝑗 )

Furthermore, Katz and Plotkin [8] have shown that, after applying one of Green’s identities to the area of flow between Sb , Sw and S∞ , the potential induced at point P can be written as Φ=

( ) ( ) 1 1 1 1 d𝑠 − d𝑠 𝜇𝑏 𝐧 ⋅ ∇ 𝜎𝑏 ∫ ∫ 4𝜋 𝑆𝑏 𝑟 4𝜋 𝑆𝑏 𝑟 ( ) 1 1 + 𝜇 𝐧⋅∇ d 𝑠 + Φ∞ 4𝜋 ∫𝑆𝑤 𝑤 𝑟

(9)

where A, B, C are influence coefficient matrices and can be given as

2. Far field: The disturbance created by the singularities on the body and wake must disappear at infinity, i.e. |𝐫|→∞

(7)

(5)

where 𝜇 b is the strength of the body doublet distribution, 𝜎 b the strength of the body source distribution and 𝜇 w the strength of the wake doublet distribution, while 𝑟 = |𝐫| is distance between the point P and the body panel. In the present formulation internal Dirichlet boundary condition is used, therefore, by equating the inner potential Φi = onset potential Φ∞ , the Eq. (5) will reduce to ( ) ( ) ( ) 1 1 1 1 1 1 𝜇 𝐧⋅∇ 𝜎 𝜇 𝐧⋅∇ d𝑠 − d𝑠 + d𝑠 = 0 4𝜋 ∫𝑆𝑏 𝑏 𝑟 4𝜋 ∫𝑆𝑏 𝑏 𝑟 4𝜋 ∫𝑆𝑤 𝑤 𝑟

where j is the index of 𝑗 th collocation point at which the velocity is being calculated and Δl(j) is length of 𝑗 th panel. The velocity component

(6) 330

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 3. Schematic diagram of pitching airfoil with unsteady PM based BEM.

normal to the surface will be zero, due to zero normal flow BC. The unsteady pressure coefficient is calculated using instantaneous Bernoulli’s equation as proposed by Katz and Plotkin [8] and can be given as 𝑐̃𝑝 (𝑥, 𝑦, 𝑡) = 1 −

𝑄2𝑡𝑜𝑡𝑎𝑙 𝑄2∞

−

2 𝜕𝜇 𝑄2∞ 𝜕𝑡

Wake vorticies tracking and propagation: In real flow field these wake vortices significantly effects the unsteady aerodynamic forces on the upstream body. Therefore, it is crucial to accurately model it and include their effect to accurately estimate the unsteady aerodynamic forces using PM. In theory wake length can be measured from starting end point of body moving through undisturbed fluid. Also once theses wake vortices are shed from the body’s trailing edge their circulation strength remains unchanged as mentioned earlier and since they do not carry any force they tend to propagate with stream local velocity. The local velocity is a result of the kinematic motion of the body and the velocity components induced by the wake and body and is usually measured in the inertial frame of reference X, Z, (Fig. 3) at each panel’s end point’s (for 2D line panel). To achieve the vortex wake roll-up, at each time step the induced velocity (u, w)i at each vortex wake point is calculated and then the vortex elements are moved by Eq. (13)

𝑤ℎ𝑒𝑟𝑒 𝑄𝑡𝑜𝑡𝑎𝑙 = 𝑉𝑙𝑜𝑐𝑎𝑙−𝑘𝑖𝑛𝑒𝑚𝑎𝑡𝑖𝑐 + 𝑞𝑖𝑛𝑑 (11)

It is worth mentioning that getting to the solution for the cp is not so straightforward. In unsteady PM which is based on underlying principle of circulatory based unsteady aerodynamics, it, crucial to correctly deal with shed vortices at each time step. In the present flow modeling the strength of nascent vortices is achieved with the help of unsteady Kutta condition implementation which is accordance with Kelvin’s theorem. Kutta condition and enforcement of Kelvin’s theorem: To determine the strength of nascent wake vortices at each time step and to avoid the abrupt pressure jump and keep the flow tangency at airfoil’s trailing edge, 2D unsteady Kutta condition can be applied along the trailing edge of the airfoil. Kutta condition is accordance with Kelvin’s theorem, which says the rate of change of bounded circulation around the continuous body (here it is airfoil + wake) is zero, thus, any change in total circulation around the airfoil due its motion is compensated by the nascent wake circulation shed from the trailing edge of the airfoil. Therefore, to comply with the Kelvin’s theorem, unsteady Kutta condition can be evoked. In unsteady Kutta condition the latest wake panel strength is set such that the over all circulation around airfoil and wake become zero, and it also ensures smooth wake separation from trailing edge without any abrupt jump in pressure or velocity. The strength of nascent doublet wake panel is calculated by making it equal to the difference of airfoil’s trailing edge upper and lower panel’s doublet strength, and this also makes the total vorticity at trailing edge zero. In the present formulation the nascent trailing edge wake’s doublet strength is the difference of two trailing edge panel induced potential and can be given as 𝜇(𝑡)𝑤 = 𝜇(𝑡)𝑢𝑡𝑒 − 𝜇(𝑡)𝑙𝑡𝑒

(Δ𝑥, Δ𝑧)𝑖 = (𝑢, 𝑤)𝑖 .Δ𝑡

(13)

In the present formulation of unsteady PM modeling the effect of trailing edge wakes on the body surface at each time step is included according to Eq. (9), where matrices C represents the wake induced potential at each body collocation point and vector 𝝁w the vectors of know wake doublet strength at each time step. The reasonable results can be obtained by tracking and including the wake far from the trailing edge. In the present unsteady PM model this distance is taken to be 20 chord length behind the trailing edge in downstream and it is obtained by wake convergence study. A simulated unsteady wake shape model of 2D 11 blade cascade using PM is presented in the Fig. 12. A more detailed discussion on numerical implementation and determining the aerodynamic loads can be found in [8]. In the present work 2D constant strength source (𝜎) and doublet (𝜇) singularity element are used to construct the flow field around both isolated oscillating airfoil Fig. 3 and blade cascade Fig. 8. The analytical expression of induced potential Φind at any arbitrary point (x, z) by constant strength 2D source element in panel coordinate system is given by Eq. 14

(12) 𝑢 𝜇𝑡𝑒

where 𝜇 w wake doublet strength, doublet strength of upper trailing 𝑙 doublet strength of lower trailing edge panel. Thus, the edge panel, 𝜇𝑡𝑒 vorticity lost by the airfoil at any time step become equivalent to the circulation of this new wake vortex. Furthermore, soon after the nascent wake vortex is formed it is cast off the airfoil and left behind spinning in the air where the airfoil left it, and it’s strength do not change with time [8]. Again a new nascent vortex with new vortex strength is released in the next time step and so on. This process continues till the end of time step. Therefore, at each time step the new nascent wake panels are released in to the wake with varying strength, whereas, once these nascent wake panel are separated from the trailing edge, they remain spinning in the air and propagate with local velocity, but their strength remain unchanged in time. Readers are encouraged to read more in details about vortex strength estimation using Kelvin’s theorem and Kutta condition in Katz [8].

𝜎 Φ𝑖𝑛𝑑 (𝑥, 𝑧) = 4𝜋

{ (𝑥 − 𝑥1 ) ln [(𝑥 − 𝑥1 )2 + 𝑧2 ] − (𝑥 − 𝑥2 ) ln [(𝑥 − 𝑥2 )2 + 𝑧2 ]

( + 2𝑧 tan−1

𝑧 𝑧 − tan−1 𝑥 − 𝑥2 𝑥 − 𝑥1

)} (14)

where (x1 , z1 ) and (x2 , z2 ) is the two ends of source element and 𝜎 is source strength. Similarly the induced potential by doublet panel is given by Eq. (15) ] [ 𝜇 𝑧 𝑧 Φ𝑖𝑛𝑑 (𝑥, 𝑧) = − tan−1 − tan−1 (15) 2𝜋 𝑥 − 𝑥2 𝑥 − 𝑥1 where (x1 , z1 ) and (x2 , z2 ) is the two ends of doublet element and 𝜇 is doublet strength of the panel. The entire numerical code is implemented in MATLAB-2017b environment. This unsteady PM is first applied to 331

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 4. Grid convergence test for the pitching airfoil PM model.

Fig. 5. No stall: Cyclic Lift coefficient for pitching NACA0012 airfoil, at 𝑘 = 0.1, 𝛼𝑒𝑓 𝑓 = 3◦ + 10◦ sin(𝜔𝑡).

model the flow around isolated pitching airfoil and the results are compared against experimental results. The step wise produce is presented in following section.

edge at each time step. The latest shed wake panel have unknown doublet strength and can be estimated in each time step by using unsteady Kutta condition, Fig. 3. Since the wakes do not carry any aerodynamic loads, therefore, only doublet panels will be used to discretize the wake panels. Furthermore, the doublet strength of wake remains unchanged once they are shed and therefore, allowed to propagate downstream freely with local free stream velocity. Dirichlet BC condition is applied at each collocation point and unsteady Kutta condition is imposed at the trialing edge of the pitching airfoil to estimated the unknown induced potential over the airfoil surface at each time step. The time varying velocity, pressure and other aerodynamic loads can be obtained by differentiating the induced potential. The airfoil’s pitching axis is mid chord with pitching frequency 𝜔 Fig. 3.

2.1. Pitching airfoil test case In this section flow modeling around pitching airfoil using PM based BEM is described. In order to model the flow around oscillating pitching airfoil using unsteady panel method, at the first step entire airfoil boundary surface is discretized into 2D panels, and constant strength source and doublet singularity elements are distributed over it. The tangent (𝜏) and normal (𝜂) are defined at the collocation points of each panel element Fig. 3, time varying wake panels are shed from trailing 332

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 6. Deep stall: Cyclic Lift coefficient for pitching NACA0012 airfoil, at 𝑘 = 0.1, 𝛼𝑒𝑓 𝑓 = 10◦ + 10◦ sin(𝜔𝑡).

2.1.1. Grid convergence test of PM model Prior to use the 2D PM model to estimate the unsteady aerodynamic parameters for the pitching airfoil test case, a grid convergence study is performed to determine the optimal numbers of panel discretizetion. The grid convergence is carried out in unsteady mode without pitch oscillation and lift coefficient (Cl ) at 6∘ angle of attack for NACA0012 airfoil at wind speed 103 m/s (Mach = 0.3) is estimated for increasing number of panels. The grid convergence results are presented in the Fig. 4 In the Fig. 4 blue line denotes the Cl value w.r.t to numbers surface panels and the red line denotes the CPU time (min) taken w.r.t number of panels. It clear from the Fig. 4 that the Cl start to converge at 60 numbers of panel, hence, 64 numbers of panel is used as optimal number of panels in the simulation for the pitching airfoil test case.

cept at the starting of upstroke and at the end of down stroke. However, the large discrepancy between results are not unexpected at this points. In the deep stall flow test case the mean flow angle is very high 10∘ , therefore, flow start of separate earlier stage of upstroke and became fully separated making the airfoil stalled when 𝛼 eff approaches to 15∘ . The flow remains separated till almost the end of down stroke cycle consequently PM estimation of lift for separated flow dominated regime is very poor and not acceptable. Therefore, it is clear form above comparisons that the present method (classical PM) can be applicable with good accuracy for low amplitude oscillation where no significant flow separation occurs on the airfoil surface, but not suitable for large amplitude or incident angle oscillation case where flow field is dominated by separated flow. Since the present study is carried out on classical flutter (or low-pitch magnitude vibrations), which is a low amplitude and low reduced frequency flutter phenomena by definition and involves no flow separation, thus the use of PM based BEM for flutter stability study is self justified and can be used successfully.

2.1.2. Results For the test case a pitching NACA0012 airfoil is modeled and simulated results are compared with experimental data presented by McCroskey et al. [17]. There are two different flow conditions are simulated 1) no stall flow case 2) deep stall flow case. At first no stall flow case is simulated and estimated cyclic lift coefficient is compared with experimental data and presented in the Fig. 5 for reduced frequency k = 0.1, 𝛼𝑒𝑓 𝑓 = 3◦ + 10◦ sin(𝜔𝑡), Pitch axis = c/4, Aspect ratio (AR) = ∞, Mach = 0.3, Chord = 0.62 m. In the Fig. 5 the simulated cyclic lift coefficient shows good agreement with experimental results for both upstroke and down-stroke cycle. However, there is slight over estimation of the Cl at starting of the down-stroke cycle, but, the discrepancy in the Cl gradually converge with the experimental data by approximately mid span of down-stroke cycle. For the no stall test case at the peak of upstroke and the very beginning of down stroke cycle, the 𝛼 eff is 13∘ which not the stalling angle for the NACA0012 airfoil, however, due to start of pitch down motion of the airfoil there is limited flow separation with strong vortices forms even 𝛼 eff start to decrease Fig. 5 (top airfoil), this abrupt flow separation for short period of time called “stall onset” phenomena [17]. The stall onset causes decrease in lift in the experimental result during starting of down stroke cycle. Since classical PM which is adopted here can not simulated the flow separation, therefore, overestimates the lift coefficient value. This limited separation soon disappears as shown in the Fig. 5 (second airfoil from top) and PM estimation become more accurate. However, the results can be further improved by adopting PM with flow separation model as presented by Prasad and Dimitriadis [5]. In the second step deep stall flow case is simulated and estimated cyclic lift coefficient is compared with experimental data and presented in the Fig. 6 for reduced frequency 𝑘 = 0.1, 𝛼𝑒𝑓 𝑓 = 10◦ + 10◦ sin(𝜔𝑡), Pitch axis = c/4, Aspect ratio (AR) = ∞, Mach = 0.3, Chord = 0.62 m. In the Fig. 6 simulated Cl have very poor match with experimental result through out the entire pitch cycle, ex-

3. Classical flutter phenomena in turbomachinery cascade In this section brief discussion about classical flutter is presented. Classical flutter as discussed earlier in the introduction section is the type of aeroelastic instability which generally occurs in the turbomachinery cascade due to phase difference between oscillating blades row and the induced unsteady aerodynamic force. Depending upon the time lag which can be appear due to different physical reason between oscillating blades and the unsteady harmonic aerodynamic forces acting on the blades. The oscillating blade either supply energy into the flow field or absorbs energy from the flow field or remains neutral during the specific vibration cycle. A schematic illustration is presented in the Fig. 7 for better understanding of the phenomena, where h(𝜔t), 𝑝̄(𝜔𝑡) and Ξ(𝜔t) are harmonic oscillatory motion (pitch or plunge) of the blade, harmonic pressure and work done by harmonic aerodynamic forces respectively, and 𝜔 is angular frequency and t is the time. Therefore, if oscillating blades gives energy to the flow then any arising oscillation will be damped and in this case it is called damped vibration in aeroelastic point of view and it is presented in sub-figure (b) of Fig. 7, whereas in the opposite case where oscillating blades absorbs energy from the flow then the vibration amplitude of oscillation will gradually increase and can cause flutter sub-figure (a) of Fig. 7. Furthermore, if there is no exchange of energy between flow and structure the vibration is called neutral during specific vibration cycle. The behavior of structure at certain period of time, thus depends upon the unsteady characteristics of fluid flow field and the structure overall vibration cycle. In the sub-figure (c) and sub-figure(d) in the Fig. 7, the stall flutter and the classical flutter 333

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 7. Schematic illustration of flow interaction with oscillating structure and types of flutter. (a) flutter: aeroelastically unstable (b) Damped vibration: aeroelastically stable (c) Stall flutter (d) Classical (or low-pitch magnitude vibrations) flutter.

Fig. 8. Schematic diagram of flow field construction for cascade row using unsteady PM.

334

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 9. UTRC Oscillating Cascade Wind Tunnel (OCWT); Carta [16].

vibrating in the TWM with IBPA (𝜑), the unsteady pressure coefficient 0,𝜑 can be expressed in complex pressure coefficient form 𝑐̂𝑝,𝐴,𝑡𝑤𝑚 (𝑥, 𝑦, 𝑡) and it represents the combined influence of reference blade on itself and the neighbouring blades and the wakes. A more details about complex form of pressure coefficient is presented in the following section.

phenomena is presented respectively, where M1 and M2 are inlet and outlet Mach number of the flow. According to Fung [18] [page.161] the classical flutter is the “Oscillatory instability in a low amplitude vibration, in which neither flow separation nor the strong shock waves are involved.” Therefore, the classical flutter is well suited to underlying assumptions in potential flow. Unlike classical flutter case in stall flutter flow separation and/or shock waves can appear. Now from the above discussion about classical flutter it is clear that the potential flow based method can be adopted to estimate the classical flutter parameters and therefore, 2D unsteady PM is used in the this research work to model the unsteady harmonic aerodynamic flow field around the oscillating blade cascade.

4. Aerodynamic damping calculation In this section estimation of aerodynamic damping coefficient for classical flutter case is described in details. The balance between the structure and aerodynamic forces in the system is represented by the aeroelastic Eq. (16) [𝑀]{𝑋̈ } + [𝐺]{𝑋̇ } + [𝐾]{𝑋} = {𝐹𝑎𝑑 (𝑡)}

3.1. Panel method for 2D cascade flow in turbomachinery

= {𝐹𝑑𝑖𝑠𝑡𝑢𝑟𝑏𝑎𝑛𝑐𝑒 (𝑡) + 𝐹𝑑𝑎𝑚𝑝𝑖𝑛𝑔 (𝑡)}

In this section modeling of unsteady flow field around 2D cascade using PM based BEM is discussed. In Fig. 8 the PM discretization of a blade cascade is presented, where m represents the number of surface panels, and the indexing start from lower trailing edge marching till upper trailing edge e.g. 1 to 2 m. The tangent (𝜏) and normal (𝜂) are defined at each collocation points and BCs are imposed at the same. All the blades in the cascade are discretized in the same fashion like isolated airfoil case. Now in order to determine the distribution of unknown doublet strength (or perturbation velocity potential) on the reference blade (middle blade 0 in Fig. 8) Eq. (6) is evaluated on each collocation point and method of superimposition is applied. The influence coefficient matrices A0 , B0 , C0 are evaluated for the reference blade (0 subscript for reference blade) for each time step, these influence coefficient matrices contain the induced potential by reference blade’s doublet and source singularity elements on itself and also by the neighbouring blade’s doublet and source elements as well as from wake vortices at each time step. Once the influence coefficient matrices are evaluated, the unknown time varying strength of doublet panel can be estimated in the similar fashion as it is done for the isolated pitching airfoil case. Once the perturbation potential or unknown 𝝁b is known, the unsteady pressure coefficient is calculated using instantaneous Bernoulli’s equation same manner as pitching airfoil case. More details about evaluation influence coefficient for the multiple airfoil case presented in the [8]. For the blade cascade

(16)

where [M], [G], and [K] are the modal mass, modal damping and modal stiffness matrices respectively,{X} represents the modal coordinate vector which comprise the bending and torsion part, and finally {Fad (t)} is the modal unsteady aerodynamic force vector which contains two terms: Fdisturbance (t) that includes the aerodynamic disturbances upstream and downstream the blade and Fdamping (t) which represents the aerodynamic damping resulting from the interaction between the blade itself and the flow. For flutter analysis, only the aerodynamic forces due to the vibration of the blade are considered, resulting in 𝐹𝑑𝑖𝑠𝑡𝑢𝑟𝑏𝑎𝑛𝑐𝑒 (𝑡) = 0 and simplifying the above equation to: Eq. (17) [𝑀]{𝑋̈ } + [𝐺]{𝑋̇ } + [𝐾]{𝑋} = {𝐹𝑑𝑎𝑚𝑝𝑖𝑛𝑔 (𝑡)}

(17)

The solution of the Eq. (17) allows us to determination the aerodynamic damping and thus, the aeroelastic stability of the system. There are different numerical methods which can be used to approximate the aerodynamic characteristics of the system. One of them is “Time Linearized Method” and it assumes that the unsteady perturbations in the flow are small compared with the mean flow, then the unsteady flow can be approximated by small harmonic perturbations around a mean value. Therefore, in this work PM based BEM is used to estimate the harmonic aerodynamic forces on the vibrating blade cascade. Verdon [13] has shown that for the small perturbation the unsteady pressure due to harmonic motion can be represented as harmonic oscillation with respect 335

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 10. Cascade geometry used for experiment; Carta [16].

to the blade around a time averaged steady value and can be given as 𝑝(𝑥, 𝑦, 𝑡) = 𝑝̄(𝑥, 𝑦) + 𝑝̃(𝑥, 𝑦, 𝑡) = 𝑝̄(𝑥, 𝑦) + 𝑝̂.𝑒

𝑖(𝜔𝑡+𝜑𝑝→ℎ )

ds is equal to the surface panel length i.e. dl. Generally the complex unsteady pressure coefficient and the force are normalized w.r.t to the amplitude (A) of motion (e.g. for torsional flutter the amplitude is the amplitude of oscillation) and the dynamic pressure head which generally upstream of the flow field and can be given by Eq. (23)

(18)

where 𝑝̄(𝑥, 𝑦) is the steady mean pressure at a arbitrary location (x,y), 𝑝̃(𝑥, 𝑦, 𝑡) is the unsteady perturbation pressure and 𝑝̂ is the amplitude of complex pressure due to harmonic oscillation and it is estimated here using unsteady PM. Further the unsteady pressure coefficient is given by Eq. (19) 𝑐̃𝑝 (𝑥, 𝑦, 𝑡) =

𝑝̃(𝑥, 𝑦, 𝑡) 𝑃𝑑𝑦𝑛,𝑟𝑒𝑓

𝑤ℎ𝑒𝑟𝑒

𝑃𝑑𝑦𝑛,𝑟𝑒𝑓 =

1 2 𝜌𝑄 2 ∞

𝑐𝑝,𝐴 ̂ =

𝑤ℎ𝑒𝑟𝑒

= (𝑐̃𝑝 (𝑥, 𝑦, 𝑡)𝑅 + 𝑖.𝑐̃𝑝 (𝑥, 𝑦, 𝑡)𝐼𝑚 ).𝑒𝑖(𝜔𝑡)

𝑐̂𝑝 (𝑥, 𝑦, 𝑡) =

𝑝̂ 𝑃𝑑𝑦𝑛,𝑟𝑒𝑓 (20)

where 𝑐̂𝑝 (𝑥, 𝑦, 𝑡) is the amplitude of the unsteady complex pressure at any time instance at point (x, y, t), 𝑐̃𝑝 (𝑥, 𝑦, 𝑡)𝑅 and 𝑐̃𝑝 (𝑥, 𝑦)𝐼𝑚 are real and imaginary components of the complex unsteady pressure coefficient respectively. Both real and imaginary components can further be expressed in the amplitude and IBPA terms as

𝑐𝑓 =

̂ 𝑓⃗ = (21)

∮

𝑑 𝑝̂.𝑒𝑖(𝜔𝑡+𝜑𝑝→ℎ ) .𝑛. 𝑑𝑠

𝐴.𝑃𝑑𝑦𝑛,𝑟𝑒𝑓

(24)

∮

⃗̂ . 𝑑 𝑠 𝑑𝑓

(25)

⃗̂ per unit and the three components of infinitesimal unsteady force 𝑑𝑓 surface span is given according the axis definition in the Fig. 8

The harmonic nature of the perturbation unsteady pressure consequently results in harmonic unsteady aerodynamic forces and it can be obtained by integrating harmonic unsteady pressure around the airfoil, therefore, can be expressed by Eq. (22) ̃ ̂ 𝐹⃗ = 𝐹⃗ .𝑒𝑖(𝜔𝑡+𝜑𝑝→ℎ ) =

̃ 𝐹⃗

and normalized unsteady force per unit span can be given by Eq. (25)

𝑐̃𝑝 (𝑥, 𝑦, 𝑡)𝑅 = 𝑐̂𝑝 (𝑥, 𝑦, 𝑡) cos(𝜑𝑝→ℎ ) 𝑐̃𝑝 (𝑥, 𝑦, 𝑡)𝐼𝑚 = 𝑐̂𝑝 (𝑥, 𝑦, 𝑡) sin(𝜑𝑝→ℎ )

(23)

In the present research work potential flow based PM method is used as a linearized flow solver assuming the small unsteady perturbation. Once the unsteady pressure coefficient is determined, then the pressure coefficient is integrated over the surface to get the unsteady harmonic forces and the work done for each cycle is obtained by multiplying the total force by respective transnational or rotational displacement. The unsteady aerodynamic force obtained from pressure interrogation ̃ is represented by 𝐹⃗ . The unsteady force coefficient cf is obtained by ̃ normalizing 𝐹⃗ and can be given by Eq. (24)

(19)

the Eq. (19) can be rewritten either in complex exponential form or in component form as 𝑐̃𝑝 (𝑥, 𝑦, 𝑡) = 𝑐̂𝑝 (𝑥, 𝑦, 𝑡).𝑒𝑖(𝜔𝑡+𝜑𝑝→ℎ )

𝑝̂ 𝐴.𝑃𝑑𝑦𝑛,𝑟𝑒𝑓

(22)

where 𝑑 𝑝̂.𝑒𝑖(𝜔𝑡+𝜑𝑝→ℎ ) is the infinitesimal unsteady harmonic pressure and ds is area of surface element, and n is surface normal. For the 2D case 336

𝑑 𝑓̂𝜉 = 𝑐𝑝,𝐴 ̂ .𝑛⃗𝜉 . 𝑑𝑠

(26)

𝑑 𝑓̂𝜂 = 𝑐𝑝,𝐴 ̂ .𝑛⃗𝜂 . 𝑑𝑠

(27)

𝑑 𝑚̂ 𝜁 = (𝑟⃗ × 𝑐𝑝,𝐴 ̂ ).𝑒⃗𝜁 . 𝑑𝑠

(28)

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 11. Indexing of blades in 11 blade cascade configuration and two different mode of cascade vibration TWM and INFC.

where 𝑛⃗ local normal vector pointing outwards on the surface. In 2D flow analysis the third component Eq. (28) is represent the moment rather than force. The work per oscillation cycle (Wcycle ) can obtained by integrating the product of force and motion during the cycle (T) and can be given by 𝑊𝑐 𝑦𝑐 𝑙𝑒 =

∫𝑇

̃ ⃗̃ 𝐹⃗ .ℎ. 𝑑𝑡 =

∫𝑇

̂ ⃗̂ 𝑖𝜔𝑡 𝐹⃗ .ℎ.𝑒 𝑑𝑡

The positive value of “Ξ” indicates the flow acting in stabilizing manner whereas negative value can cause flutter, therefore, unstable. 5. Test case for classical flutter analysis using unsteady 2D PM In this section application of PM based solver to simulate the classical flutter stability parameters on two test cases is described. To construct the flow field around oscillating cascade the similar discretization scheme is adopted as described in Section 3.1. For the test case 11 blade cascade with NACA65 series airfoil profile is selected and flow condition are similar to the Carta [16]. A brief description about experimental setup and cascade geometry is presented in the following section.

(29)

⃗̂ is the complex motion of the blade. In the present case only where ℎ pitching motion of the cascade blade is analyzed therefore, it is corresponds to torsion flutter, therefore, aerodynamic moment Eq. (28) is used, After integration for the moment the Eq. (29) can be written more explicitly as 𝑊𝑐 𝑦𝑐 𝑙𝑒 = 𝜋|𝑎𝜁 |.|𝑚𝜁 |. sin(𝜑𝑚𝜁 →𝜑𝜁 )

5.1. Experimental setup description

(30)

where a𝜁 amplitude of oscillation, m𝜁 moment around torsional axis and 𝜑𝜁 IBPA of torsional flutter. Form the Eq. (30) it is oblivious that the only imaginary part of the force term contribute in the work, therefore, if the response is lagging the excitation, the imaginary part gives negative contribution into flow, consequently has stabilizing effect. More detail about calculation of work per oscillation cycle (Wcycle ) is presented in [19]. Therefore, aerodynamic damping parameter is given by work done per cycle and normalized by 𝜋 and amplitude of oscillation “A” (for bending flutter A is equivalent to “h” and for torsional flutter A is equivalent to “𝜔” and is given in Eq. (31) 𝑊𝑐 𝑦𝑐 𝑙𝑒 𝐴𝑒𝑟𝑜𝑑 𝑦𝑛𝑎𝑚𝑖𝑐 𝑑 𝑎𝑚𝑝𝑖𝑛𝑔(Ξ) = − (31) 𝜋𝐴2

The experiment on the oscillating cascade is carried out by Fraklin O. Carta in the UTRC linear subsonic Oscillating Cascade Wind Tunnel (OCWT). The schematic diagram of the OCWT is presented the Fig. 9. The test section dimension of the wind tunnel is 25.4 cm wide and 68.6 cm high. The wind tunnel is facilitated with 11 shaft-mounted blade cascade with sidewall configuration for 2D flow simulation. With the shaft mounted blades it is possible to achieve different IBPA between to simulated TWM oscillation. The experimental cascade’s geometrical configuration is presented in the Fig. 10. The test cascade consist of 11 NACA 65-series blades, each blade have chord length c = 15.24 cm, and span 24.4 cm, with a 135∘ circular arc camber and a thickness-tochord ratio of 0.06. The 11 blades are arranged in 30∘ staggered angle 337

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 12. Unsteady wake vortices of the cascade at 𝛼0 = 6◦ and 𝛼̄ = 2◦ at reduce frequency 𝑘 = 0.122 and wind speed 𝑄∞ = 61.2 m∕s, 𝜑 = −45◦ using PM model.

Fig. 13. The AD (Ξ) vs IBPA, Simulated (PM) with Exp. at 𝛼0 = 6◦ and 𝛼̄ = 0.5◦ at reduce frequency 𝑘 = 0.122 and wind speed 𝑄∞ = 61.0 m∕s .

(𝛼 eff ) at any time instance is given by Eq. (33)

with slant height 11.43 cm. The entire set of airfoils can be coherently driven in a sinusoidal pitching motion with an amplitude of 𝛼. ̄ The system is driven by a constant speed, 5.6 kW (7.5 hp) electric motor with a continuously variable speed transmission through a timing belt and pulleys. There are two sets of test cases at two different reduced frequency (𝑘 = 0.122 and 0.151 based on semi chord) values are selected for comparison. All the other flow parameters remain same for both test cases which have constant inlet velocity of 61 m/set (200 ft/s), 𝑅𝑒 = 6.15 × 105 , mean camber line incidence angles (𝛼0 = 6◦ ), pitching amplitudes (𝛼̄ = 0.5◦ ), and 8 IBPA (𝜑) = 0, ± 45∘ , ± 90∘ , ± 135∘ , ± 180∘ . More details about the experimental setup and data acquisition is presented in [16]. Furthermore, the blade indexing strategy in numerical model and schematic representation of TWM and the INFC mode oscillation is presented in the Fig. 11. The oscillatory motion of the cascade blade can be given in complex form by Eq. (32) 𝛼(𝑡) = 𝛼𝑒 ̄ 𝑖(𝜔𝑡±𝑁𝜑)

𝛼𝑒𝑓 𝑓 (𝑡) = 𝛼0 + 𝛼(𝑡) = 𝛼0 + 𝛼𝑒 ̄ 𝑖(𝜔𝑡±𝑁𝜑) = 𝛼0 + 𝛼̄ sin(𝜔𝑡 ± 𝑁𝜑)

(33)

5.2. Results and discussion The aerodynamic damping curve w.r.t IBPA is estimated using PM solver for comparison purpose. The simulated aerodynamic damping (Ξ) for 𝛼0 = 6◦ and 𝛼̄ = 0.5◦ at two reduce frequency 𝑘 = 0.122, 0.151 and wind speed 𝑄∞ = 61 m/s, 𝑅𝑒 = 6.15 × 105 (based on reference chord of the blade), 𝑀𝑎𝑐ℎ = 0.178 is compared with the experimental result as presented by Carta [16]. In the simulation IBPA is −180◦ to +180◦ with increment of 10∘ interval is set. The results are presented in the Figs. 13 and 14. In the Fig. 12 simulated unsteady wake shape for the 11 blade cascade at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122 and wind speed 𝑄∞ = 61.2 m∕s, 𝜑 = −45◦ is presented. This 3rd test case have been used to CFD results comparison in the next section. In the Fig. 12 the simulated wake shapes looks reasonable with wavy pattern caused by the

(32)

where N is the blade number (Fig. 11) and the 𝜔 is the frequency of 2.𝑘.𝑄 oscillation in rad/sec and 𝜔 = 𝑐ℎ𝑜𝑟𝑑∞ . Hence, the effective angle of attack 338

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 14. The AD (Ξ) vs IBPA, Simulated (PM) with Exp. at 𝛼0 = 6◦ and 𝛼̄ = 0.5◦ at reduce frequency 𝑘 = 0.151 and wind speed 𝑄∞ = 61.0 m∕s . Fig. 15. Grid strategy for CFD numerical model of the 11 blade cascade .

timated the aerodynamic damping magnitude beyond > +60◦ IBPA. The main cause of this overshoot in amplitude can be caused by disagreement in phase lead or acoustic resonance at those IBPA [12], however, it is to be investigated in more details. Fig. 14 which is AD curve for higher value of k = 0.151 has also good agreement with experimental result but similar to the previous test case of k = 0.122, there is overestimation of the aerodynamic damping magnitude for > +60◦ IBPA. The effect of phase lead or acoustic resonance is not accounted in the flow model. Moreover, there is not significant change in amplitude of aerodynamic

pitching motion of the blades and also wake rollup at the end is as expected. Only 100 time step wake shape is presented in the figure for better visibility, however, in calculation wake length up to 20 chord length is considered which is corresponds to 300 time steps. In the Fig. 13 simulated aerodynamic damping shows good agreement with experimental result, and simulated result captures the stable and unstable zone with great accuracy and the stability limit (Ξ = 0) is well predicted and it is highlighted in blue circle in both figures, however, in the unstable zone the PM based flow model has slightly overes339

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 16. Grid strategy and BC for CFD- inviscid numerical model of the cascade.

damping observed with increase of reduce frequency from k = 0.122 to 0.151, in both experimental and simulated cases. Note worthy that the time of execution for PM based model to simulate each of these test cases is less than 15 CPU minutes on single processor (processor Intel i7 speed 3.75 GHz) significantly less CPU time than conventional time-domain NS-based CFD models. This can be one of the main key advantage of PM based solver over CFD models for classical flutter stability analysis.

6.2. CFD simulation for the 2D cascade and grid strategy The entire flow field is discretized using triangular FV elements. For the oscillating cascade flow simulation it is very crucial to use proper grid with adequate refinement in order to capture unsteady aerodynamics forces with high accuracy. Therefore, the entire domain is divided into three parts, 1) Boundary layer close to the airfoil wall 2) non- deforming zone in the proximity of the airfoil (surrounded within blue line Fig. 15). 3) outer deformable zone for adoptive meshing. The grid strategy for the present case is shown in the Fig. 15, a similar grid strategy is adopted by Vimmr et al. [20] to estimate the AD for flat plate cascade. The actual size of entire domain is large, for the image clarity reason only the mesh around the cascade in shown in the Fig. 15.

6. Estimation of aerodynamics damping using NS based CFD method In this section the Navier Stokes (N-S) based CFD methods is employed to estimate the AD of the oscillating cascade configuration as experimented by Carta [16]. The main purpose of this study is to evaluate the efficiency and accuracy of the PM based model compared to classical CFD based model. The simulated aerodynamic damping (Ξ) for 𝛼0 = 6◦ and 𝛼̄ = 2◦ at reduce frequency 𝑘 = 0.122 and wind speed 𝑄∞ = 61.2 m∕s, 𝑅𝑒 = 6.3 × 105 , 𝑀𝑎𝑐ℎ = 0.178 is compared with the experimental result as presented by Carta [16]. In the CFD simulation 8 IBPA (𝜑) = 0, ± 45∘ , ± 90∘ , ± 135∘ , ± 180∘ are simulated to get the stability s-curve, whereas, in PM based model IBPA −180◦ to +180◦ with increment of 10∘ is used. The brief produce for CFD simulation is presented in the following subsections.

6.2.1. Grid convergence test on CFD mesh A grid convergence study similar to that of PM method case has also been carried out on CFD mesh. For the CFD convergence steady state inviscid case is simulated. The lift coefficient of the reference (0) blade at 𝛼0 = 6◦ and wind speed 𝑄∞ = 61.2 m∕s is estimate for different mesh density. The graph of lift coefficient and CPU time (h) vs grid size is present in the Fig. 16. It can be noticed in the Fig. 16 the lift coefficient start to converge around 0.5 million mesh size. Therefore, 0.55 million mesh size is used in the simulation as optimal size w.r.t accuracy and CPU time. During the mesh refinement special attention is paid at the leading and trailing edge of the airfoil and very fine grid size is used in those region Fig. 15. The global mesh specification is given in the Table 1. Once the boundary conditions are defined in ANSYS/Fluent, the simulation can be launched. The oscillatory motion of the blade is defined by user define function (UDF) in ANSYS/Fluent. In the first step steady state simulation is carried out, once the convergence achieved, the transient simulation with time marching procedure is carried out. The time step size is equal to the one 100th of period of oscillation (i.e 0.000645 s), the transient simulation is carried out for 300 time steps. The unsteady pressure and the aerodynamics moment coefficient (Cm)

6.1. CFD numerical model description For the CFD numerical modeling purpose of blade cascade inviscid numerical model is adopted here, since PM is inviscid method, therefore, it is more appropriate to compare it with inviscid CFD model for classical flutter stability analysis case. The numerical implementation of the CFD fluid model based on the finite volume (FV) discretization of inviscid version of Navier-Stokes equation (or Euler equation) of compressible flow. The CFD numerical analysis is carried out using commercial CFD code ANSYS/Fluent 17.2. 340

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 17. The AD (Ξ) vs IBPA, Simulated (PM & CFD) against Exp. at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122 and wind speed 𝑄∞ = 61.2 m∕s.

Table 1 CFD grid specifications. Parameters

Values

Discretization type Element type Number of Elements Fluid type 𝑌 + value

FV 2D Triangular 0.55 million Ideal gas ≈ 2.7

which agrees well with experimental data. Furthermore, in the velocity contour plot at different time step in Fig. 22 no significant flow separation is observed on suction side of the cascade airfoils and the flow remains attached through out oscillation period, thus the PM based result are arguably acceptable for this flow condition which is the case of classical flutter. However, the simulated AD values are slightly underestimated (highlighted by the circle) for the stability limit at 𝜑 = 0◦ compare to experimental data, but still remain positive. Furthermore, in both the simulated results the peak magnitude values of AD have poor agreement with experimental data. This is mainly caused by the disagreement of moment coefficient and also due to disparagement in phase lead at those IBPA (𝜑) [12]. The other reason for scattering in the magnitude of AD around maximum and minimum values in the cascade in Fig. 17 is caused by acoustic resonance [12], which is not considered in any of the numerical model. Moreover, the discrepancy in magnitude of AD of the simulated results at maximum and minimum values can also be caused by flow separation on suction side due to higher amplitude oscillation in this case 2∘ which is way higher than in the previous test case 0.5∘ at 𝑘 = 0.122, and for the same test case PM results are more close to the experimental results because of complete attached flow. Furthermore, the simulated AD value using inviscid CFD model is very close to PM model’s estimations, this is because of invscid nature of the both solvers. From the results in Fig. 17 it can be said that the CFD based models are also failed to estimate the AD very accurately at maximum and minimum points, however, the CFD results are slightly better compare to PM. In the Figs. 18 and 19, the variation of real and imaginary part of time averaged ΔCp on reference blade is compared with PM model and CFD model. The simulated results shows good agreement with ex-

around the mid of the reference blade saved as main results of simulation. The AD coefficient is estimated by multiplying Cm for one complete cycle by blade motion in one complete cycle. The AD is normalized by dividing with pitching amplitude (𝛼̄ = 2◦ ). The CFD simulation in carried out in both viscous and inviscid models. 6.3. Results and discussion The AD and the difference of pressure coefficient (ΔCp ) between suction and pressure side on the reference blade is compared with the experimental data from Carta [16]. The AD estimated by CFD inviscid model and PM for the 𝛼0 = 6◦ and 𝛼̄ = 2◦ at reduce frequency 𝑘 = 0.122 and wind speed 𝑈∞ = 61.2 m∕s is compared against experimental value and presented in Fig. 17 and variation of real (ΔCp (R)) in Fig. 18 and imaginary part (ΔCp (Im)) in Fig. 19 In the Fig. 17 estimated AD using PM and CFD inviscid model have good agreement with the experimental data. Both PM and CFD inviscid model can capture the typical s-curve of AD vs IBPA, and the prediction of AD for 𝜑 = 0◦ in stable zone boundary 341

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 18. Variation of Real ΔCp , against Exp. at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122, wind speed 𝑄∞ = 61.2 m∕s and 𝜑 = −45◦ .

Fig. 19. Variation of Img ΔCp , against Exp. at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122, wind speed 𝑄∞ = 61.2 m∕s and 𝜑 = −45◦ .

342

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 20. Comparison of steady state Cp at reference blade with Exp. at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122, wind speed 𝑄∞ = 61.2 m∕s and 𝜑 = −45◦ .

Table 2 Numerical model execution time comparison.

perimental value, however, there is jump in simulated ΔCp close to the leading edge peak values, but is hard to very the accuracy at this point due to lack of experimental data at those points. In the Fig. 20 the simulated steady state pressure coefficient distribution over reference blade for the same flow condition is presented and the results show good agreement with experimental data, however, the overestimation of Cp at the leading edge by PM and CFD inviscid model can not be verified due to unavailability of measured data at very leading edge, it is difficult to install measuring sensors at these positions. The steady state pressure contour (CFD-inviscid) profile in all 11 blades is presented in the Fig. 21 and there is no flow separation is observed on the reference blade and the immediate 4 neighbouring blades. In TWM effect of neighbouring blades −2 to +2 are most significant [12] and the aerodynamic coupling effect become negligible as we move further. In the Fig. 22 velocity profile for different time steps of one complete cycle are presented and the time instance are indicated by the solid circle on the sin curve of oscillation in the sub-figure (a, b and c) of Fig. 22, only velocity contour around −2, −1, 0, +1, +2 is presented for better clarity. The change in velocity contour with blade position can be observed in the three sub-figure (a, b and c) of Fig. 22, the difference in not very large due lower amplitude of oscillation, it is more clear in the animation. In the sub-figure (d) the velocity profile a 0 time instance around all 11 blade is presented, and there is not significant flow separation is predicted by CFD-inviscid model. The execution time comparison between PM model and CFD inviscid model to estimate complete the scurve (in Fig. 17) for 11 blade cascade model is presented in the Table 2. Once again, since PM model is a inviscid method, therefore, it is more reasonable to compare it with CFD inviscid model. It is clear in the table that the PM method have big advantage over CFD inviscid model in the term of execution time, yet the CFD prepossessing and post prepossessing time in not included in table, which can be make over all time for

Model

No. processor

CPU time

IBPA interval (−180◦ to +180◦ )

PM CFD-Inviscid

1-core, Intel i7 4-core, Intel i7

15 min 13 h

36 (with 10∘ increment.) 9 (with 30∘ increment)

CFD mode even higher. It worth mentioning that generally CFD viscous simulation takes even longer time compare to inviscid simulations. The close observation from above results, it is clear that the both PM and CFD can estimate the AD and pressure coefficient with good accuracy, however, CFD result are more close to experimental value compared to PM based model in large amplitude case, but the computational time for CFD model significantly higher than PM based model. Therefore, PM method clearly have advantage over CFD on computation which makes it good compromise of speed and accuracy over CFD model for various industrial or academic application where flow field can be assumed as potential flow. 7. Conclusion The paper explored the issue of modeling unsteady flow around 2D cascade and estimation of aerodynamics damping for classical flutter case using computationally less expensive boundary element methods e.g. Panel method. The aerodynamic damping and pressure coefficient for 11 blade 2D cascade have been successfully estimated using medium fidelity potential flow based PM model and the results are compared against experimental as well as CFD-inviscid models. The careful review of the AD s-curve the Figs. 13, 14 and 17 and the pressure coefficient results, shows that for classical flutter or low-amplitude vibration of blade cascade, the flow field can be modeled with reasonable accuracy 343

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

Fig. 21. Steady state Pressure contour at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122, wind speed 𝑄∞ = 61.2 m∕s and 𝜑 = −45◦ for CFDInviscid numerical model.

Fig. 22. Velocity contour at three different time at 𝛼0 = 6◦ , 𝛼̄ = 2◦ , 𝑘 = 0.122, wind speed 𝑄∞ = 61.2 m∕s and 𝜑 = −45◦ of CFD-Inviscid numerical model of the cascade.

using PM based boundary element method. However, the peak magnitude values of AD damping s-curve show poor estimation by PM model compared to experimental result due to influence of acoustic resonance and phase lead at peak IBPA in experimental data. Nevertheless, the accuracy of the present versions of PM can further be improved by incorporating acoustic resonance model. This can be an interesting topic of research in this field in future. Further, for a large amplitude or stall flutter cases the present version of PM is not suitable and it have been shown in Fig. 6 for the deep stall pitching airfoil cyclic lift simulation, for such condition CFD models will continue to dominate. Now at this point in the light of the research carried out in this paper it can be said that, in the turbomachinery cascade if the flow field is incompressible and free from significant flow separation, boundary element methods e.g.

PM model, can be the best option over CFD models to model the flow field and to estimate the aerodynamics/aeroelastic parameters, due to following key characteristics, • Computationally less expensive compared to CFD models. • The accuracy of PM based BEM solver’s results are well within acceptable range for low amplitude oscillation and until the flow field resembles the underlying principles of potential flow, keeping it in mind as medium fidelity method. Therefore, the present method can be applied to estimate the aeroelastic stability parameters for low frequency and low amplitude oscillation cases e.g. classical flutter analysis. Further, more advance version of PM with flow separation model including viscous effect can be 344

C.S. Prasad and L. Pešek

Engineering Analysis with Boundary Elements 113 (2020) 328–345

implemented to further increase the accuracy of the results including the effect of acoustic resonance model. As a over all conclusion it can be said that the present PM based aeroelastic model for cascade flow is a versatile design tool and it can have significant impact on the computational time reduction and provide researchers and engineering in this field a freedom to iterate different blade profile in very short period of time at preliminary design stage of LP turbine and compressor design and for blade shape optimize. However, the use of CFD based aeroelastic model are inevitable for the flow conditions where significant amount of flow non-linearity can exist e.g. in transonic flows, separated flows or large amplitude high frequency flow in the cascade (stall flutter case).

[6] Prasad C, Chen Q-Z, Bruls O, D’Ambrosio F, Dimitriadis G. Advanced aeroservoelastic modeling for horizontal axis wind turbines. In: Proceedings of the 9th international conference on structural dynamics, EURODYN 2014, Porto, Portugal; 2014. p. 3097–104. [7] Prasad C, Chen Q-Z, Bruls O, D’Ambrosio F, Dimitriadis G. Aeroservoelastic simulations for horizontal axis wind turbines. Proc Inst Mech Eng Part A 2017;231(2):103–17. [8] Katz J, Plotkin A. Low-speed aerodynamics. 2nd ed. New York, USA: Cambridge University Press; 2001. [9] McFariand E. Solution of plane cascade flow using improved surface singularity methods. J Eng Power July 1982;104:669. [10] Vahala LL, Lei ZQ, Liang GZ. Solution of turbine blade cascade flow using an improved panel method. Int J Aerosp Eng 2015;2015 6 pages. Article ID: 312430 [11] Barbarossa F, Parry AB, Green JS, di Mare L. An aerodynamic parameter for low– pressure turbine flutter. J Turbomach 2016;138(5):51001. [12] Fransson T. Aeroelasticity in axial flow turbomachines. Brussels, Belgium: Von Karman Institute for Fluid Dynamics; 1999. [13] Verdon JM. Review of unsteady aerodynamic methods for turbomachinery aeroelastic and aeroacoustic applications. AIAA J 1993;31(2):235–50. [14] Vogel K, Naidu A, Fischer M. Comparison of the influence coefficient method and travelling wave mode approach for the calculation of aerodynamic damping of centrifugal compressors and axial turbines. In: ASME turbo expo 2017: turbomachinery technical conference and exposition. American Society of Mechanical Engineers; 2017. V07BT36A022–V07BT36A022 [15] Camara E. Validation of time domain flutter predictiontool with experimental results. 2015. [16] Carta FO. An Experimental Investigation of Gapwise Periodicity and Unsteady Aerodynamic Response in an Oscillating Cascade I -Experimental and Theoretical Results. Tech. Rep. NASA; 1982. [17] McCroskey WJ, McAlister K, Carr L, Pucci S, Lambert O, Indergrand R. Dynamic stall on advanced airfoil sections. J Am Helicopter Soc 1981;26(3):40–50. [18] Hansen YC. An introduction to the theory of aeroelasticity. 7th ed. Mineola, New York: Dover Publication Inc; 2008. [19] Vogt D. Experimental Investigation of Three-Dimensional Mechanisms in Low-Pressure Turbine Flutter; 2005. Ph.D. Thesis. [20] Vimmr J, Bublík O, Pecka A, Pešek L, Procházka P. Numerical and experimental study of fluid flow in simplified blade cascade with prescribed harmonic motion. In: EPJ web of conferences, 180. EDP Sciences; 2018. p. 02116.

Acknowledgment This research is supported by the research project of Czech Science Foundation no. 16-04546S “Aeroelastic couplings and dynamic behavior of rotational periodic bodies”. References [1] Atassi H, Akai T. Effect of blade loading and thickness on the aerodynamics of oscillating cascades. In: 16th aerospace sciences meeting; 1978. p. 277. [2] Panovsky J, Kielb R. A design method to prevent low pressure turbine blade flutter. J Eng Gas Turbine Power 2000;122(1):89–98. [3] Fransson T. Udine: lecture notes/introduction to blade flutter in axial flow turbomachinery. 1993. [4] Hess JL. Calculation of potential flow about arbitrary three-dimensional lifting bodies. Final Technical Report MDC J5679-01. Naval Air Systems Command, Department of the Navy; 1972. [5] Prasad CS, Dimitriadis G. Application of a 3D unsteady surface panel method with flow separation model to horizontal axis wind turbines. J Wind Eng Ind Aerodyn 2017;166:74–89.

345