International Journal of Rotating Machinery
Volume 2016 (2016), Article ID 4293414, 12 pages
http://dx.doi.org/10.1155/2016/4293414

Research Article

Comparison of Frequency Domain and Time-Domain Methods for Aeromechanical Analysis

M. T. Rahmati

Department of Engineering Science, Oxford University, Parks Road, Oxford OX1 3PJ, UK

Received 19 December 2015; Revised 9 March 2016; Accepted 7 April 2016

Academic Editor: Gerard Bois

Copyright © 2016 M. T. Rahmati. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Unsteady flow around an oscillating plate cascade and that through a single compressor rotor subject to vibration have been computationally studied, aimed at examining the predictive ability of two low fidelity frequency methods compared with a high fidelity time-domain solution method for aeroelasticity. The computational solutions demonstrate the capabilities of the frequency domain methods compared with the nonlinear time-domain solution method in capturing small perturbations in the unsteady flow. They also show the great advantage of significant CPU time saving by the frequency methods over the nonlinear time method. Comparisons of two different frequency methods, nonlinear harmonic and phase solution method, show that these methods can produce different results due to the differences in numeric and physical conditioning. The results obtained using phase solutions method are in better agreement with the nonlinear time-domain solution. This is because the same numeric and physical conditioning are used in both the nonlinear time-domain method and phase solution frequency domain method.

1. Introduction

The blading aerodynamic design has long been employing steady flow methods as they are highly efficient and robust [1–3]. The steady solution methods can be simply automated for design applications using optimization technique or inverse approaches [4, 5]. Even at a detailed design stage, the solution to the steady flow equations rather than the unsteady flow equations is sought widely for blading aerodynamics predictions. This is because in many aerodynamics design applications the unsteadiness in the flow field particularly at its design operating point is usually small. As a consequence, the time-averaged flow solution is not greatly influenced by the unsteadiness perturbations. In fact, the current aerodynamic blade and airfoil designs have achieved high aerodynamic performance by using the steady flow model. The solution of unsteady flow equations is much more costly, requiring substantial computer resources [6–8].

There are many cases in which the unsteady forces have significant effects that cannot be ignored. For example, the effects on blade vibration, noise generation, fatigue, or failures of blades are very important. As modern multirow turbomachinery components are typically designed with high loading and more compact structural configurations, those unsteady effects are intensified in them. As aeroloading is increased, the blade mechanical integrity, aerodynamic noise, and vibration stress levels need to be carefully examined as the blades will be more vulnerable to flow-induced vibration problems. Recently a great advancement has been made in computing power and numerical methods for unsteady flows. However, within a multidisciplinary design environment there is always a need for developing efficient and fast prediction methods. This is because it is impractical to perform an aeromechanic related design optimization based upon costly unsteady flow solution methods, considering that a design optimization is carried out in an iterative process. As a result, significant research has been devoted in recent years to developing efficient numerical approaches capable of capturing the major unsteady features which are relevant to the engineering problems of interest while reducing the solution time to an acceptable level for use in routine design. One of the earliest methods is the time-linearized harmonic frequency domain methods which have been widely used for turbomachinery aeromechanical applications [9, 10]. In these methods the unsteady flow equation is considered as one steady equation and one perturbation equation. The unsteady perturbation of the flow solution is represented in a Fourier series. The validity of these methods is limited by the linear assumption which often exhibit solution divergence behaviours for highly nonlinear flows.

Recently advanced frequency domain such as nonlinear harmonic and phase solution methodologies have been developed which on the one hand is efficient for an unsteady flow solution for aeromechanics and on the other hand can be used for blading design optimization without much extra effort. One recent method, phase solution method [11–13], provides a straightforward simple method of modeling unsteady perturbations. In these methods the unsteady flow equations with a single periodic unsteadiness are solved at two or three distinctive phases of a period of unsteadiness. Using this approach the same computational efficiency as a conventional time-linearized method can be obtained. While the other time-linearized harmonic methods express the whole flow solution in a Fourier series, this method is based on casting the unsteady flow equations into a set of steady-like equations at a series of phases of a period of unsteadiness. By using this method the same steady flow solution method used for aerodynamic design can be used for aeromechanical design. Therefore the aeroelasticity of a blade can be optimized in a design optimization process to check both the aerodynamics and aeromechanics simultaneously. Recently, Valero et al. [14] developed a concurrent design optimization method on a phase solution method developed by Rahmati et al. [15] in conjunction with a commercial FEA solver for blade structural dynamics.

Validations and verifications are crucial in establishing our confidence in the numerical methods and thus the conclusions drawn upon those numerical results. It is not the aim of this study to apply the frequency domain method for very complex configurations. The application of various frequency domain methods for complex multirow turbomachines has been already reported [16, 17]. This paper focuses on two relatively simple cases, flow around a 3D flat plate oscillating cascade and through a compressor rotor subject to vibration, to demonstrate the validity and compare the effectiveness of three computational methods. The aim is to highlight the importance of the numerical diffusions and physical conditioning on the computational results. Although in the phase solutions method the same numeric and physical conditioning as the nonlinear time-domain method are used, this is not the case in the nonlinear harmonic method.

A brief description of the governing equations and the time and frequency domain methods is given in the next section; this is followed by the prediction of flow over an oscillating plate and flow through a single compressor rotor using various methods. These cases are provided to verify the capability of various numerical methods based on frequency domain method in predicting unsteady flow in comparison with time-domain. It is also used to show that it is important to develop both frequency and time-domain solution methods using the same numeric and physical conditioning. By doing so, when the high fidelity nonlinear time-domain solution is compared with the low fidelity frequency domain, one can consistently find the difference in the modeling fidelity without interference from differences in numeric.

2. Flow Modeling Formulation

2.1. Flow Governing Equations

The unsteady Reynolds averaged Navier-Stokes equations in a cylindrical coordinate system integral form can be described by the following form:

, , and are the convective flux vectors and , , and are the viscous flux vectors while is the source term due to rotating effects.

An absolute frame of reference is used for the solution of the flow equations. The standard one-equation Spalart-Allmaras model is implemented to calculate the eddy viscosity for the turbulence model. The equations are discretized in the finite-volume form with a cell-centred variable storage. The structured -mesh is used in all calculations.

The following form can be used to describe the semidiscrete integral form of the equations:

In the above equation the summation is over the boundary faces of a mesh cell. The normal projected areas in the three coordinate directions are , , and , respectively. In this formulation , , and are the flux terms, including the moving grid terms, the nonlinear convective fluxes, and viscous terms in the momentum and energy equations. The right-hand side of (2) can be considered together as a residual term. So, the equation can be simply expressed as

A 2nd-order cell-centre based finite-volume scheme with a blend of 2nd- and 4th-order numerical dissipations is used for spatial discretization (see Jameson et al. [18]). The time-marching solution is based on a 4-step Runge-Kutta integration in the pseudo time. The temporal change of flow variables at each fractional time step () is

The multigrid and local time stepping techniques are used to accelerate the time-marching solutions [15].

2.2. Time-Domain Direct Method for Unsteady Flow

The time-domain solution is obtained using the basic time-marching solver. It is important to note that a time-consistent multigrid is used for accelerating time-integration for an unsteady flow. A simple implementation of the time-consistent multigrid is to use just two grids. In this implementation, the spatial resolution is governed by the fine mesh while the coarse block is for enlarging the time step in a time-consistent manner [19]. Several levels of intermediate-meshes between the basic fine mesh and the coarse mesh are introduced. If is the levels of intermediate-meshes, the general multigrid formulation for the temporal change of flow variables for a fixed mesh can be written as

The subscripts “” denote the fine mesh while “” represent the intermediate-mesh (level ) and the coarse mesh is represented by “”. is the net flux; is the cell or block volume. and are the allowable time step lengths on the fine mesh and the intermediate-mesh of level (), respectively.

In order to maintain a uniform global time step length for unsteady calculations, the time step length for the coarse block needs to satisfy the time-consistence condition:

2.3. Harmonic Method

For simplicity, the flow governing equation for unsteady flow is rewritten as

Since the grid movement is usually prescribed and these terms can be calculated prior to the main flow field solution, in the above equation, the unsteady terms corresponding to the temporal change of the mesh volume are moved to the right-hand side. In a frequency domain method the unsteady flow is assumed to be composed of a steady and a fluctuating part:

where is the steady part of the flow variables while is the fluctuating part. The fluctuating part is harmonic in time for many flows of engineering interest. So the unsteady perturbation flow can be expressed as a Fourier series. Since a linear assumption has been made, the behaviour of each Fourier component can be analysed individually and then summed together to form the total solution. Therefore, for a single periodic disturbance, the fluctuating part can be represented as

where is the vector of complex amplitudes of perturbations while is the complex conjugate of .

Substituting the relationships (8) into (9) yields the following equation:

Substituting the relationships (10) into (7) and collecting the zeroth- and first-order terms lead to two equations: one an equation for the steady and the other one to a linearized perturbation equation. The steady equation is solved by ignoring the time-dependent term and has the following form:

A similar equation can be found for the linearized perturbation equation [20]. So basically in this method, first, a steady flow solution is obtained by solving the steady Navier-Stokes equation. Then, for a given frequency and interblade phase angles, the coefficients of the time-linearized equations are formed from the steady flow solution. Finally the time-linearized perturbation equations are solved. Thus, a time-dependent unsteady problem in the time-domain is effectively transformed to solving two equations in the frequency domain. This linear method can be extended to nonlinear method based on the work of Adamczyk [21] who showed that time-averaging the Navier-Stokes equations resulted in the inclusion of the effect of periodic perturbations on the mean flow through stress terms. The nonlinear harmonic method is similar to RANS equations in which the effect of turbulence is included by introducing the Reynolds stresses. Extra closure models are required to work out these deterministic stress terms similar to turbulence modeling for Reynolds stress terms. In this approach, a significant modification is that a time-averaged value instead of a steady flow value is used as the basis for the harmonic perturbations. So, it is assumed that the flow field is composed of a time-averaged flow with a small perturbation unsteady flow:

Here is the vector of the time-averaged conservative variables and is the vector of perturbation conservative variables. By substituting expressions (12) into (7) and taking a time-averaging, the resultant time-averaged Navier-Stokes equations are obtained. In comparison with the steady form used in linear methods the time-averaging produces extra terms because of the nonlinearity of the flow governing equations. The extra stress terms due to the velocity fluctuations are the result of nonlinearity of the flow governing equations. One assumption made here is that the unsteady perturbation is dominated by first-order terms so the first-order harmonic perturbation equation has the same form as the unsteady perturbation equation in the time-linearized method. Obviously if the time-averaged flow is the same as the steady flow, the above first harmonic perturbation equation reduces to the conventional time-linearized perturbation equation. The time-averaging procedure produces extra unsteady stress terms in the time-averaged equations which are evaluated from unsteady perturbations. While the unsteady perturbations are obtained from solving the harmonic perturbation equation, the coefficients of disturbance equations resulted from the solution of time-averaged equation and this interaction is achieved through a strong coupling procedure. More details regarding the governing equations can be found in [20].

2.4. Phase-Shift Solution Method

In this formulation similar to nonlinear harmonic method, by considering only one harmonic, the flow solution is given by

By substituting (13) into (7) the following equation is obtained:

At three different temporal phases (), (13) can be written as follows:

In (15a)–(15c) the three unknowns , , and are expressed in terms of the flow solutions at the three phases. Substituting these unknowns in terms of flow solutions at three phases into (7) yields the following equations:

These equations are very similar to the steady equations. The difference is that there are three sets of equations instead of only one and there is an extra source term in each equation. The equations are instantaneously solved in a similar way to that of the steady RANS equations. The source terms are simply calculated using the flow solution at the previous iteration. One advantage of this phase solution method in comparison with the conventional harmonic solution method or linear harmonic method is that the nonlinearity is automatically retained in the nonlinear convection terms in the discretized equations. This feature leads to much improved solution convergence for harmonic disturbances, particularly for cases with flow separation. More details regarding this method can be found in [16].

2.5. Boundary Conditions

Similar boundary conditions are used in all methods. The slip wall condition is applied and the log-law is applied on solid blade and end wall surfaces to determine the surface shear stress. For inlet and exit to the computational domain boundaries, the same reflective treatment as those for the steady flow solutions can be used by specifying inlet stagnation pressure, stagnation temperature, inlet flow angle, and exit static pressure. Alternatively, the local characteristic disturbance-based nonreflective approach by Giles [22] can be used for inlet and outlet boundary treatments. The wall and far field inlet/exit conditions are applied in exactly the same way for both the time-domain and frequency domain solutions.

For the periodic boundaries, the direct periodic (repeating) condition is applied for the time-domain method. Similarly, a direct periodic condition is applied to a steady flow solution. For the frequency domain solution method using a single passage domain the harmonic components are phase-shifted between the upper and lower periodic boundaries by a given interblade phase angle, (the phase leads to the upper blade vibration relative to its lower neighbouring one).

3. Flow over an Oscillating Plate

A three-dimensional flat plate cascade test case is used for the comparison of the frequency domain and time-domain method.

The geometric parameters and flow conditions are shown in Figure 1 and Table 1. The blades are oscillating in a three-dimensional mode. Each two-dimensional section is subject to torsion mode around its leading edge with the torsion amplitude varying linearly along span from at the lower end to at the upper end.

Table 1: Geometric parameters and flow conditions for the flat plate test case.

True chord length
0.1 m

Stagger angle
= 45°

Pitch/chord ratio
1

Span/chord ratio
3

Steady incidence
0

Steady Mach number
0.7

Reduced frequency
1

Interblade phase angle
0, 180

Figure 1: Three-dimensional flat plate cascade test case geometry.

The present calculations were carried out using both the frequency domain and the time-domain method with a mesh density of in the streamwise, pitchwise, and radial directions, respectively.

Calculated unsteady pressure jumps at each two-dimensional section are presented in the form of :

where is the first harmonic pressure jump across the blade; is the torsion amplitude at the tip in radian. is a complex number with its real and imaginary parts, which can be used to work out amplitude and phase angle. The calculations were performed at the spanwise position of various normalized distance measured from the hub section ().

Figures 2–5 show the unsteady pressure on the lower and upper plate. The phase solution analysis results are virtually identical to the time-domain analysis results. Similar results are obtained for the linear and nonlinear harmonic method. Figures 6 and 7 show the chordwise distributions of the real and imaginary parts of the unsteady pressure jump at six spanwise sections for the interblade phase angle of 180 and 0 deg. The numerical results are compared with the analytical results for this case provided by Namba [23, 24] and He and Denton [25]. The linear and nonlinear harmonic results are very similar. This is something within expectation considering that the same numeric is used in these two different fidelity modes and no extra unsteadiness is available in the models. However, compared to nonlinear time-domain method, the phase solution method gives better agreement in particular for interblade phase angle of 180 degrees. It is important to note that in the phase solution formulation no extra numerical method is used for calculations of flow perturbations. So, basically, both phase solution and time-domain solution methods use the same numeric and physical conditioning. So, a comparison between the phase solution method and the time-domain method shows the difference in the modeling fidelity, which is none here, without interference from differences in numeric.

(a) Upper plate

(b) Lower plate

Figure 2: Unsteady pressure (Pa), frequency method (interblade phase angle = 180 deg).

(a) Upper plate

(b) Lower plate

Figure 3: Unsteady pressure (Pa), time-domain method (interblade phase angle = 180 deg).

(a)

(b)

Figure 4: Real and imaginary parts of unsteady pressure, frequency domain method (interblade phase angle = 180 deg).

(a)

(b)

Figure 5: Real and imaginary parts of unsteady pressure, time-domain method (interblade phase angle = 180 deg).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Comparison of real and imaginary parts of unsteady pressure (interblade phase angle = 0 deg).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: Comparison of real and imaginary parts of unsteady pressure (interblade phase angle = 180 deg).

A very small discrepancy exists between the calculation and analytical data in all methods around the inlet due to reflection effects at inlet. A typical time-domain solver run for this case takes about 37 hours on a single CPU using a 2 GHz Intel Xeon CPU. The CPU time saved using the frequency domain solution method is quite significant as the solution in all methods takes about thirty to forty minutes using the same computer.

4. DLR Compressor Rotor

Three-dimensional unsteady flow through a transonic compressor from the German Aerospace Center (DLR, Koln, Germany) subject to blade vibration is computed to evaluate the capabilities of the frequency domain and time-domain methods. The rotor rotational speed is 20,260 RPM and it gives a stage pressure ratio of 1.59 at a mass flow rate of 17.15 kg/s. The compressor has 28 rotor blades and 31 stator blades. The detailed information of the compressor and test rig can be found in [3]. In previous studies [9, 11] multistage effects on the aerodamping of this DLR compressor have been investigated using both frequency domain and time-domain methods. In this study, unsteady flow through the rotor row has been computationally studied. The aim is to examine the predictive ability of two low fidelity frequency methods for aeroelasticity when the only source of unsteadiness is the rotor vibration. To make sure that no other source of unsteadiness (e.g., reflections from the downstream outlet) influences the unsteady flow through the compressor, the outlet in the computational domain is placed far away from the trailing edge. Figure 8 shows the 3D computational model of the rotor while Figure 9 shows the blade-to-blade view of the computational mesh at the mid-span section used in current computation. A total mesh point of 162500 is used in a single row.

Figure 8: The geometry of the DLR compressor rotor.

Figure 9: Computational mesh for the DLR compressor rotor.

In the modal analysis to get more realistic mode shape and fundamental vibration frequencies, the chosen material for the rotor blade has a density of 4428.8 kg/m³, Young’s modulus of 1251011 Pa, and Poisson ratio of 0.27. These properties are close to the property of titanium at the ambient temperature.

The unsteady flow through the rotor blade at various Nodal Diameters (ND) at one modal frequency is investigated. The ND of 4 with IBPA of 51.428° at mode 1 vibration frequency of 1350 Hz is examined in this paper. The nonlinear time-domain method is used for validations of two frequency domain methods: nonlinear harmonic and phase solution method. A typical time-domain solver run for this case takes about 3 days on a single CPU with a 2 GHz Intel Xeon CPU. The CPU time required for the frequency domain solution is only 2 hrs using the same computer.

Figure 10 shows the time-averaged pressure values at the mid-span of the blade calculated using time-domain and frequency domain method. Figure 11 shows the unsteady part of pressure values at the mid-span of the rotor. No noticeable difference can be observed in those pressure contours. Detailed comparisons of time-averaged pressure, 1st harmonic pressure amplitude, and phase angle are shown in Figures 12–14. Compared to the time-domain method both nonlinear harmonic and phase solution method predict the same time-averaged values for pressure at the mid-span. The comparisons of the 1st harmonic pressure amplitudes and phases show that the phase solution method predicts the same pressure amplitude and phase at the mid-span as those of the time-domain one. The results obtained using the nonlinear harmonic method, however, are slightly different. This method predicts lower values for the 1st harmonic pressure amplitudes at both suction side and pressure side and a higher value for phase at pressure side. As the only differences between these two frequency domain methods are the formulation and numeric used for predicting the unsteady perturbation, and there is no other source of unsteadiness in the flow, it can be concluded that these differences are because of the numerical diffusion in this model. The phase solution method predicts the same results as those of the time-domain method because in this method the same numeric and physical conditioning are used and no extra formulation is used to take into account the unsteady perturbation. Indeed, in this method, the same equations are solved at three phases. Since there is only one unsteady harmonic in the model, one expects that these two models predict the same results. This case highlights the importance of using the same numeric and physical conditioning for both low fidelity and high fidelity solution methods. When the high fidelity nonlinear time-domain solution is compared with the low fidelity phase solution method, one can consistently find the difference in the modeling fidelity without interference from differences in numeric. However, the same cannot be said for the nonlinear harmonic method due to the difference in numeric and the subsequent differences in the numerical diffusions compared with the time-domain method.

(a) Frequency domain

(b) Time-domain

Figure 10: Time-averaged pressure at mid-span.

(a) Frequency domain

(b) Time-domain

Figure 11: Amplitude of the 1st harmonic unsteady pressure at mid-span.

Figure 12: The time-averaged pressure coefficient at mid-span.

Figure 13: Phase angle of the 1st harmonic unsteady pressure at mid-span.

Figure 14: The 1st harmonic unsteady pressure coefficient at mid-span.

5. Conclusion

A linear oscillating plate cascade and a transonic compressor rotor are investigated computationally to assess the validity, merits, and limitations of two low fidelity frequency domains and a high fidelity time-domain solution method. The frequency domain methods can produce one or two orders of magnitude saving in computational effort. On practical relevance and implications, this low computational cost makes the frequency domain method a valuable tool in an aeroelasticity design environment in which a design optimization is carried out in an iterative process. Comparisons of the nonlinear harmonic and phase solution results demonstrated good agreement with the time-domain method over a wide range of spanwise positions. The results using phase solution method are virtually identical to the time-domain analysis results. The outcomes of the present work suggest that in aeroelasticity applications, in which the perturbation magnitudes are usually very small compared to the mean flow, the solution of a nonlinear phase solution method can be used for analysing unsteady flow. By comparing the time-domain solution with the phase solution method, one can consistently find the difference in the modeling fidelity without interference from differences in numeric. This is because in the phase solution method the same numeric and physical conditioning as the high fidelity time-domain are used.

Nomenclature

: | Torsion amplitude at the lower end
: | Torsion amplitude at the upper end
: | Stagger angle
: | First harmonic pressure jump, Pa
: | Complex pressure jump coefficient
: | Real part of pressure jump coefficient
: | Imaginary part of pressure jump coefficient
: | Vibration frequency, Hz
: | Radial direction, m
: | Spanwise position
: | Axial direction, m
: | Spanwise direction, m
: | Chordwise position
: | Phase angle, deg
: | Density, kg/m³
: | Steady flow component
: | Fluctuating flow component
: | Time-averaged flow component
: | Unsteady flow variable.

Disclosure

Present address of M. T. Rahmati is Department of Mechanical, Aerospace, and Civil Engineering, Brunel University London, Uxbridge UB8 3PH, UK.

Competing Interests

The author declares having no competing interests.

Acknowledgments

First and foremost thanks are due to Professor L. He from Oxford University for Providing Namba’s analytical data and valuable discussion and help throughout the course of this research at Oxford. The author would like to thank Siemens Industrial Turbomachinery Ltd., Lincoln, for providing funding for this research at Oxford University. Thanks are also due to Professor M. Namba from Sojo University, Japan, for valuable discussion regarding his analytical results for the three-dimensional flat plate cascade used in this study.

References

P. Hu, M. Zangeneh, B. Choo, and M. T. Rahmati, “On design of a transonic fan rotors by 3D inverse design method,” in Proceedings of the ASME Turbo Expo, GT2006-91173, pp. 447–458, Barcelona, Spain, May 2006.

M. T. Rahmati, “Application of a pressure correction method for modeling incompressible flow through turbomachines,” International Journal of Computational Methods, vol. 6, no. 3, pp. 399–411, 2009.

N. A. Cumpsty and E. M. Greitzer, “Ideas and methods of turbomachinery aerodynamics: a historical view,” Journal of Propulsion and Power, vol. 20, no. 1, pp. 15–26, 2004.

M. T. Rahmati, “An inverse approach for turbomachinery blade design,” AIAA Journal, vol. 47, pp. 703–709, 2009.

M. T. Rahmati, “A new Navier-Stokes inverse method based on mass-averaged tangential velocity for blade design,” International Journal for Numerical Methods in Fluids, vol. 60, no. 3, pp. 323–336, 2009.

J. G. Marshall and M. Imregun, “A review of aeroelasticity methods with emphasis on turbomachinery applications,” Journal of Fluids and Structures, vol. 10, no. 3, pp. 237–267, 1996.

M. Vahdati, A. I. Sayma, M. Imregun, and G. Simpson, “Multibladerow forced response modeling in axial-flow core compressors,” Journal of Turbomachinery, vol. 129, no. 2, pp. 412–422, 2007.

M. T. Rahmati, L. He, and D. Wang, “Non-linear time and frequency domain methods for multi-row aeromechanical analysis,” in Proceeding of the ASME Turbo Expo 2012: Power for Land, Sea and Air, GT2012-68723, Copenhagen, Denmark, June 2012.

H. D. Li and L. He, “Blade aerodynamic damping variation with rotor-stator gap: a computational study using single-passage approach,” Journal of Turbomachinery, vol. 127, no. 3, pp. 573–579, 2005.

K. C. Hall and C. B. Lorence, “Calculation of three-dimensional unsteady flows in turbomachinery using the linearized harmonic Euler equations,” Journal of Turbomachinery, vol. 115, no. 4, pp. 800–809, 1993.

L. He, “Harmonic solution of unsteady flow around blades with separation,” AIAA Journal, vol. 46, no. 6, pp. 1299–1307, 2008.

K. C. Hall, J. P. Thomas, and W. S. Clark, “Computation of unsteady nonlinear flows in cascades using a harmonic balance technique,” AIAA Journal, vol. 40, no. 5, pp. 879–886, 2002.

M. T. Rahmati, L. He, and R. G. Wells, “Interface treatment for harmonic solution in multi-row aeromechanic analysis,” in Proceedings of the ASME Turbo Expo 2010: Power for Land, Sea and Air, GT2010-23376, Glasgow, UK, June 2010.

O. Valero, L. He, and Y. S. Li, “Concurrent blade aerodynamic-aeroelastic design optimization with re-scaled response surface,” in Proceedings of the ASME Turbo Expo: Turbine Technical Conference and Exposition (GT '14), Paper no. GT2014-25854, pp. 1–15, Düsseldorf, Germany, June 2014.

M. T. Rahmati, L. He, D. X. Wang, Y. S. Li, R. G. Wells, and S. K. Krishnababu, “Nonlinear time and frequency domain methods for multirow aeromechanical analysis,” ASME Journal of Turbomachinery, vol. 136, no. 4, Article ID 041010, 2014.

M. T. Rahmati, L. He, and Y. S. Li, “The blade profile orientations effects on the aeromechanics of multirow turbomachines,” Journal of Engineering for Gas Turbines and Power, vol. 138, no. 6, Article ID 062606, 2015.

M. T. Rahmati, L. He, and Y. S. Li, “Multi-row interference effects on blade aeromechanics in compressor and turbine stages,” in Proceedings of the 13th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines (ISUAAAT '13), Tokyo, Japan, September 2012.

A. Jameson, W. Schmidt, and E. Turkel, “Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes,” in Proceedings of the 14th Fluid and Plasma Dynamics Conference, Palo Alto, Calif, USA, 1981.

L. He, “Three-dimensional unsteady navier-stokes analysis of stator-rotor interaction in axial-flow turbines,” Journal of Power and Energy, vol. 214, pp. 13–22, 2000.

W. Ning and L. He, “Computation of unsteady flows around oscillating blades using linear and nonlinear harmonic Euler methods,” Journal of Turbomachinery, vol. 120, no. 3, pp. 508–514, 1998.

J. J. Adamczyk, “Model equations for simulating flows m multistage turbomachinery,” ASME Paper 85-GT-226, 1985.

M. B. Giles, “Stator/rotor interaction in a transonic turbine,” AIAA Journal of Propulsion and Power, vol. 6, no. 5, pp. 621–627, 1990.

M. Namba, “Lifting surface theory for unsteady flows in rotating, annular cascade,” in Proceedings of the 1st IUTAM Symposium on Aero-Elasticity in Turbomachines, Paris, France, 1976.

M. Namba, “Three-dimensional flows,” in AGARD Manual on Aeroelasticity in Axial-Flow Turbomachines, M. F. Platzer and F. O. Carta, Eds., vol. 1 of Unsteady Turbomachinery Aerodynamics, chapter 4, AGARD-AG-298, 1987.

L. He and J. D. Denton, “Three-dimensional time-marching inviscid and viscous solutions for unsteady flows around vibrating blades,” Journal of Turbomachinery, vol. 116, no. 3, pp. 469–476, 1994.