Milind A. Bakhle

Glenn Research Center, Cleveland, Ohio

T.S.R. Reddy

University of Toledo, Toledo, Ohio

Harmonic Balance Computations of

Fan Aeroelastic Stability

NASA/TM—2010-216222

August 2010

IFASD–2009–101
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Milind A. Bakhle

Glenn Research Center, Cleveland, Ohio

T.S.R. Reddy

University of Toledo, Toledo, Ohio

Harmonic Balance Computations of

Fan Aeroelastic Stability

NASA/TM—2010-216222

August 2010

IFASD–2009–101

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

International Forum on Aeroelasticity and Structural Dynamics (IFASD 2009)

cosponsored by CEAS and AIAA

Seattle, Washington, June 21–25, 2009
Acknowledgments

Support for this work from the NASA Fundamental Aeronautics Program is gratefully acknowledged. The authors would like 

to thank Dr. Nateri Madavan and Dr. James Heidmann (Associate Principal Investigators for the Subsonic Fixed Wing Project) 

for their support of this work. The authors would also like to thank Mr. Brian Fite for providing us the opportunity to work with 

this confi guration. The authors would like to acknowledge Dr. Jeffrey Thomas of Duke University for his support of this work.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216222 1

Harmonic Balance Computations of Fan Aeroelastic Stability

Milind A. Bakhle 

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

T.S.R. Reddy

University of Toledo

Toledo, Ohio 43606

Abstract

A harmonic balance (HB) aeroelastic analysis, which has been recently developed, was used to 

determine the aeroelastic stability (flutter) characteristics of an experimental fan. To assess the numerical 

accuracy of this HB aeroelastic analysis, a time-domain aeroelastic analysis was also used to determine 

the aeroelastic stability characteristics of the same fan. Both of these three-dimensional analysis codes 

model the unsteady flowfield due to blade vibrations using the Reynolds-averaged Navier-Stokes (RANS)

equations. In the HB analysis, the unsteady flow equations are converted to a HB form and solved using a 

pseudo-time marching method. In the time-domain analysis, the unsteady flow equations are solved using 

an implicit time-marching approach. Steady and unsteady computations for two vibration modes were 

carried out at two rotational speeds: 100 percent (design) and 70 percent (part-speed). The steady and 

unsteady results obtained from the two analysis methods compare well, thus verifying the recently 

developed HB aeroelastic analysis. Based on the results, the experimental fan was found to have no 

aeroelastic instability (flutter) at the conditions examined in this study.

1.0 Introduction

Aircraft engine fan and compressor blades are susceptible to flutter, which can lead to catastrophic 

failure. Flutter problems will likely be encountered more often in newer engine designs that are being 

developed to reduce noise, improve performance, and reduce weight—using thinner blades, more 

complex blade geometry, and higher blade loading. To avoid unexpected flutter during engine 

development and testing, it is important to develop and validate numerical tools that can be used to verify 

aeroelastic stability using high-fidelity physics-based models. Research has been on-going in the 

development, validation and application of high fidelity models for aeroelastic vibrations in aircraft 

engine fan, compressor and turbine blades (Ref. 1). Recent work has included time-domain solution of the 

Reynolds-averaged Navier-Stokes (RANS) equations to provide the unsteady flowfield and unsteady 

aerodynamic loads on vibrating blades. An example of such work is the TURBO aeroelastic analysis code 

(Refs. 2 and 3), which has predicted aeroelastic stability characteristics that agreed well with experimental 

data. However, such time-domain solutions require large computational times due to the long startup 

transients that must decay before solutions become periodic. To reduce computational cost while 

maintaining the same level of prediction accuracy, the harmonic balance (HB) method has recently been 

developed and used (Ref. 4). In this approach, each flow variable is represented by a Fourier series in 

time, leading to a HB form of the RANS equations. The resulting equations can be solved using methods 

that are typically used for steady flow problems such as pseudo-time marching and local time stepping. 

The HB approach may thus provide a high fidelity modeling capability with a reduced computational cost

compared to the typical time-domain solution method.

In the present study, the aeroelastic stability (flutter) characteristics of an experimental fan are 

analyzed using a HB RANS aeroelastic analysis and computer code developed by Hall et al. (Refs. 4 and 

5). Since the focus here is to verify the recently developed HB aeroelastic code by direct comparison, a 

time-domain RANS aeroelastic code (Ref. 2) (TURBO) is also used to analyze the same experimental 
NASA/TM—2010-216222 2

fan. The objective is to validate the accurate computation of unsteady aerodynamics due to vibrating 

blades. The experimental fan configuration used in this study is described in detail in Reference 6. This 

fan was designed to demonstrate a noise-reduction technology based on trailing edge flow injection. Note, 

however, that in this paper, the fan is analyzed only in the mode of operation with no trailing edge flow 

injection.

Both analyses used in this study solve the three-dimensional unsteady, RANS equations with the 

ability to model a rotating blade row with harmonic blade vibrations or incoming periodic distortions. For 

flutter calculations, the blade vibration is prescribed to be the modal deflection and a frequency, both of 

which are calculated from a separate structural dynamics analysis. Computations are performed in a 

single blade passage for both steady and unsteady analyses. In the HB analysis (Refs. 4 and 5), the 

resulting unsteady pressures on the blade are used to calculate a generalized aerodynamic force and then 

an eigenvalue problem is solved to calculate the aerodynamic damping. The damping value is needed to 

determine the aeroelastic stability in the selected vibration mode and traveling wave pattern (Ref. 7). In 

the time-domain analysis (Ref. 2), the unsteady pressures are used to calculate a work-per-cycle that is 

then used to calculate the aerodynamic damping, which is then used to determine aeroelastic stability. 

This paper describes the aeroelastic stability calculations and compares the results from the HB and timedomain analyses. Results are presented at design speed and part-speed for aeroelastic analysis with the 

first and second vibration modes, followed by conclusions.

2.0 Analyses Methods

2.1 Aeroelastic Model

The equations of motion for a fan blade (with all blades assumed to be identical) can be written as

[M ]{q}+ [K]{q}= [A]{q} (1)

where [M] and [K] are generalized mass and stiffness matrices, {q} is the generalized displacement 

vector, and [A] is the blade vibration-dependent generalized aerodynamic force matrix. The matrices [M], 

[K], and [A] are of size nm×nm; {q} is of size nm×1; nm is the number of modes.

The elements of [M] and [K] are obtained from a free-vibration analysis using commercial structural 

dynamics analysis software. The matrices [M] and [K] are diagonal and their non-zero elements are 

related as

K = M ω ( + iζ) i i i 1 2 2 (2)

where ωi is the natural frequency of the i

th mode, and ζ is the structural damping ratio; usually the mode 

shapes are mass-normalized and therefore Mi = 1.

Since all the blades are identical (that is, a tuned rotor), the aeroelastic modes consist of individual 

blades vibrating with equal amplitudes at a fixed interblade phase angle between adjacent blades. Hence, 

the motion of a blade in r

th interblade phase angle mode can be written as

{ } { } i t i s r r q q e e ω σ = (3)

where ω is the vibration frequency, σr is the interblade phase angle related to nodal diameter (ND) pattern 

of the traveling wave and number of blades Nblades as

blades, blades blades 2 / / 2 / 2 σ=π r ND N N ND N − ≤ ≤ (4)
NASA/TM—2010-216222 3

Thus, the equations of motion for a blade become

−ω[M ]{qr}+ [K]{qr}= [Ar ]{qr} (5)

The following subsections describe the flutter calculation method used with the HB code and separately 

with the time-domain TURBO code.

2.2 Harmonic Balance Flutter Analysis

The HB flutter analysis requires calculation of elements of the generalized aerodynamic force matrix 

[Ar]. Unsteady flowfield computations are carried out for each vibration mode and assumed frequency. 

For a selected value of interblade phase angle, the HB code is used to calculate the unsteady pressure 

distribution on the blade surface, which is further used to calculate the (complex-valued) elements of the 

generalized aerodynamic force matrix [Ar]. This calculation is repeated for Nblades interblade phase angles 

given by Equation (4).

To calculate flutter stability, Equation (5) is written in a standard eigenvalue form as:

[[P]− γ[Q]]{qr}= {0} (6)

where [ ] ([ ] [ ]) [ ] [ ] ( ) 2 2 P K Ar o Q M = ω ωo = − / ω ; = ; γ / and ωo is the assumed frequency used in the 

calculation of the elements of the aerodynamic force matrix [Ar].

The solution of the above eigenvalue problem results in nm complex eigenvalues of the form

i(ω ωo )= i γ= μ ± iν (7)

The real part of the eigenvalue ( μ ) represents the damping ratio, and the imaginary part ( v ) represents 

the damped frequency; flutter occurs if μ ≥ 0 for any eigenvalue.

In the present work, structural damping is set to zero and therefore damping is referred to as 

aerodynamic damping. Also, note that μ is opposite in sign to damping.

2.3 TURBO Work-Per-Cycle Method

The work-per-cycle approach is used in the TURBO analysis to determine flutter stability. First, the 

flowfield through the blade row is calculated with no prescribed blade vibration. Starting with this 

converged steady flowfield, blade vibrations are prescribed in a selected mode, frequency, and nodal 

diameter pattern or phase angle. After the transients in the flowfield decay, and a periodic flowfield is 

obtained, the work done on the vibrating blade is calculated for a cycle of blade vibration as follows: 

W= − pdA⋅(∂X ∂t)dt ∫ ∫

surface

(8)

where p is the blade surface pressure, A is the surface area vector, X is the position vector, and

t is time. For harmonic vibration, the work-per-cycle of oscillation, can be re-written as

W= − pdA⋅δq ω (ωt)dt ∫ ∫ cos 0

surface

(9)
NASA/TM—2010-216222 4

where δ is the modal displacement vector, q0 is the amplitude of the generalized displacement, and ω is 

the vibration frequency. The aerodynamic damping ratio (ζ ) associated with blade vibration is related to 

the work-per-cycle (W) and the average kinetic energy (KE) of the blade over one cycle of vibration 

through the following equation from Reference 8

2 1

8

− ζ

πζ = − KE

W (10)

where

∫ = mV dt

T

KE

2

2

1 1 (11)

=C Ccr ζ ; Ccr = 2mω (12)

In the preceding equations, C is the damping, Ccr is the critical damping, m is the mass of the

blade, V is the surface velocity due to blade vibration, and T is the time period. 

For small values of damping ratio which typically occur in aeroelastic calculations of interest,

ζ << , the aerodynamic damping ratio can be approximated as 

≈−W πKE ζ 8 (13)

If aerodynamic damping is negative, flutter can occur. Note that the structural damping (material and 

mechanical damping) has not been considered. Also, note that aerodynamic damping (ζ ) is opposite in 

sign to μ .

3.0 Results

In this section, results of steady and unsteady computations are presented for an experimental fan 

described in Reference 6. This fan has 18 blades with a tip diameter of 22 in.; the design parameters are 

rotational speed of 8750 rpm, pressure ratio of 1.284, tip speed of 840 ft/s, and mass flow rate of

85.7 lbm/s. Steady and unsteady computations were carried out at two rotational speeds: 100 percent

(design) and 70 percent (part-speed). The latter rotational speed was selected to investigate, generically, 

the possibility of part-speed flutter for this fan, with the understanding that multiple speeds must be 

examined in practice to establish flutter stability of the fan. Unsteady computations were carried out for 

the first two structural vibration modes (bending and torsion).

Note that in the present work, the HB aeroelastic code was run with very small amplitude vibrations 

of the fan blade to calculate linearized unsteady aerodynamics for a conventional linear flutter analysis. 

Also, the amplitude of vibration for the time-domain (TURBO) computations was in a linear range.

3.1 Steady Computational Results

The airfoil geometry used in this study was calculated based on static deflections obtained from 

structural analysis at design speed. The static deflections included the effects of rotational speed and 

aerodynamic pressure at design condition; nonlinear geometry effects were included in the structural 

analyses. The same blade geometry was used for both design and part-speed computations. Previous 

calculations for a different fan (Ref. 2) have shown that changes from the nominal blade geometry due to 

changes in rotational speed were not significant and therefore a single geometry was used in the present 

work for computations at both speeds and all operating conditions.
NASA/TM—2010-216222 5

The computational grids used here were generated using commercial software. The grid for the HB

computations is shown in Figure 1(a); the grid size is 193 by 33 by 49 for the O-grid block that wraps

around the blade airfoil with 193 points around the airfoil, 33 points in the circumferential direction, and 

49 points in the spanwise direction. The H-grid blocks in the inlet and exit sections are each 17 by 33 by 

49 with 17 points in the streamwise direction, 33 grid points in the circumferential direction, and 49 grid 

points in the spanwise direction. The tip clearance present in the grid was based on design intent and is 

modeled using 9 points between the blade tip and casing. For the TURBO time-domain computations, the 

H-grid, shown in Figure 1(b), has a size of 111 by 41 by 39 with 111 points in the axial direction,

41 points in the spanwise direction, and 39 points in the circumferential direction; 4 cells are used in the 

tip clearance region.

In the computations, inlet flow conditions consisted of uniform total pressure and total temperature 

(uniform inflow). The exit flow conditions consisted of a static pressure value prescribed at the hub and 

calculated at other radial locations using radial equilibrium.

Azimuthal view Meridional view

Figure 1(a).—Computational mesh for HB computations.

Azimuthal view Meridional view

Figure 1(b).—Computational mesh for time-domain TURBO computations.
NASA/TM—2010-216222 6

Figure 2.—Steady solution convergence from HB 

method.

Figure 3.—Fan map showing performance predictions 

from HB and time-domain codes.

Steady computations were carried out for two rotational speeds (100 and 70 percent). The exit static 

pressure on the hub was varied to change the mass flow rate and the operating point on the fan map. Thus, 

two speed lines were generated using each analysis code. Figure 2 shows the convergence of the steady 

HB computations as a plot of non-dimensional torque with iteration counter (100 iterations per counter). 

These computations were carried out at 100 percent speed for various values of imposed back-pressure at 

the exit boundary. As can be noted from Figure 2, excellent convergence was obtained for all operating 

conditions. Similar convergence was obtained for the computations at 70 percent speed as well as for the 

computations with the TURBO code at both speeds. All steady results presented here are from wellconverged solutions.

The fan map, total pressure ratio versus mass flow rate, is shown in Figure 3 with the results of HB

steady computations denoted as HB and the results of the time-domain computations denoted as TURBO. 

Results are presented for computations at several operating conditions for 100 and 70 percent speeds. The 

computed results from the two analysis codes correlated well at both speeds. At 100 percent speed, the 

difference in pressure ratio is <3 percent and the difference is about 1 percent at 70 percent speed. This

level of agreement between the two analysis codes is consistent with previous results for another fan 

configuration as reported in Reference 7. The stall line shown is based on design intent and was not part 

of the current computations.

Overall, the HB results for pressure ratio show good correlation with the TURBO results, except for a

small shift in the speed line towards lower pressure ratio. This shift is slightly larger for 100 percent speed 

than for 70 percent speed. Also, for both speeds, the HB results show lower mass flow rate than the 

TURBO results for the same imposed back-pressure although this is not discernable from Figure 3.

Figure 4 is a plot of calculated corrected flow rate for various imposed back-pressure values that 

illustrates the differences in mass flow between TURBO and HB results. The differences between the 

calculated mass flow rates are not very large, except near the stall line (high back-pressure) at 100 percent 

speed. Further investigation is required to understand these differences in the results from the two 

analyses. For the unsteady computations and results to follow, it should be noted that a comparison for the 

same value of imposed back-pressure would be expected to show differences simply based on the 

different steady operating conditions that would result from the TURBO and HB analyses.

The efficiency computed from the HB and TURBO computations is plotted in Figure 5 and shows

excellent correlation between the two sets of results for both speeds, with some differences near the stall 

line. It should be noted that significant differences exist between the numerical method used in the 

TURBO computations and the HB computations. These differences include time-domain versus 

frequency-domain, numerical discretization, grid, and algorithm used to solve the RANS equations, 

turbulence modeling, and others. Therefore, the differences in results noted in Figures 3 to 5 are not 

entirely surprising.
NASA/TM—2010-216222 7

Figure 4.—Variation of mass flow rate with back-pressure 

from HB and time-domain codes.

Figure 5.—Predicted efficiency from HB and

time-domain codes.

(a) bending, 502 Hz (b) torsion, 1492 Hz

Figure 6.—Un-scaled mode shapes for experimental fan. Figure 7.—HB unsteady solution convergence for first 

vibration mode.

3.2 Unsteady Computational Results

Unsteady computations were performed using the HB code to predict the aeroelastic stability of the 

fan. These computations were carried out at both 100 and 70 percent speeds for the first and second 

structural vibration modes (bending and torsion) and several nodal diameter (ND) patterns (see Eq. (4)). 

The blade mode shapes are shown in Figure 6 as un-scaled contour plots of modal displacements at

100 percent speed. The mode shapes at 70 percent speed are substantially similar to those at 100 percent 

speed. HB and TURBO unsteady computations were carried out at various steady conditions for both 

speeds, which enabled trends along the speed line to be assessed.

3.2.1 100 Percent Speed, First Vibration Mode

At 100 percent speed, computations were carried out for the first blade vibration mode and several 

ND of the travelling wave covering the entire possible range. Figure 7 shows the convergence of the 

unsteady HB computations as a plot of non-dimensional complex generalized force with iteration counter, 

for various values of nodal diameter (–9 to 9). These HB computations were for an operating point near 

stall with a non-dimensional back-pressure of 1.11. As can be noted, although the rates of convergence 

vary somewhat with nodal diameter, all the results are very well converged. Similar convergence 

characteristics were observed for other operating conditions or back-pressure values.

Figure 8 shows the variation of the converged generalized force with nodal diameter for vibration in 

the first mode at a back-pressure value of 1.11 at 100 percent speed. For a single vibration mode, the 

stability is determined by the imaginary part of the generalized force and Figure 8 shows that the 

imaginary part of generalized force is negative for all nodal diameters, indicating flutter stability. Also, 

nodal diameters near zero have the lowest stability.
NASA/TM—2010-216222 8

Figure 8.—HB generalized force variation with nodal 

diameter for 100 percent speed, back-pressure 1.11, 

first mode.

Figure 9.—Variation of HB aerodynamic damping with 

nodal diameter for 100 percent speed, first mode.

Figure 10.—TURBO unsteady solution convergence 

for first vibration mode.

Figure 11.—Variation of TURBO aerodynamic 

damping with nodal diameter for 100 percent 

speed, first mode.

The generalized force values, shown in Figure 8, were used to calculate the aerodynamic damping, 

which is plotted in Figure 9. The results are presented for back-pressure values near of 1.11 (near stall) 

and 1.04 (away from stall). The aerodynamic damping, which is positive for all nodal diameters, is seen 

to vary by about 1.5 percent over the range, in contrast to previous results (Ref. 2) for a different 

configuration where the aerodynamic damping varied by an order of magnitude with nodal diameter. 

Also, for any given nodal diameter, the aerodynamic damping increases or decreases by only a small 

amount as back-pressure is changed, indicating only a weak dependence on operating condition. The 

lowest value of aerodynamic damping at a back-pressure of 1.11 is about 0.6 percent and occurs at 1ND.

For evaluating flutter stability at a fixed rotational speed, the trend in aerodynamic damping is 

relevant. Thus, in Figure 9, nodal diameters for which aerodynamic damping decreases with decreasing 

mass flow rate (increasing back-pressure) are particularly interesting. Nodal diameters of 3ND and 6ND 

are not of particular interest since the aerodynamic damping is increasing as we approach the stall line. 

Also, –3ND, –6ND, and ±9ND are not of particular interest because of the relatively large value of 

aerodynamic damping and relatively small variation with back-pressure. That leaves the low nodal 

diameters, which show low aerodynamic damping. But, they do not show much variation with backpressure in this case.

The time-domain TURBO code was used for similar unsteady computations for 100 percent speed 

and first structural vibration mode. Figures 10 and 11 show the corresponding results. Figure 10 shows 

the convergence of the unsteady computations in the form of aerodynamic damping calculated at the end 

of each vibration cycle. For the various nodal diameters considered, the unsteady results are well 

converged. Figure 11 shows the variation of aerodynamic damping for various nodal diameters for two 

values of back-pressure: 1.11 and 1.13, which are both near the stall line. The lowest value of 

aerodynamic damping at a back-pressure of 1.13 is about 0.6 percent at 1ND. The range of variation with
NASA/TM—2010-216222 9

Figure 12.—Variation of aerodynamic damping with 

mass flow rate for 100 percent speed, first mode.

nodal diameter is about 2 percent. The trend with increasing back-pressure is one of increasing 

aerodynamic damping (more stable) for large nodal diameters >3 and <–3. For low nodal diameters, there 

is only a small variation with back-pressure, which is the same trend noted in the HB results. 

The TURBO results in Figure 11 are very similar to the HB results in Figure 9. Recall that the 

operating points (mass flow rate) calculated from HB and TURBO are not the same for a given backpressure, as shown in Figure 4. Hence, the aerodynamic damping from Figures 9 and 11 is not expected to 

match for the same values of back-pressure. Nonetheless, the observed trends are very similar. One 

notable difference is seen in the comparison of the trends at large positive and large negative nodal 

diameters. The TURBO results show the same trend, increasing stability with increasing back-pressure, at 

both positive and negative high nodal diameters. However, the HB results show a smaller variation with 

back-pressure and the trends for large positive nodal diameters are opposite to those for large negative 

nodal diameters.

Unsteady computations were performed at additional back-pressure values for 1ND and 2ND 

traveling wave patterns. The resulting variation of aerodynamic damping with mass flow rate, computed 

from both the HB and TURBO codes, is shown in Figure 12. The results correlate very well between the 

two analyses in terms of the calculated values of damping, the relative values of damping for the 1ND and 

2ND traveling wave patterns, and the relatively small change with mass flow rate. The only notable 

difference is that the TURBO results show some weak trend with increasing and decreasing mass flow 

rate, whereas the HB results are largely unchanged for the mass flow rates considered. Based on the levels 

of aerodynamic damping calculated from HB and TURBO analyses, and the absence of any significant 

decrease along the speedline, no flutter is expected in the first mode at design speed.

3.2.2 100 Percent Speed, Second Vibration Mode

Unsteady computations were carried out for the second vibration mode in a manner similar to that 

described in the previous sub-section for the first vibration mode. Selected results are presented in this 

section. Figures 13 and 14 show the HB and TURBO variations of aerodynamic damping with nodal 

diameter for the two values of back-pressure used for the first mode calculations. A comparison between 

the two sets of results shows the similarities in the variation with nodal diameter. The amount of variation 

is significantly smaller than that obtained for the first mode. The lowest values of aerodynamic damping 

are smaller than those for the first mode. Also, Figures 13 and 14 show a similar trend with changing 

back-pressure. Results for 3ND and 6ND show significant variation with backpressure and a decrease 

with mass flow rate towards the stall line. Additional computations (with both HB and TURBO codes) are 

performed for these nodal diameter patterns and the results are presented in Figure 15 as variations of 

aerodynamic damping with mass flow rate. The results correlate very well between the two analyses in 

terms of the calculated values of damping and the relative values of damping for the 3ND and 6ND 

traveling wave patterns. Some difference is noted in terms of the slope of the variations with mass flow 

rate. The TURBO results show a clear linear trend, whereas the HB results show a smaller slope for 6ND 

and very little variation for 3ND.
NASA/TM—2010-216222 10

Figure 13.—Variation of HB aerodynamic damping 

with nodal diameter for 100 percent speed, second

mode.

Figure 14.—Variation of TURBO aerodynamic 

damping with nodal diameter for 100 percent

speed, second mode.

Figure 15.—Variation of aerodynamic damping with 

mass flow rate for 100 percent speed, second 

mode.

Since the aerodynamic damping drops slightly along the speedline, it is possible to extrapolate to a 

condition where the aerodynamic damping goes to zero. However, for the present configuration, it is clear 

that such calculated flutter points would fall considerably beyond the stall line (64 to 68 lbm/s depending 

on pressure ratio) shown in Figure 3. Thus, no flutter is expected in the second mode at design speed 

within the normal operating region. The computational results for 100 percent speed correlate well with 

the experimental observation of no flutter during wind-tunnel testing.

3.2.3 70 Percent Speed

Unsteady computations were carried out at 70 percent speed for the first and second vibration modes. 

The blade geometry and computational grids used were the same as those used for the 100 percent 

analyses. Selected results are presented here. Figure 16 shows the variation of HB aerodynamic damping 

with nodal diameter for the first vibration mode at 100 and 70 percent speeds. The same back-pressure 

value was prescribed for both sets of computations. The trend in variation with nodal diameter is 

essentially the same for both speeds; the lowest aerodynamic damping value changes from 0.6 percent at 

100 percent speed to 0.4 percent at 70 percent speed. Figure 17 shows the corresponding results for the 

second vibration mode, for which some change in the variation with nodal diameter is noted, with a peak 

occurring near –3ND for 70 percent speed.

Figure 18 shows the variation of aerodynamic damping with back-pressure for the second mode. In 

terms of the variation along the speedline (with back-pressure), the notable trend is the large drop at 6ND 

and 9ND when back-pressure is increased from 1.04 to 1.11. Based on the low aerodynamic damping at 

9ND, additional computations were done for this nodal diameter at 70 percent speed. The results are 

presented in Figure 19. Both the HB and TURBO show very similar values and trends, with a notable 

difference in the slope.
NASA/TM—2010-216222 11

Figure 16.—Variation of HB aerodynamic damping 

with nodal diameter at 100 and 70 percent speeds, 

first mode.

Figure 17.—Variation of HB aerodynamic damping with 

nodal diameter at 100 and 70 percent speeds, second 

mode.

Figure 18.—Variation of HB aerodynamic damping 

with nodal diameter and back-pressure; 70 percent

speed, second mode.

Figure 19.—Variation of aerodynamic damping with 

mass flow rate for 70 percent speed; second mode, 

9ND pattern.

Figure 20.—Fan map showing flutter points.

Based on the nearly linear variation, it is possible to extrapolate to zero aerodynamic damping; the 

mass flow rate at the calculated flutter point is approximately 44.4 lbm/s based on HB results, and 

approximately 34.5 lbm/s based on TURBO results. Notably, both these fall beyond the stall line

(46 to 47 lbm/s depending on pressure ratio), although the TURBO results indicate a far larger margin 

than indicated by the HB results, as shown in Figure 20. Note that this trend in flutter results from the two 

analysis codes is consistent with previous part-speed results for another fan configuration as reported in 

(Ref. 7), namely, that the HB result points to a flutter point at a higher mass flow rate as compared to 

TURBO. Further work is needed to understand the differences in the rate of change of aerodynamic 

damping along the speedline from the two analyses. Note that both results are consistent with the 

observation of no flutter during wind-tunnel testing.
NASA/TM—2010-216222 12

4.0 Conclusions

Steady and unsteady computations have been carried out for an experimental fan configuration using 

two turbomachinery aeroelastic analysis codes that solve the RANS equations with blade vibrations using 

a time-marching method (TURBO) and a HB. The computational results were summarized and compared. 

Overall, the TURBO and HB results are in good agreement. The steady results were compared on the 

performance map and correlate well. The steady HB results show a slightly lower pressure ratio than the 

steady TURBO results for the same mass flow rate, which is consistent with previous results for a 

different configuration. The difference in calculated pressure ratio is higher for the design speed than for 

the selected part-speed. For the same imposed back-pressure, the HB computations resulted in lower mass 

flow rates. The variation of efficiency with mass flow rate was compared and found to agree well with 

small difference near stall (HB calculations yield slightly higher values near stall). A detailed look at the 

flowfields predicted by HB and TURBO may provide additional information regarding the source of these 

differences. 

Unsteady computations have been performed at design speed for the first two structural vibration 

modes for various nodal diameter patterns. Both HB and TURBO computational results show very similar 

levels and trends of aerodynamic damping variations with nodal diameter. The range of variation with 

nodal diameter is less than 2 percent, in contrast to results previously obtained for a different fan. The 

second mode is seen to have a slightly lower aerodynamic damping than the first mode. The variations of 

aerodynamic damping with mass flow rate (along a speedline) are examined for the least stable nodal 

diameter patterns. The levels and trends correlate quite well between the two analyses. For the first 

vibration mode, both analyses show no significant decrease in aerodynamic damping along the speedline, 

leading to the conclusion that no flutter is expected, which correlates with the experimental observation of 

flutter stability. For the second mode, both analyses show that the aerodynamic damping for the +6ND 

nodal diameter pattern decreases slowly towards stall. However, with an aerodynamic damping value 

more than 0.5 percent, no flutter is expected. The computational results correlate well with the 

experimental observation of no flutter during wind-tunnel testing.

For 70 percent speed, the results for the first vibration mode are similar to those for 100 percent speed 

except for lower levels of aerodynamic damping. For the second mode, the 9ND nodal diameter pattern is 

seen to have a low level of aerodynamic damping along with a large rate of decrease in it towards the stall 

line. Based on the TURBO results, the extrapolated flutter point is well beyond the stall line, indicating 

no flutter in the normal operating region. The HB results also show that the extrapolated flutter point is 

beyond the stall line but very close to it in contrast to the TURBO results. Note that these results are 

consistent with the observation of no flutter during wind-tunnel testing. Further work is needed to 

understand the differences in the rate of change of aerodynamic damping along the speedline calculated 

by the two analyses, which resulted in different extrapolated flutter points.

In the present study, very small amplitude vibrations were prescribed to calculate linearized unsteady 

aerodynamics for a conventional linear flutter analysis as a first step to understand the characteristics of 

the HB code and of the flutter of an experimental fan. Future work will investigate the non-linear 

amplitude-dependent effects using the same HB aeroelastic code and the time-domain TURBO codes.

References

1. Marshall, J.G., and M. Imregun, “A review of aeroelasticity methods with emphasis on 

turbomachinery applications,” Journal of Fluids and Structures, Vol. 10, 1996, pp. 237–267.

2. Bakhle, M.A., Srivastava, R., Panovsky, J., Keith, T.G., Jr., Stefko, G.L., “Flutter Calculations of an 

Experimental Fan,” Paper US-2, CEAS/AIAA/NVvL International Forum on Aeroelasticity and 

Structural Dynamics, Amsterdam, June 4–6, 2003. 

3. Srivastava, R., Bakhle, M.A., and Keith, T.G., Jr., “Numerical Simulation of Aerodynamic Damping 

for Flutter Analysis of Turbomachinery Blade Rows,” Journal of Propulsion and Power, Vol. 19,

No. 2, March–April 2003. 
NASA/TM—2010-216222 13

4. Hall, K.C., Thomas, J.P., and Clark, W.S., “Computation of Unsteady Nonlinear Flows in Cascades 

Using a Harmonic Balance Technique,” AIAA Journal, Vol. 40, No. 5, May 2002.

5. Ekici, K., Hall, K.C., “Nonlinear Analysis of Unsteady Flows in Multistage Turbomachines Using 

Harmonic Balance,” AIAA Journal, Vol. 45, No. 5, May 2007, pp. 1047–1057. 

6. Fite, E.B., Woodward, R.P., and Podboy, G.G., “Effect of Trailing Edge Flow Injection on Fan Noise 

and Aerodynamic Performance,” AIAA Paper 2006–2844, June 2006. 

7. Bakhle, M.A., Thomas, J.P., and Reddy, T.S.R., “Fan Flutter Computations Using the Harmonic 

Balance Method,” AIAA Paper 2008–4743, July 2008.

8. Carta, F.O., “Coupled Blade-Disk-Shroud Flutter Instabilities in Turbojet Engine Rotors,” Journal of 

Engineering for Power, July 1967, pp. 419–426.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-08-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Harmonic Balance Computations of Fan Aeroelastic Stability 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Bakhle, Milind, A.; Reddy, T.S.R. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 561581.02.08.03.21.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17203 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216222 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 07 and 05 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A harmonic balance (HB) aeroelastic analysis, which has been recently developed, was used to determine the aeroelastic stability (flutter) 

characteristics of an experimental fan. To assess the numerical accuracy of this HB aeroelastic analysis, a time-domain aeroelastic analysis 

was also used to determine the aeroelastic stability characteristics of the same fan. Both of these three-dimensional analysis codes model the 

unsteady flowfield due to blade vibrations using the Reynolds-averaged Navier-Stokes (RANS) equations. In the HB analysis, the unsteady 

flow equations are converted to a HB form and solved using a pseudo-time marching method. In the time-domain analysis, the unsteady 

flow equations are solved using an implicit time-marching approach. Steady and unsteady computations for two vibration modes were 

carried out at two rotational speeds: 100 percent (design) and 70 percent (part-speed). The steady and unsteady results obtained from the two 

analysis methods compare well, thus verifying the recently developed HB aeroelastic analysis. Based on the results, the experimental fan 

was found to have no aeroelastic instability (flutter) at the conditions examined in this study.

15. SUBJECT TERMS

Fan blades; Flutter; Reynolds-averaged Navier-Stokes equation; Turbomachinery; Aeroelasticity; Aeroelastic stability 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

19 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18