0 Introduction

Flutter is a typical dynamic aeroelastic phenomenon, which could cause catastrophic results. Transonic flutter prediction is greatly important for the design of transport aircrafts and fighters. Nowadays, Computational Fluid Dynamics (CFD) has been widely applied in transonic flutter predictions^[1]. However, it is still very time-consuming to calculate instantaneous response of aeroelastic system by time-marching methods.

However, it is still very time-consuming to calculate instantaneous response of aeroelastic system by time-marching methods.

To improve analysis efficiency, reduced order model has been widely applied in flutter analysis.Weiwei Zhang et al^[2] applied ARX reduced order model to construct aeroelastic equation in state space form, which shows similar results with CFD/CSD simulation, except for transonic condition. Walter A. Silva^[3] created multiple-input multiple-output reduced order aeroelastic system based on Volterra series, which established the relationship between modal inputs and generalized aerodynamic forces (GAFs). However, the precision of unsteady aerodynamic forces based on Volterra series is very sensitive to amplitude of input signals and time steps^[4]. Weiwei Zhang et al^[5] also predicted flutter boundary by nonlinear ROM based on RBF neural network model.

Except for the time domain method, frequency domain method is another important approach of flutter analysis^[6]. Aerodynamic describing functions can be obtained by exciting aeroelastic system with harmonic disturbance of small amplitudes. Then, flutter velocity and frequency of different structural parameters can be obtained by V-g or P-K method, which shows great efficiency, comparing with time-marching methods^[7]. However, it takes excessive computational costs to calculate harmonic responses of aeroelastic system, because it needs to carry on unsteady flow simulation for each structural mode every harmonic frequency.

In addition, harmonic balance (HB) method shows superior efficiency in the simulation of periodic unsteady flow^[8-10], which is widely applied in the simulation of turbomachines^[11-12], helicopter rotors^[13-14], dynamic derivations^[15], et al. Therefore, a frequency domain flutter analysis method based harmonic balance is proposed in this paper, which takes efficiency and precision into account.

1 Time-marching method

Reynolds-Averaged Navier-Stokes (RANS) equation is the highest-fidelity flow solver in current engineering applications, which can be written as:

(1)

where Q and R are the vectors of flow variables and residuals respectively. Third-order Roe scheme is applied in spatial discretization, while approximate factorization scheme is used in implicit time-marching method. And S-A turbulence model is adopted.

Two-degree-of-freedom structural dynamic equation can be written as:

(2)

where

Eq.(2) can be transferred into state space form:

(3)

where

Eq.(3) can be solved by second-order Euler prediction-correction method:

(4)

The flow chart of time-marching approach of aeroelastic system is shown in Fig. 1.

2 Harmonic balance method

Harmonic balance method is an efficient approach of periodic unsteady flow field. Because of the periodicity, the vector of flow variables and residuals can be deployed by Fourier series:

(5)

The time derivation of Q is:

(6)

Substituting Eq.(5) into Eq.(1) and rearranging, we can get a pair of algebraic equations:

(7)

Eq.(7) can be written into matrix form:

(8)

where

Eq.(8) is the governing equation of original harmonic balance method. In this method, the calculation of residual terms by Fourier transformation is very complex and tedious, which is as follows:

(9)

To resolve this problem, sub-time levels were introduced by Hall et al^[11] to convert Q back to time domain through Fourier transformation:

(10)

where

{t_i, i=1, 2, ..., 2N_H} are equally distributed sub-time levels in a cycle, and t_i=2π(i-1)/2N_H.

Substituting Eq.(10) into Eq.(8):

(11)

Both sides in Eq.(11) are multiplied by E ^-1:

(12)

To achieve an easy approach, the residual terms are simplified by:

(13)

where

In this step, nonlinear residual terms can produce terms whose frequencies are higher that N_Hω, which are included in low-frequency terms by Eq.(13). This rough simplification may produce non-physical solutions in strong nonlinear problems^[16]. Therefore, the number of harmonics should be chose adequately to guarantee that high-frequency harmonic terms approach zero. Finally the governing equation of harmonic balance method can be obtained:

(14)

In order to solve Eq.(14) by time-marching method, pseudo time derivation is added:

(15)

Eq.(15) is the same as original RANS Eq.(1), except for its middle term (which is also called harmonic balance source term). Therefore, harmonic balance solver can be obtained from original steady RANS solver.

Then the harmonic balance method is validated by NACA64A010 airfoil under forced harmonic motion^[17]. The angle of attack changes as follows:

(16)

The flow conditions of two test cases are shown in Table 1. The comparison of unsteady aerodynamic forces between harmonic balance and time-marching methods are illustrated in Fig. 1 and Fig. 3. It is demonstrated that the unsteady forces calculated by harmonic balance and time-marching methods coincide with each other, which validates the precision of harmonic balance method. The time costs of time-marching and harmonic balance methods are shown in Table 2.

3 Frequency domain solver

When the aeroelastic system is undergoing single harmonic motion, the plunge and pitch displacements can be written as:

(17)

The corresponding lift and pitch moment coefficients can be written as:

(18)

Therefore, the structural Eq.(2) can be transferred into frequency domain:

(19)

Through aerodynamic describing function matrix F, the aerodynamic forces can be easily obtained by:

(20)

where

The aerodynamic describing functions can be obtained through exciting aeroelastic system by harmonic disturbance of small amplitude. Harmonic Balance method is applied to acquired aerodynamic describing functions efficiently.

However, the aeroelastic systemmay exhibit different behaviors (convergent oscillation, simple harmonic or divergent oscillation) in different freestream conditions. Therefore, an artificial damping coefficient j g K can be added into the aeroelastic system (19) to maintain harmonic motion:

(21)

The stability of aeroelastic system can be judged by the sign of artificial damping coefficient. Eq.(21) can be rearranged into a generalized complex eigenvalue problem:

(22)

where

At a fixed frequency ω, the matrix A will only change with freestream velocity. Different eigenvalue λ can be obtained in various velocity. When the real part of eigenvalue λ_Re is equal to 1.0, the imaginary part of eigenvalue λ_Im is the damping coefficient at a certain velocity. Thus we can get the relationship between damping coefficients and frequencies (or velocities). When damping coefficient is equal to 0, the frequency and velocity is flutter frequency and velocity respectively.

The procedure of frequency domain flutter analysis is as follows:

1) Calculate aerodynamic forces when the system are undergoing harmonic motionsat different frequencies;

2) Obtain aerodynamic transfer function matrix;

3) Construct generalized complex eigenvalue problem by aerodynamic transfer function matrix and structural dynamic equation;

4) Get the freestream velocity which makes the real part of eigenvalue λ_Re=1;

5) Obtain the relationship between damping coefficients and frequencies/velocities;

6) Acquire the flutter velocity and frequency.

4 Flutter analysis of Isogai wing

The Isogai Case A wing^[18-19] is a standard configuration of flutter analysis. The structural parameters of this case are: x_α=1.8, r_α²=3.48, a=-2.0, ω_h=ω_α=100rad/s, μ=60. The inviscid flutter velocity boundary shows a "S" shape in transonic Mach number, while the viscous flutter velocity boundary does not.

4.1 Aerodynamic describing function matrix

The aerodynamic describing function matrix is calculated byharmonic balance method to improve efficiency. Take h -/b=0.001 and ω=0.10 as an example. Tables 3 and 4 are the lift describing functions ( C_{L h/b}= C_LRe+j C _LIm ) and pitch moment describing functions ( C_{mḣ /b}= C_mRe+j C _mIm ) at two Mach number (Ma=0.60 and Ma=0.87) respectively.

It is illustrated that first-order harmonic is enough for subsonic case and more harmonics are necessary for transonic case to ensure the precision of aerodynamic describing functions.

4.2 Flutter boundary

The inviscid and viscous flutter boundaries are shown in Figs. 4 and 5. It is demonstrated that the difference between time-marching and frequency domain methods is insignificant, even in trmansonic range. The inviscid damping coefficient at Ma=0.88 is exhibited in Fig. 6, from which we can obtain the reason of "S" shape flutter boundary. The damping coefficient of plunge mode exceeds zero at a low flow velocity. But as the velocity increases, the damping coefficient returns back and falls below zero at a moderate velocity. Then the damping coefficient of pitch mode exceeds zero at a high velocity.

4.3 Influence of structural parameters

If the structural parameters undergo some change, flutter boundary can be easily obtained by generalized eigenvalue analysis in frequency domain method, other than time-marching method, which must carry on numerous CFD/CSD simulations. Therefore, the influences of mass ratio and frequency ratio on viscous flutter boundary are investigated below, based on the frequency domain method.

The influences of mass ratio μ and frequency ratio ω_h /ω_α are shown in Fig. 6 and 8 respectively. It is illustrated that higher mass ratio is preferable in subsonic region, while lower mass ratio is better in transonic region of Isogai A case. The frequency ratio also has a great impact on flutter boundary. When the frequency ratio is close to 1.0, the flutter velocity decreases significantly.

5 Conclusion

1) Comparing with time-marching method, harmonic balance method is an efficient approach of periodic unsteady flow;

2) A frequency domain flutter analysis method based on harmonic balance is proposed, which shows high efficiency and precision;

3) When the structural parameters undergo some change, new flutter boundary can be obtained quickly by generalized eigenvalue analysis in the proposed frequency domain method, other than time-marching method, which must carry on a series of CFD/CSD coupled simulations.